From ce255157d3935f80813ecf695cbd92b21f190454 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 8 Nov 2022 14:24:39 +0100 Subject: [PATCH 1/5] Fix static MHD convergence test --- .../mhd_convergence/mhd_convergence.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/tst/regression/test_suites/mhd_convergence/mhd_convergence.py b/tst/regression/test_suites/mhd_convergence/mhd_convergence.py index 18fea9fc..9d3de08e 100644 --- a/tst/regression/test_suites/mhd_convergence/mhd_convergence.py +++ b/tst/regression/test_suites/mhd_convergence/mhd_convergence.py @@ -98,6 +98,7 @@ def Prepare(self, parameters, step): mb_nx1 //= 2 parameters.driver_cmd_line_args = [ + "job/problem_id=linear_wave_mhd", "parthenon/mesh/nx1=%d" % (2 * res), "parthenon/meshblock/nx1=%d" % mb_nx1, "parthenon/mesh/nx2=%d" % res, @@ -164,7 +165,7 @@ def Analyse(self, parameters): ) # quick and dirty test - if data[47, 4] > 6.14e-12: + if data[47, 4] > 4.05e-10: print("QUICK AND DIRTY TEST FAILED") analyze_status = False @@ -182,21 +183,18 @@ def Analyse(self, parameters): ), ) - plt.plot([32, 512], [7e-7, 7e-7 / (512 / 32)], "--", label="first order") + plt.plot([32, 256], [1.1e-6, 1.1e-6 / (256 / 32)], "--", label="first order") plt.plot( - [32, 512], [1.7e-7, 1.7e-7 / (512 / 32) ** 2], "--", label="second order" + [32, 256], [1.7e-7, 1.7e-7 / (256 / 32) ** 2], "--", label="second order" ) plt.plot( - [32, 512], [3.7e-8, 3.7e-8 / (512 / 32) ** 2], "--", label="second order" + [32, 256], [3.1e-8, 3.1e-8 / (256 / 32) ** 2], "--", label="second order" ) plt.plot( - [32, 512], [5.6e-8, 5.6e-8 / (512 / 32) ** 3], "--", label="third order" - ) - plt.plot( - [32, 512], [3.6e-9, 3.6e-9 / (512 / 32) ** 3], "--", label="third order" + [32, 256], [8.3e-8, 8.3e-8 / (256 / 32) ** 3], "--", label="third order" ) - plt.ylim(1e-12, 5e-6) + plt.ylim(3e-10, 2e-6) plt.legend(bbox_to_anchor=(1, 1), loc="upper left") plt.xscale("log") From fdbc4d0b2c9b9bd03f5c5506e3ab196dd9c524ee Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 8 Nov 2022 15:38:46 +0100 Subject: [PATCH 2/5] Fix lin wave 2d setup --- src/pgen/linear_wave_mhd.cpp | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/src/pgen/linear_wave_mhd.cpp b/src/pgen/linear_wave_mhd.cpp index 8a7a6a80..96708bb6 100644 --- a/src/pgen/linear_wave_mhd.cpp +++ b/src/pgen/linear_wave_mhd.cpp @@ -157,6 +157,8 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) { Eigensystem(d0, u0, v0, w0, h0, bx0, by0, bz0, xfact, yfact, ev, rem, lem); + std::cout << "rem=" << rem[0][wave_flag] << std::endl; + // TODO(pgrete) see how to get access to the SimTime object outside the driver // if (pin->GetOrAddBoolean("problem/linear_wave", "test", false) && ncycle == 0) { if (pin->GetOrAddBoolean("problem/linear_wave", "test", false)) { @@ -366,9 +368,6 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm) //======================================================================================== 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); // Initialize the magnetic fields. Note wavevector, eigenvectors, and other variables // are set in InitUserMeshData @@ -392,9 +391,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto &coords = pmb->coords; // Initialize components of the vector potential - for (int k = kb.s - 1; k <= kb.e + 1; k++) { - for (int j = jb.s - 1; j <= jb.e + 1; j++) { - for (int i = ib.s - 1; i <= ib.e + 1; i++) { + IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire); + IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire); + IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire); + for (int k = kb.s; k <= kb.e; k++) { + for (int j = jb.s; j <= jb.e; j++) { + for (int i = ib.s; i <= ib.e; i++) { a1(k, j, i) = A1(coords.x1v(i), coords.x2v(j), coords.x3v(k)); a2(k, j, i) = A2(coords.x1v(i), coords.x2v(j), coords.x3v(k)); a3(k, j, i) = A3(coords.x1v(i), coords.x2v(j), coords.x3v(k)); @@ -407,6 +409,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { auto &u_dev = rc->Get("cons").data; // initializing on host auto u = u_dev.GetHostMirrorAndCopy(); + ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior); + jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior); + kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior); + + const bool three_d = pmb->cellbounds.ncellsk(IndexDomain::entire) > 1; for (int k = kb.s; k <= kb.e; k++) { for (int j = jb.s; j <= jb.e; j++) { for (int i = ib.s; i <= ib.e; i++) { @@ -422,10 +429,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) { u(IM2, k, j, i) = mx * cos_a2 * sin_a3 + my * cos_a3 - mz * sin_a2 * sin_a3; u(IM3, k, j, i) = mx * sin_a2 + mz * cos_a2; - u(IB1, k, j, i) = (a3(k, j + 1, i) - a3(k, j - 1, i)) / coords.dx2v(j) / 2.0 - - (a2(k + 1, j, i) - a2(k - 1, j, i)) / coords.dx3v(k) / 2.0; - u(IB2, k, j, i) = (a1(k + 1, j, i) - a1(k - 1, j, i)) / coords.dx3v(k) / 2.0 - - (a3(k, j, i + 1) - a3(k, j, i - 1)) / coords.dx1v(i) / 2.0; + u(IB1, k, j, i) = (a3(k, j + 1, i) - a3(k, j - 1, i)) / coords.dx2v(j) / 2.0; + u(IB1, k, j, i) -= + three_d ? (a2(k + 1, j, i) - a2(k - 1, j, i)) / coords.dx3v(k) / 2.0 : 0.0; + u(IB2, k, j, i) = -(a3(k, j, i + 1) - a3(k, j, i - 1)) / coords.dx1v(i) / 2.0; + u(IB2, k, j, i) += + three_d ? (a1(k + 1, j, i) - a1(k - 1, j, i)) / coords.dx3v(k) / 2.0 : 0.0; u(IB3, k, j, i) = (a2(k, j, i + 1) - a2(k, j, i - 1)) / coords.dx1v(i) / 2.0 - (a1(k, j + 1, i) - a1(k, j - 1, i)) / coords.dx2v(j) / 2.0; From e79f7e782d916596928edcfef7126e91fd472c4d Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 9 Nov 2022 09:34:41 +0100 Subject: [PATCH 3/5] Bump parth and update interfaces --- external/parthenon | 2 +- src/hydro/hydro_driver.cpp | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/external/parthenon b/external/parthenon index 22f277ac..3930f9f8 160000 --- a/external/parthenon +++ b/external/parthenon @@ -1 +1 @@ -Subproject commit 22f277ac80a155ca81c740a14d44bc027264aba3 +Subproject commit 3930f9f897257802164378c01be885ca61093631 diff --git a/src/hydro/hydro_driver.cpp b/src/hydro/hydro_driver.cpp index 4b67a2bf..5d86dfd4 100644 --- a/src/hydro/hydro_driver.cpp +++ b/src/hydro/hydro_driver.cpp @@ -12,7 +12,7 @@ // Parthenon headers #include "bvals/cc/bvals_cc_in_one.hpp" -#include "mesh/refinement_cc_in_one.hpp" +#include "mesh/refinement_in_one.hpp" #include "refinement/refinement.hpp" #include // AthenaPK headers @@ -275,8 +275,7 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) { tl.AddTask(none, parthenon::cell_centered_bvars::ReceiveBoundBufs, mu0); auto set = tl.AddTask(recv, parthenon::cell_centered_bvars::SetBounds, mu0); if (pmesh->multilevel) { - tl.AddTask(set, parthenon::cell_centered_refinement::RestrictPhysicalBounds, - mu0.get()); + tl.AddTask(set, parthenon::cell_centered_bvars::RestrictGhostHalos, mu0, false); } } From 969f3480b11a9675accea86295f341f469a47386 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 9 Nov 2022 10:22:15 +0100 Subject: [PATCH 4/5] Add skel for custom refine op --- src/CMakeLists.txt | 1 + src/hydro/hydro.cpp | 3 + src/hydro/prolongation/custom_ops.hpp | 108 ++++++++++++++++++++++++++ 3 files changed, 112 insertions(+) create mode 100644 src/hydro/prolongation/custom_ops.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2066812d..76d5c69e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,6 +12,7 @@ add_executable( hydro/hydro_driver.cpp hydro/hydro.cpp hydro/glmmhd/dedner_source.cpp + hydro/prolongation/custom_ops.hpp hydro/srcterms/gravitational_field.hpp hydro/srcterms/tabular_cooling.hpp hydro/srcterms/tabular_cooling.cpp diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index f5ba44c6..905bfb24 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -30,6 +30,7 @@ #include "glmmhd/glmmhd.hpp" #include "hydro.hpp" #include "outputs/outputs.hpp" +#include "prolongation/custom_ops.hpp" #include "rsolvers/rsolvers.hpp" #include "srcterms/tabular_cooling.hpp" #include "utils/error_checking.hpp" @@ -541,6 +542,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { Metadata m( {Metadata::Cell, Metadata::Independent, Metadata::FillGhost, Metadata::WithFluxes}, std::vector({nhydro + nscalars}), cons_labels); + m.RegisterRefinementOps(); pkg->AddField("cons", m); m = Metadata({Metadata::Cell, Metadata::Derived}, std::vector({nhydro + nscalars}), diff --git a/src/hydro/prolongation/custom_ops.hpp b/src/hydro/prolongation/custom_ops.hpp new file mode 100644 index 00000000..b8a53ee6 --- /dev/null +++ b/src/hydro/prolongation/custom_ops.hpp @@ -0,0 +1,108 @@ +//======================================================================================== +// AthenaPK - a performance portable block structured AMR astrophysical MHD code. +// Copyright (c) 2022, Athena-Parthenon Collaboration. All rights reserved. +// Licensed under the BSD 3-Clause License (the "LICENSE"). +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2020-2022 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 2021. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do +// so. +//======================================================================================== + +#ifndef HYDRO_PROLONGATION_CUSTOM_OPS_HPP_ +#define HYDRO_PROLONGATION_CUSTOM_OPS_HPP_ + +#include +#include + +#include "coordinates/coordinates.hpp" // for coordinates +#include "kokkos_abstraction.hpp" // ParArray +#include "mesh/domain.hpp" // for IndesShape +#include "mesh/mesh_refinement_ops.hpp" + +namespace Hydro { +namespace refinement_ops { + +using parthenon::Coordinates_t; +using parthenon::ParArray6D; + +template +struct ProlongateCellMinModMultiD { + KOKKOS_FORCEINLINE_FUNCTION static void + Do(const int l, const int m, const int n, const int k, const int j, const int i, + const IndexRange &ckb, const IndexRange &cjb, const IndexRange &cib, + const IndexRange &kb, const IndexRange &jb, const IndexRange &ib, + const Coordinates_t &coords, const Coordinates_t &coarse_coords, + const ParArray6D *pcoarse, const ParArray6D *pfine) { + using namespace parthenon::refinement_ops::util; + auto &coarse = *pcoarse; + auto &fine = *pfine; + + const Real fc = coarse(l, m, n, k, j, i); + + int fi; + Real dx1fm, dx1fp, dx1m, dx1p; + GetGridSpacings<1>(coords, coarse_coords, cib, ib, i, &fi, &dx1m, &dx1p, &dx1fm, + &dx1fp); + const Real gx1c = GradMinMod(fc, coarse(l, m, n, k, j, i - 1), + coarse(l, m, n, k, j, i + 1), dx1m, dx1p); + + int fj = jb.s; // overwritten as needed + Real dx2fm = 0; + Real dx2fp = 0; + Real gx2c = 0; + if constexpr (DIM > 1) { + Real dx2m, dx2p; + GetGridSpacings<2>(coords, coarse_coords, cjb, jb, j, &fj, &dx2m, &dx2p, &dx2fm, + &dx2fp); + gx2c = GradMinMod(fc, coarse(l, m, n, k, j - 1, i), coarse(l, m, n, k, j + 1, i), + dx2m, dx2p); + } + int fk = kb.s; + Real dx3fm = 0; + Real dx3fp = 0; + Real gx3c = 0; + if constexpr (DIM > 2) { + Real dx3m, dx3p; + GetGridSpacings<3>(coords, coarse_coords, ckb, kb, k, &fk, &dx3m, &dx3p, &dx3fm, + &dx3fp); + gx3c = GradMinMod(fc, coarse(l, m, n, k - 1, j, i), coarse(l, m, n, k + 1, j, i), + dx3m, dx3p); + } + + // KGF: add the off-centered quantities first to preserve FP symmetry + // JMM: Extraneous quantities are zero + fine(l, m, n, fk, fj, fi) = fc - (gx1c * dx1fm + gx2c * dx2fm + gx3c * dx3fm); + fine(l, m, n, fk, fj, fi + 1) = fc + (gx1c * dx1fp - gx2c * dx2fm - gx3c * dx3fm); + if constexpr (DIM > 1) { + fine(l, m, n, fk, fj + 1, fi) = fc - (gx1c * dx1fm - gx2c * dx2fp + gx3c * dx3fm); + fine(l, m, n, fk, fj + 1, fi + 1) = + fc + (gx1c * dx1fp + gx2c * dx2fp - gx3c * dx3fm); + } + if constexpr (DIM > 2) { + fine(l, m, n, fk + 1, fj, fi) = fc - (gx1c * dx1fm + gx2c * dx2fm - gx3c * dx3fp); + fine(l, m, n, fk + 1, fj, fi + 1) = + fc + (gx1c * dx1fp - gx2c * dx2fm + gx3c * dx3fp); + fine(l, m, n, fk + 1, fj + 1, fi) = + fc - (gx1c * dx1fm - gx2c * dx2fp - gx3c * dx3fp); + fine(l, m, n, fk + 1, fj + 1, fi + 1) = + fc + (gx1c * dx1fp + gx2c * dx2fp + gx3c * dx3fp); + } + } +}; +} // namespace refinement_ops +} // namespace Hydro + +#endif // HYDRO_PROLONGATION_CUSTOM_OPS_HPP_ From 3be50508da6b2a542148eabb5e8c9550295fc676 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 9 Nov 2022 15:39:04 +0100 Subject: [PATCH 5/5] Limit multi-d prolong slopes --- src/hydro/prolongation/custom_ops.hpp | 42 +++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/src/hydro/prolongation/custom_ops.hpp b/src/hydro/prolongation/custom_ops.hpp index b8a53ee6..d03ca9d1 100644 --- a/src/hydro/prolongation/custom_ops.hpp +++ b/src/hydro/prolongation/custom_ops.hpp @@ -56,8 +56,8 @@ struct ProlongateCellMinModMultiD { Real dx1fm, dx1fp, dx1m, dx1p; GetGridSpacings<1>(coords, coarse_coords, cib, ib, i, &fi, &dx1m, &dx1p, &dx1fm, &dx1fp); - const Real gx1c = GradMinMod(fc, coarse(l, m, n, k, j, i - 1), - coarse(l, m, n, k, j, i + 1), dx1m, dx1p); + Real gx1c = GradMinMod(fc, coarse(l, m, n, k, j, i - 1), coarse(l, m, n, k, j, i + 1), + dx1m, dx1p); int fj = jb.s; // overwritten as needed Real dx2fm = 0; @@ -82,6 +82,44 @@ struct ProlongateCellMinModMultiD { dx3m, dx3p); } + // Max. expected total difference. (dx#fm/p are positive by construction) + Real dqmax = std::abs(gx1c) * std::max(dx1fm, dx1p); + int jlim = 0; + int klim = 0; + if constexpr (DIM > 1) { + dqmax += std::abs(gx2c) * std::max(dx2fm, dx2fp); + jlim = 1; + } + if constexpr (DIM > 2) { + dqmax += std::abs(gx3c) * std::max(dx3fm, dx3fp); + klim = 1; + } + // Min/max values of all coarse cells used here + Real qmin = fc; + Real qmax = fc; + for (int koff = -klim; koff <= klim; koff++) { + for (int joff = -jlim; joff <= jlim; joff++) { + for (int ioff = -1; ioff <= 1; ioff++) { + qmin = std::min(qmin, coarse(l, m, n, k + koff, j + joff, i + ioff)); + qmax = std::max(qmax, coarse(l, m, n, k + koff, j + joff, i + ioff)); + } + } + } + + // Scaling factor to limit all slopes simultaneously + Real alpha = 1.0; + if (dqmax * alpha > (qmax - fc)) { + alpha = (qmax - fc) / dqmax; + } + if (dqmax * alpha > (fc - qmin)) { + alpha = (fc - qmin) / dqmax; + } + + // Ensure no new extrema are introduced in multi-D + gx1c *= alpha; + gx2c *= alpha; + gx3c *= alpha; + // KGF: add the off-centered quantities first to preserve FP symmetry // JMM: Extraneous quantities are zero fine(l, m, n, fk, fj, fi) = fc - (gx1c * dx1fm + gx2c * dx2fm + gx3c * dx3fm);