From 6e0151e6aa57f054be9c0383605b9a8b2204850a Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 5 Oct 2023 12:00:31 +0200 Subject: [PATCH] Add linwave3d decay diffusion test and fix parabolic dt --- inputs/diffusion.in | 2 +- src/hydro/diffusion/conduction.cpp | 4 +- src/hydro/diffusion/resistivity.cpp | 3 +- src/hydro/diffusion/viscosity.cpp | 3 +- src/hydro/hydro.cpp | 4 + src/main.cpp | 1 + src/pgen/linear_wave_mhd.cpp | 34 ++++ src/pgen/pgen.hpp | 1 + tst/regression/CMakeLists.txt | 3 + .../test_suites/diffusion/diffusion.py | 11 +- .../diffusion_linwave3d/__init__.py | 0 .../diffusion_linwave3d.py | 183 ++++++++++++++++++ 12 files changed, 239 insertions(+), 10 deletions(-) create mode 100644 tst/regression/test_suites/diffusion_linwave3d/__init__.py create mode 100644 tst/regression/test_suites/diffusion_linwave3d/diffusion_linwave3d.py diff --git a/inputs/diffusion.in b/inputs/diffusion.in index d7ce3da8..71eeae67 100644 --- a/inputs/diffusion.in +++ b/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 file_type = hdf5 diff --git a/src/hydro/diffusion/conduction.cpp b/src/hydro/diffusion/conduction.cpp index d6e57d30..40256024 100644 --- a/src/hydro/diffusion/conduction.cpp +++ b/src/hydro/diffusion/conduction.cpp @@ -179,8 +179,8 @@ Real EstimateConductionTimestep(MeshData *md) { }, Kokkos::Min(min_dt_cond)); } - - return fac * min_dt_cond; + const auto &cfl_diff = hydro_pkg->Param("cfl_diff"); + return cfl_diff * fac * min_dt_cond; } //--------------------------------------------------------------------------------------- diff --git a/src/hydro/diffusion/resistivity.cpp b/src/hydro/diffusion/resistivity.cpp index d351d210..ab10bc61 100644 --- a/src/hydro/diffusion/resistivity.cpp +++ b/src/hydro/diffusion/resistivity.cpp @@ -84,7 +84,8 @@ Real EstimateResistivityTimestep(MeshData *md) { PARTHENON_THROW("Needs impl."); } - return fac * min_dt_resist; + const auto &cfl_diff = hydro_pkg->Param("cfl_diff"); + return cfl_diff * fac * min_dt_resist; } //--------------------------------------------------------------------------------------- diff --git a/src/hydro/diffusion/viscosity.cpp b/src/hydro/diffusion/viscosity.cpp index cb23d85a..aa552f14 100644 --- a/src/hydro/diffusion/viscosity.cpp +++ b/src/hydro/diffusion/viscosity.cpp @@ -88,7 +88,8 @@ Real EstimateViscosityTimestep(MeshData *md) { PARTHENON_THROW("Needs impl."); } - return fac * min_dt_visc; + const auto &cfl_diff = hydro_pkg->Param("cfl_diff"); + return cfl_diff * fac * min_dt_visc; } //--------------------------------------------------------------------------------------- diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 06a9e39f..e2062d6c 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -665,6 +665,10 @@ std::shared_ptr Initialize(ParameterInput *pin) { } if (diffint != DiffInt::none) { pkg->AddParam("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("cfl")); + pkg->AddParam<>("cfl_diff", cfl_diff); } pkg->AddParam<>("diffint", diffint); diff --git a/src/main.cpp b/src/main.cpp index 93f984bb..fa0c4dec 100644 --- a/src/main.cpp +++ b/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; diff --git a/src/pgen/linear_wave_mhd.cpp b/src/pgen/linear_wave_mhd.cpp index 7f1b549b..75f69663 100644 --- a/src/pgen/linear_wave_mhd.cpp +++ b/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 *md) { + auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro"); + + const auto &prim_pack = md->PackVariables(std::vector{"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>( + 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(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::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 diff --git a/src/pgen/pgen.hpp b/src/pgen/pgen.hpp index 182d8497..2aeee63c 100644 --- a/src/pgen/pgen.hpp +++ b/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 { diff --git a/tst/regression/CMakeLists.txt b/tst/regression/CMakeLists.txt index 5404dec4..f4fa246e 100644 --- a/tst/regression/CMakeLists.txt +++ b/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") diff --git a/tst/regression/test_suites/diffusion/diffusion.py b/tst/regression/test_suites/diffusion/diffusion.py index d0e7dfc8..4927ebe1 100644 --- a/tst/regression/test_suites/diffusion/diffusion.py +++ b/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) diff --git a/tst/regression/test_suites/diffusion_linwave3d/__init__.py b/tst/regression/test_suites/diffusion_linwave3d/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tst/regression/test_suites/diffusion_linwave3d/diffusion_linwave3d.py b/tst/regression/test_suites/diffusion_linwave3d/diffusion_linwave3d.py new file mode 100644 index 00000000..07f2484a --- /dev/null +++ b/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