Add linwave3d decay diffusion test and fix parabolic dt

inputs/diffusion.in

@@ -70,7 +70,7 @@ mom_diff_coeff_code = 0.25
 resistivity = isotropic     # none (disabled) or isotropic
 resistivity_coeff = fixed
 ohm_diff_coeff_code = 0.25
-rkl2_max_dt_ratio = 200.0
+rkl2_max_dt_ratio = 400.0
 
 <parthenon/output0>
 file_type = hdf5

src/hydro/diffusion/conduction.cpp

@@ -179,8 +179,8 @@ Real EstimateConductionTimestep(MeshData<Real> *md) {
         },
         Kokkos::Min<Real>(min_dt_cond));
   }
-
-  return fac * min_dt_cond;
+  const auto &cfl_diff = hydro_pkg->Param<Real>("cfl_diff");
+  return cfl_diff * fac * min_dt_cond;
 }
 
 //---------------------------------------------------------------------------------------

src/hydro/diffusion/resistivity.cpp

@@ -84,7 +84,8 @@ Real EstimateResistivityTimestep(MeshData<Real> *md) {
     PARTHENON_THROW("Needs impl.");
   }
 
-  return fac * min_dt_resist;
+  const auto &cfl_diff = hydro_pkg->Param<Real>("cfl_diff");
+  return cfl_diff * fac * min_dt_resist;
 }
 
 //---------------------------------------------------------------------------------------

src/hydro/diffusion/viscosity.cpp

@@ -88,7 +88,8 @@ Real EstimateViscosityTimestep(MeshData<Real> *md) {
     PARTHENON_THROW("Needs impl.");
   }
 
-  return fac * min_dt_visc;
+  const auto &cfl_diff = hydro_pkg->Param<Real>("cfl_diff");
+  return cfl_diff * fac * min_dt_visc;
 }
 
 //---------------------------------------------------------------------------------------

src/hydro/hydro.cpp

@@ -665,6 +665,10 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     }
     if (diffint != DiffInt::none) {
       pkg->AddParam<Real>("dt_diff", 0.0, true); // diffusive timestep constraint
+      // As in Athena++ a cfl safety factor is also applied to the theoretical limit.
+      // By default it is equal to the hyperbolic cfl.
+      auto cfl_diff = pin->GetOrAddReal("diffusion", "cfl", pkg->Param<Real>("cfl"));
+      pkg->AddParam<>("cfl_diff", cfl_diff);
     }
     pkg->AddParam<>("diffint", diffint);
 

src/main.cpp

@@ -55,6 +55,7 @@ int main(int argc, char *argv[]) {
     pman.app_input->InitUserMeshData = linear_wave_mhd::InitUserMeshData;
     pman.app_input->ProblemGenerator = linear_wave_mhd::ProblemGenerator;
     pman.app_input->UserWorkAfterLoop = linear_wave_mhd::UserWorkAfterLoop;
+    Hydro::ProblemInitPackageData = linear_wave_mhd::ProblemInitPackageData;
   } else if (problem == "cpaw") {
     pman.app_input->InitUserMeshData = cpaw::InitUserMeshData;
     pman.app_input->ProblemGenerator = cpaw::ProblemGenerator;

src/pgen/linear_wave_mhd.cpp

@@ -708,4 +708,38 @@ void Eigensystem(const Real d, const Real v1, const Real v2, const Real v3, cons
   left_eigenmatrix[6][6] = left_eigenmatrix[0][6];
 }
 
+// For decaying wave with diffusive processes test problem, dump max V_2
+Real HstMaxV2(MeshData<Real> *md) {
+  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
+
+  const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
+
+  IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
+  IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
+
+  Real max_v2 = 0.0;
+
+  Kokkos::parallel_reduce(
+      "HstMaxV2",
+      Kokkos::MDRangePolicy<Kokkos::Rank<4>>(
+          parthenon::DevExecSpace(), {0, kb.s, jb.s, ib.s},
+          {prim_pack.GetDim(5), kb.e + 1, jb.e + 1, ib.e + 1},
+          {1, 1, 1, ib.e + 1 - ib.s}),
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lmax) {
+        lmax = Kokkos::fmax(lmax, Kokkos::fabs(prim_pack(b, IV2, k, j, i)));
+      },
+      Kokkos::Max<Real>(max_v2));
+
+  return max_v2;
+}
+
+void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg) {
+  if (pin->GetOrAddBoolean("problem/linear_wave", "dump_max_v2", false)) {
+    auto hst_vars = pkg->Param<parthenon::HstVar_list>(parthenon::hist_param_key);
+    hst_vars.emplace_back(parthenon::HistoryOutputVar(
+        parthenon::UserHistoryOperation::max, HstMaxV2, "MaxAbsV2"));
+    pkg->UpdateParam(parthenon::hist_param_key, hst_vars);
+  }
+}
 } // namespace linear_wave_mhd

src/pgen/pgen.hpp

@@ -25,6 +25,7 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin);
 void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
 void UserWorkAfterLoop(Mesh *mesh, parthenon::ParameterInput *pin,
                        parthenon::SimTime &tm);
+void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg);
 } // namespace linear_wave_mhd
 
 namespace cpaw {

tst/regression/CMakeLists.txt

@@ -52,6 +52,9 @@ setup_test_both("aniso_therm_cond_gauss_conv" "--driver ${PROJECT_BINARY_DIR}/bi
 setup_test_both("diffusion" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \
   --driver_input ${PROJECT_SOURCE_DIR}/inputs/diffusion.in --num_steps 12" "convergence")
 
+  setup_test_both("diffusion_linwave3d" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \
+  --driver_input ${PROJECT_SOURCE_DIR}/inputs/linear_wave3d.in --num_steps 2" "convergence")
+
 setup_test_both("field_loop" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \
   --driver_input ${PROJECT_SOURCE_DIR}/inputs/field_loop.in --num_steps 12" "convergence")
 

tst/regression/test_suites/diffusion/diffusion.py

@@ -24,7 +24,7 @@
 int_cfgs = ["unsplit", "rkl2"]
 res_cfgs = [256, 512, 1024]
 tlim = 2.0
-nu = 0.25
+diff_coeff = 0.25
 
 all_cfgs = list(itertools.product(diff_cfgs, res_cfgs, int_cfgs))
 
@@ -79,9 +79,10 @@ def Prepare(self, parameters, step):
             f"diffusion/viscosity={viscosity_}",
             f"diffusion/resistivity={resistivity_}",
             # we can set both as their activity is controlled separately
-            f"diffusion/mom_diff_coeff_code={nu}",
-            f"diffusion/ohm_diff_coeff_code={nu}",
+            f"diffusion/mom_diff_coeff_code={diff_coeff}",
+            f"diffusion/ohm_diff_coeff_code={diff_coeff}",
             f"diffusion/integrator={int_cfg}",
+            f"diffusion/rkl2_max_dt_ratio=200",
         ]
 
         return parameters
@@ -104,8 +105,8 @@ def Analyse(self, parameters):
         def get_ref(x):
             return (
                 1e-6
-                / np.sqrt(4.0 * np.pi * nu * (0.5 + tlim))
-                * np.exp(-(x**2.0) / (4.0 * nu * (0.5 + tlim)))
+                / np.sqrt(4.0 * np.pi * diff_coeff * (0.5 + tlim))
+                * np.exp(-(x**2.0) / (4.0 * diff_coeff * (0.5 + tlim)))
             )
 
         num_diff = len(diff_cfgs)

tst/regression/test_suites/diffusion_linwave3d/diffusion_linwave3d.py

@@ -0,0 +1,183 @@
+# ========================================================================================
+# AthenaPK - a performance portable block structured AMR MHD code
+# Copyright (c) 2023, Athena Parthenon Collaboration. All rights reserved.
+# Licensed under the 3-clause BSD License, see LICENSE file for details
+# ========================================================================================
+
+# Modules
+import math
+import numpy as np
+import matplotlib
+
+matplotlib.use("agg")
+import matplotlib.pylab as plt
+import sys
+import os
+import utils.test_case
+from numpy.polynomial import Polynomial
+
+""" To prevent littering up imported folders with .pyc files or __pycache_ folder"""
+sys.dont_write_bytecode = True
+
+# if this is updated make sure to update the assert statements for the number of MPI ranks, too
+lin_res = [16, 32]  # resolution for linear convergence
+
+# Upper bound on relative L1 error for each above nx1:
+error_rel_tols = [0.22, 0.05]
+# lower bound on convergence rate at final (Nx1=64) asymptotic convergence regime
+rate_tols = [2.0]  # convergence rate > 3.0 for this particular resolution, sovler
+
+method = "explicit"
+
+_nu = 0.01
+_kappa = _nu * 2.0
+_eta = _kappa
+_c_s = 0.5  # slow mode wave speed of AthenaPK linear wave configuration
+
+
+class TestCase(utils.test_case.TestCaseAbs):
+    def Prepare(self, parameters, step):
+        res = lin_res[(step - 1)]
+        # make sure we can evenly distribute the MeshBlock sizes
+        err_msg = "Num ranks must be multiples of 2 for test."
+        assert parameters.num_ranks == 1 or parameters.num_ranks % 2 == 0, err_msg
+        # ensure a minimum block size of 4
+        assert 2 * res / parameters.num_ranks >= 8, "Use <= 8 ranks for test."
+
+        mb_nx1 = (2 * res) // parameters.num_ranks
+        # ensure that nx1 is <= 128 when using scratch (V100 limit on test system)
+        while mb_nx1 > 128:
+            mb_nx1 //= 2
+
+        parameters.driver_cmd_line_args = [
+            f"parthenon/mesh/nx1={2 * res}",
+            f"parthenon/meshblock/nx1={mb_nx1}",
+            f"parthenon/mesh/nx2={res}",
+            f"parthenon/meshblock/nx2={res}",
+            f"parthenon/mesh/nx3={res}",
+            f"parthenon/meshblock/nx3={res}",
+            "parthenon/mesh/nghost=2",
+            "parthenon/time/integrator=vl2",
+            "parthenon/time/tlim=3.0",
+            "hydro/reconstruction=plm",
+            "hydro/fluid=glmmhd",
+            "hydro/riemann=hlld",
+            # enable history dump to track decay of v2 component
+            "parthenon/output2/file_type=hst",
+            "parthenon/output2/dt=0.03",
+            "problem/linear_wave/dump_max_v2=true",
+            f"parthenon/job/problem_id={res}",  # hack to rename parthenon.hst to res.hst
+            # setup linear wave (L slow mode)
+            "job/problem_id=linear_wave_mhd",
+            "problem/linear_wave/amp=1e-4",
+            "problem/linear_wave/wave_flag=2",
+            "problem/linear_wave/compute_error=false",  # done here below, not in the pgen
+            # setup diffusive processes
+            "diffusion/integrator=unsplit",
+            "diffusion/conduction=isotropic",
+            "diffusion/conduction_coeff=fixed",
+            f"diffusion/thermal_diff_coeff_code={_kappa}",
+            "diffusion/viscosity=isotropic",
+            "diffusion/viscosity_coeff=fixed",
+            f"diffusion/mom_diff_coeff_code={_nu}",
+            "diffusion/resistivity=isotropic",
+            "diffusion/resistivity_coeff=fixed",
+            f"diffusion/ohm_diff_coeff_code={_eta}",
+        ]
+
+        return parameters
+
+    def Analyse(self, parameters):
+        analyze_status = True
+
+        # Following test evaluation is adapted from the one in Athena++.
+        # This also includes the limits/tolerances set above, which are identical to Athena++.
+
+        # Lambda=1 for Athena++'s linear wave setups in 1D, 2D, and 3D:
+        L = 1.0
+        ksqr = (2.0 * np.pi / L) ** 2
+        # Equation 3.13 from Ryu, et al. (modified to add thermal conduction term)
+        # fast mode decay rate = (19\nu/4 + 3\eta + 3(\gamma-1)^2*kappa/gamma/4)*(2/15)*k^2
+        # Equation 3.14 from Ryu, et al. (modified to add thermal conduction term)
+        # slow mode decay rate = (4\nu + 3\eta/4 + 3(\gamma-1)^2*kappa/gamma)*(2/15)*k^2
+        slow_mode_rate = (
+            (4.0 * _nu + 3.0 * _eta / 4.0 + _kappa * 4.0 / 5.0) * (2.0 / 15.0) * ksqr
+        )
+
+        # Equation 3.16
+        re_num = (4.0 * np.pi**2 * _c_s) / (L * slow_mode_rate)
+        analyze_status = True
+        errors_abs = []
+
+        for nx, err_tol in zip(lin_res, error_rel_tols):
+            print(
+                "[Decaying 3D Linear Wave]: "
+                "Mesh size {} x {} x {}".format(2 * nx, nx, nx)
+            )
+            filename = os.path.join(parameters.output_path, f"{nx}.hst")
+            hst_data = np.genfromtxt(filename, names=True, skip_header=1)
+
+            tt = hst_data["1time"]
+            max_vy = hst_data["11MaxAbsV2"]
+            # estimate the decay rate from simulation, using weighted least-squares (WLS)
+            yy = np.log(np.abs(max_vy))
+            plt.plot(tt, yy)
+            plt.show()
+            p, [resid, rank, sv, rcond] = Polynomial.fit(
+                tt, yy, 1, w=np.sqrt(max_vy), full=True
+            )
+            resid_normal = np.sum((yy - p(tt)) ** 2)
+            r2 = 1 - resid_normal / (yy.size * yy.var())
+            pnormal = p.convert(domain=(-1, 1))
+            fit_rate = -pnormal.coef[-1]
+
+            error_abs = np.fabs(slow_mode_rate - fit_rate)
+            errors_abs += [error_abs]
+            error_rel = np.fabs(slow_mode_rate / fit_rate - 1.0)
+            err_rel_tol_percent = err_tol * 100.0
+
+            print(
+                "[Decaying 3D Linear Wave {}]: Reynolds number of slow mode: {}".format(
+                    method, re_num
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: R-squared of WLS regression = {}".format(
+                    method, r2
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Analytic decay rate = {}".format(
+                    method, slow_mode_rate
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Measured decay rate = {}".format(
+                    method, fit_rate
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Decay rate absolute error = {}".format(
+                    method, error_abs
+                )
+            )
+            print(
+                "[Decaying 3D Linear Wave {}]: Decay rate relative error = {}".format(
+                    method, error_rel
+                )
+            )
+
+            if error_rel > err_tol:
+                print(
+                    "WARN [Decaying 3D Linear Wave {}]: decay rate disagrees"
+                    " with prediction by >{}%".format(method, err_rel_tol_percent)
+                )
+                analyze_status = False
+            else:
+                print(
+                    "[Decaying 3D Linear Wave {}]: decay rate is within "
+                    "{}% of analytic value".format(method, err_rel_tol_percent)
+                )
+                print("")
+
+        return analyze_status