From 5fea7126e5e4ac0babcafe86212d5f6043ca6972 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 28 Sep 2023 16:48:29 +0200 Subject: [PATCH] Add viscosity convergence test --- .github/workflows/ci.yml | 1 + inputs/viscosity.in | 67 +++++++ src/hydro/diffusion/viscosity.cpp | 4 +- src/pgen/visc.cpp | 5 +- tst/regression/CMakeLists.txt | 3 + .../test_suites/viscosity/__init__.py | 0 .../test_suites/viscosity/viscosity.py | 166 ++++++++++++++++++ 7 files changed, 241 insertions(+), 5 deletions(-) create mode 100644 inputs/viscosity.in create mode 100644 tst/regression/test_suites/viscosity/__init__.py create mode 100644 tst/regression/test_suites/viscosity/viscosity.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 05049ba7..f48c0ce8 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -54,6 +54,7 @@ jobs: build/tst/regression/outputs/cluster_tabular_cooling/convergence.png build/tst/regression/outputs/aniso_therm_cond_ring_conv/ring_convergence.png build/tst/regression/outputs/aniso_therm_cond_gauss_conv/cond.png + build/tst/regression/outputs/viscosity/visc.png build/tst/regression/outputs/field_loop/field_loop.png build/tst/regression/outputs/riemann_hydro/shock_tube.png build/tst/regression/outputs/turbulence/parthenon.hst diff --git a/inputs/viscosity.in b/inputs/viscosity.in new file mode 100644 index 00000000..6c652937 --- /dev/null +++ b/inputs/viscosity.in @@ -0,0 +1,67 @@ +# AthenaPK - a performance portable block structured AMR MHD code +# Copyright (c) 2021, Athena Parthenon Collaboration. All rights reserved. +# Licensed under the BSD 3-Clause License (the "LICENSE"); + + +problem = Thermal diffusion setup + + +problem_id = visc + + +t0 = 0.5 # Temporal offset for initial Gaussian profile +amp = 1e-6 # Amplitude of Gaussian profile + + +refinement = none +nghost = 2 + +nx1 = 256 +x1min = -6.0 +x1max = 6.0 +ix1_bc = outflow +ox1_bc = outflow + +nx2 = 1 +x2min = -1.0 +x2max = 1.0 +ix2_bc = periodic +ox2_bc = periodic + +nx3 = 1 +x3min = -1.0 +x3max = 1.0 +ix3_bc = periodic +ox3_bc = periodic + + +nx1=32 +nx2=1 +nx3=1 + + +integrator = rk2 +cfl = 0.8 +tlim = 2.0 +nlim = -1 + + +fluid = euler +riemann = none +eos = adiabatic +reconstruction = dc +gamma = 1.4 +calc_dt_hytp = true # need in combination with rkl2_max_dt_ratio to limit number of stages in RKL2 + + +integrator = unsplit +viscosity = isotropic +viscosity_coeff = fixed +mom_diff_coeff_code = 0.25 +rkl2_max_dt_ratio = 200.0 + + +file_type = hdf5 +dt = 2.0 +variables = prim +id = prim diff --git a/src/hydro/diffusion/viscosity.cpp b/src/hydro/diffusion/viscosity.cpp index 8ef3cbec..cb23d85a 100644 --- a/src/hydro/diffusion/viscosity.cpp +++ b/src/hydro/diffusion/viscosity.cpp @@ -295,6 +295,4 @@ void MomentumDiffFluxIsoFixed(MeshData *md) { //! TODO(pgrete) Calculate thermal conduction, general case, i.e., anisotropic and/or with //! varying (incl. saturated) coefficient -void MomentumDiffFluxGeneral(MeshData *md) { - PARTHENON_THROW("Needs impl."); -} +void MomentumDiffFluxGeneral(MeshData *md) { PARTHENON_THROW("Needs impl."); } diff --git a/src/pgen/visc.cpp b/src/pgen/visc.cpp index 7144c341..dd7e7ea5 100644 --- a/src/pgen/visc.cpp +++ b/src/pgen/visc.cpp @@ -27,10 +27,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto &mbd = pmb->meshblock_data.Get(); auto &u = mbd->Get("cons").data; - Real gm1 = pin->GetReal("hydro", "gamma") - 1.0; + auto gamma = pin->GetReal("hydro", "gamma"); + auto gm1 = gamma - 1.0; Real v1 = 0., v2 = 0., v3 = 0.; Real d0 = 1.; - Real p0 = 1.; + Real p0 = d0 / gamma; // i.e., c_s = 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"); diff --git a/tst/regression/CMakeLists.txt b/tst/regression/CMakeLists.txt index c43b4de5..422251de 100644 --- a/tst/regression/CMakeLists.txt +++ b/tst/regression/CMakeLists.txt @@ -49,6 +49,9 @@ setup_test_both("aniso_therm_cond_ring_multid" "--driver ${PROJECT_BINARY_DIR}/b setup_test_both("aniso_therm_cond_gauss_conv" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \ --driver_input ${PROJECT_SOURCE_DIR}/inputs/diffusion.in --num_steps 24" "convergence") +setup_test_both("viscosity" "--driver ${PROJECT_BINARY_DIR}/bin/athenaPK \ + --driver_input ${PROJECT_SOURCE_DIR}/inputs/viscosity.in --num_steps 6" "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/viscosity/__init__.py b/tst/regression/test_suites/viscosity/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tst/regression/test_suites/viscosity/viscosity.py b/tst/regression/test_suites/viscosity/viscosity.py new file mode 100644 index 00000000..b6362793 --- /dev/null +++ b/tst/regression/test_suites/viscosity/viscosity.py @@ -0,0 +1,166 @@ +# ======================================================================================== +# 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 itertools +import utils.test_case +from scipy.optimize import curve_fit + +# To prevent littering up imported folders with .pyc files or __pycache_ folder +sys.dont_write_bytecode = True + +int_cfgs = ["unsplit", "rkl2"] +res_cfgs = [256, 512, 1024] +tlim = 2.0 +nu = 0.25 + +all_cfgs = list(itertools.product(res_cfgs, int_cfgs)) + + +def get_outname(all_cfg): + res, int_cfg = all_cfg + return f"{res}_{int_cfg}" + + +class TestCase(utils.test_case.TestCaseAbs): + def Prepare(self, parameters, step): + assert parameters.num_ranks <= 4, "Use <= 4 ranks for diffusion test." + + res, int_cfg = all_cfgs[step - 1] + + outname = get_outname(all_cfgs[step - 1]) + + parameters.driver_cmd_line_args = [ + "parthenon/mesh/nx1=%d" % res, + "parthenon/meshblock/nx1=64", + "parthenon/mesh/x1min=-6.0", + "parthenon/mesh/x1max=6.0", + "parthenon/mesh/nx2=1", + "parthenon/meshblock/nx2=1", + "parthenon/mesh/x2min=-1.0", + "parthenon/mesh/x2max=1.0", + "parthenon/mesh/nx3=1", + "parthenon/meshblock/nx3=1", + "parthenon/output0/id=%s" % outname, + "parthenon/time/tlim=%f" % tlim, + "diffusion/mom_diff_coeff_code=%f" % nu, + "diffusion/integrator=%s" % int_cfg, + ] + + return parameters + + def Analyse(self, parameters): + sys.path.insert( + 1, + parameters.parthenon_path + + "/scripts/python/packages/parthenon_tools/parthenon_tools", + ) + + try: + import phdf + except ModuleNotFoundError: + print("Couldn't find module to read Parthenon hdf5 files.") + return False + + tests_passed = True + + 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))) + ) + + num_rows = len(res_cfgs) + num_cols = len(int_cfgs) + fig, p = plt.subplots(num_rows + 1, 2, sharey="row", sharex="row") + + l1_err = np.zeros((len(int_cfgs), len(res_cfgs))) + for step in range(len(all_cfgs)): + outname = get_outname(all_cfgs[step]) + data_filename = f"{parameters.output_path}/parthenon.{outname}.final.phdf" + data_file = phdf.phdf(data_filename) + prim = data_file.Get("prim") + zz, yy, xx = data_file.GetVolumeLocations() + mask = yy == yy[0] + v2 = prim[2][mask] + x = xx[mask] + res, int_cfg = all_cfgs[step] + row = res_cfgs.index(res) + col = int_cfgs.index(int_cfg) + + v2_ref = get_ref(x) + l1 = np.average(np.abs(v2 - v2_ref)) + l1_err[ + int_cfgs.index(int_cfg), + res_cfgs.index(res), + ] = l1 + p[row, col].plot(x, v2, label=f"N={res} L$_1$={l1:.2g}") + + # Plot convergence + for j, int_cfg in enumerate(int_cfgs): + p[0, j].set_title(f"Integrator: {int_cfg}") + + p[-1, j].plot( + res_cfgs, + l1_err[j, :], + label=f"data", + ) + + # Simple convergence estimator + conv_model = lambda log_n, log_a, conv_rate: conv_rate * log_n + log_a + popt, pconv = curve_fit(conv_model, np.log(res_cfgs), np.log(l1_err[j, :])) + conv_a, conv_measured = popt + # Note that the RKL2 convergence on the plots is currently significantly better + # than expected (<-3) though the L1 errors themself are larger than the unsplit + # integrator (as expected). + # For a more reasonable test (which would take longer), reduce the RKL2 ratio to, + # say, 200 and extend the resolution grid to 1024 (as the first data point at N=128 + # is comparatively worse than at N>128). + if conv_measured > -1.98: + print( + f"!!!\nConvergence for test with {int_cfg} integrator " + f"is worse ({conv_measured}) than expected (-1.98).\n!!!" + ) + tests_passed = False + p[-1, j].plot( + res_cfgs, + np.exp(conv_a) * res_cfgs**conv_measured, + ":", + lw=0.75, + label=f"Measured conv: {conv_measured:.2f}", + ) + + p[-1, 0].set_xscale("log") + p[-1, 0].set_yscale("log") + p[-1, 0].legend(fontsize=6) + p[-1, 1].legend(fontsize=6) + + # Plot reference lines + x = np.linspace(-6, 6, 400) + for i in range(num_rows): + for j in range(num_cols): + y = get_ref(x) + p[i, j].plot(x, y, "-", lw=0.5, color="black", alpha=0.8) + p[i, j].grid() + p[i, j].legend(fontsize=6) + + fig.tight_layout() + fig.savefig( + os.path.join(parameters.output_path, "visc.png"), + bbox_inches="tight", + dpi=300, + ) + + return tests_passed