Add viscosity test problem

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") {

src/pgen/CMakeLists.txt

@@ -24,4 +24,5 @@ target_sources(athenaPK PRIVATE
     rand_blast.cpp
     sod.cpp
     turbulence.cpp
+    visc.cpp
     )

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;
 

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