From 969f3480b11a9675accea86295f341f469a47386 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 9 Nov 2022 10:22:15 +0100 Subject: [PATCH] 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_