From 1c514353b1680926c1ce7b1dcd2406f220d9941b Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Thu, 28 Sep 2023 16:01:22 +0200
Subject: [PATCH] Add viscosity test problem

---
 src/main.cpp            |  2 ++
 src/pgen/CMakeLists.txt |  1 +
 src/pgen/pgen.hpp       |  6 +++++
 src/pgen/visc.cpp       | 56 +++++++++++++++++++++++++++++++++++++++++
 4 files changed, 65 insertions(+)
 create mode 100644 src/pgen/visc.cpp

diff --git a/src/main.cpp b/src/main.cpp
index 93f984bb..81bbc7bc 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -74,6 +74,8 @@ int main(int argc, char *argv[]) {
     pman.app_input->ProblemGenerator = advection::ProblemGenerator;
   } else if (problem == "orszag_tang") {
     pman.app_input->ProblemGenerator = orszag_tang::ProblemGenerator;
+  } else if (problem == "visc") {
+    pman.app_input->ProblemGenerator = visc::ProblemGenerator;
   } else if (problem == "diffusion") {
     pman.app_input->ProblemGenerator = diffusion::ProblemGenerator;
   } else if (problem == "field_loop") {
diff --git a/src/pgen/CMakeLists.txt b/src/pgen/CMakeLists.txt
index 382b5b7b..9e566580 100644
--- a/src/pgen/CMakeLists.txt
+++ b/src/pgen/CMakeLists.txt
@@ -24,4 +24,5 @@ target_sources(athenaPK PRIVATE
     rand_blast.cpp
     sod.cpp
     turbulence.cpp
+    visc.cpp
     )
diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp
index 182d8497..ef40620a 100644
--- a/src/pgen/pgen.hpp
+++ b/src/pgen/pgen.hpp
@@ -67,6 +67,12 @@ using namespace parthenon::driver::prelude;
 void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
 } // namespace orszag_tang
 
+namespace visc {
+using namespace parthenon::driver::prelude;
+
+void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
+} // namespace visc
+
 namespace diffusion {
 using namespace parthenon::driver::prelude;
 
diff --git a/src/pgen/visc.cpp b/src/pgen/visc.cpp
new file mode 100644
index 00000000..7144c341
--- /dev/null
+++ b/src/pgen/visc.cpp
@@ -0,0 +1,56 @@
+
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR MHD code
+// Copyright (c) 2023, Athena Parthenon Collaboration. All rights reserved.
+// Licensed under the 3-Clause License (the "LICENSE")
+//========================================================================================
+//! \file visc.cpp
+//! \brief Viscous diffusion of a 1D Gaussian
+//========================================================================================
+
+// Parthenon headers
+#include "mesh/mesh.hpp"
+#include <parthenon/driver.hpp>
+#include <parthenon/package.hpp>
+
+// AthenaPK headers
+#include "../main.hpp"
+
+namespace visc {
+using namespace parthenon::driver::prelude;
+
+void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
+  IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  auto &mbd = pmb->meshblock_data.Get();
+  auto &u = mbd->Get("cons").data;
+
+  Real gm1 = pin->GetReal("hydro", "gamma") - 1.0;
+  Real v1 = 0., v2 = 0., v3 = 0.;
+  Real d0 = 1.;
+  Real p0 = 1.;
+  auto amp = pin->GetOrAddReal("problem/visc", "amp", 1.e-6);
+  auto t0 = pin->GetOrAddReal("problem/visc", "t0", 0.5);
+  auto nu_iso = pin->GetReal("diffusion", "mom_diff_coeff_code");
+  auto x1_0 = pin->GetOrAddReal("problem/visc", "x1_0", 0.);
+
+  auto &coords = pmb->coords;
+
+  pmb->par_for(
+      "ProblemGenerator: viscosity", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int k, const int j, const int i) {
+        u(IDN, k, j, i) = d0;
+        u(IM1, k, j, i) = u(IDN, k, j, i) * v1;
+        u(IM2, k, j, i) =
+            u(IDN, k, j, i) * amp / std::pow(std::sqrt(4. * M_PI * nu_iso * t0), 1.0) *
+            std::exp(-(std::pow(coords.Xc<1>(i) - x1_0, 2.)) / (4. * nu_iso * t0));
+        u(IM3, k, j, i) = u(IDN, k, j, i) * v3;
+        u(IEN, k, j, i) =
+            p0 / gm1 +
+            0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) /
+                u(IDN, k, j, i);
+      });
+}
+} // namespace visc