From 318a254e9162677dd0e634556abed0ec39508ad9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 13:49:43 -0400 Subject: [PATCH 1/7] update to 24.11 --- external/Microphysics | 2 +- external/amrex | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/external/Microphysics b/external/Microphysics index 533f712d2d..1753bb908e 160000 --- a/external/Microphysics +++ b/external/Microphysics @@ -1 +1 @@ -Subproject commit 533f712d2d217ea4f4063f3a809eb16a839c9ef6 +Subproject commit 1753bb908e2a966e66a3a79872f4c496c1df9fd3 diff --git a/external/amrex b/external/amrex index aba5f01be1..4a6e7c87e1 160000 --- a/external/amrex +++ b/external/amrex @@ -1 +1 @@ -Subproject commit aba5f01be1f0dbd93600671f6d26df228b7cca12 +Subproject commit 4a6e7c87e1494d10893b36f0d7ce6bf9df19fe0d From 32e14bff48d1bc6a11f52b422bcffaa47e8eb5b6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 13:50:26 -0400 Subject: [PATCH 2/7] update changes for 24.11 (#2985) --- CHANGES.md | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index dd315bc074..ed65cbc0b4 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,24 @@ +# 24.11 + + * a new well-balanced method was added to the CTU PPM solver. This + does the characteristic projection only on the perturbed pressure + and then adds back in the hydrostatic pressure. It can be enabled + via `castro.ppm_well_balanced` (#2945) + + * fixed a bug in the div{U} calculation for artificial viscosity on + symmetry boundaries (#2983) + + * more development on 2D spherical geometry (#2973, #2975, #2981) + + * updates to the massive star plotting scripts (#2979) + + * add some new checks to prevent running unsupported combinations + of solvers (#2978) + + * documentation updates (#2977) + + * `flame_wave` can now be run in 1D (#2976) + # 24.10 * update initial model for `subchandra` when doing ASE NSE (#2970) From b38a018a1ee36964348a686a5dd8261a3056d050 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Fri, 8 Nov 2024 09:53:38 -0500 Subject: [PATCH 3/7] name variable to geom when its amrex::Geometry to avoid confusion (#2989) --- Source/driver/Derive.H | 70 ++++++++++++------------ Source/driver/Derive.cpp | 114 +++++++++++++++++++-------------------- 2 files changed, 92 insertions(+), 92 deletions(-) diff --git a/Source/driver/Derive.H b/Source/driver/Derive.H index 7dd302da16..0533c2d824 100644 --- a/Source/driver/Derive.H +++ b/Source/driver/Derive.H @@ -12,125 +12,125 @@ extern "C" void ca_derpres (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dereint1 (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dereint2 (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derlogden (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_deruplusc (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_deruminusc (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dersoundspeed (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dergamma1 (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermachnumber (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derentropy (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #ifdef DIFFUSION void ca_dercond (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdiffcoeff (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdiffterm (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #endif void ca_derenuc (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derenuctimescale (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dervel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagvel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermaggrav (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derradialvel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dercircvel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagmom (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derangmomx (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_derangmomy (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_derangmomz (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_derkineng (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_dernull @@ -144,53 +144,53 @@ extern "C" void ca_derspec (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derabar (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derye (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagvort (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdivu (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derstate (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #ifdef MHD void ca_dermagcenx (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagceny (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagcenz (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdivb (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #endif diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index 3064fc3dd2..63d31449d1 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -21,7 +21,7 @@ extern "C" // need to explicitly synchronize after GPU kernels. void ca_derpres(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -54,7 +54,7 @@ extern "C" } void ca_dereint1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -76,7 +76,7 @@ extern "C" } void ca_dereint2(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -92,7 +92,7 @@ extern "C" } void ca_derlogden(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -107,7 +107,7 @@ extern "C" } void ca_deruplusc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -142,7 +142,7 @@ extern "C" } void ca_deruminusc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -177,7 +177,7 @@ extern "C" } void ca_dersoundspeed(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -213,7 +213,7 @@ extern "C" void ca_dergamma1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -248,7 +248,7 @@ extern "C" } void ca_dermachnumber(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -286,7 +286,7 @@ extern "C" } void ca_derentropy(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -321,7 +321,7 @@ extern "C" #ifdef DIFFUSION void ca_dercond(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -333,7 +333,7 @@ extern "C" } void ca_derdiffcoeff(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -345,7 +345,7 @@ extern "C" } void ca_derdiffterm(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -360,9 +360,9 @@ extern "C" fill_temp_cond(obx, dat, coeff_arr); - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); - const int coord_type = geomdata.Coord(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); + const int coord_type = geom.Coord(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -432,14 +432,14 @@ extern "C" #ifdef REACTIONS void ca_derenuctimescale(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); Real dd = 0.0_rt; #if AMREX_SPACEDIM == 1 @@ -494,7 +494,7 @@ extern "C" } void ca_derenuc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); @@ -512,7 +512,7 @@ extern "C" #endif void ca_dervel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -528,7 +528,7 @@ extern "C" void ca_dermagvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -549,7 +549,7 @@ extern "C" void ca_dermaggrav(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -568,7 +568,7 @@ extern "C" } void ca_derradialvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -577,8 +577,8 @@ extern "C" auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -623,7 +623,7 @@ extern "C" void ca_dercircvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -632,8 +632,8 @@ extern "C" auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -688,7 +688,7 @@ extern "C" void ca_dermagmom(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -707,13 +707,13 @@ extern "C" } void ca_derangmomx (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { int idir = 0; - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); auto const dat = datfab.array(); auto const L = derfab.array(); @@ -724,7 +724,7 @@ extern "C" Real loc[3]; //loc calculated like sum_utils.cpp - //This might be equivalent and more modular: position(i, j, k, geomdata, loc); + //This might be equivalent and more modular: position(i, j, k, geom, loc); loc[0] = problo[0] + (0.5_rt + i) * dx[0]; #if AMREX_SPACEDIM >= 2 @@ -759,13 +759,13 @@ extern "C" } void ca_derangmomy (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { int idir = 1; - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); auto const dat = datfab.array(); auto const L = derfab.array(); @@ -807,13 +807,13 @@ extern "C" } void ca_derangmomz (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { int idir = 2; - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); auto const dat = datfab.array(); auto const L = derfab.array(); @@ -856,7 +856,7 @@ extern "C" } void ca_derkineng (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -887,7 +887,7 @@ extern "C" } void ca_derspec(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -903,7 +903,7 @@ extern "C" void ca_derye(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -925,7 +925,7 @@ extern "C" } void ca_derabar(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -947,19 +947,19 @@ extern "C" } void ca_dermagvort(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); - const int coord_type = geomdata.Coord(); + const int coord_type = geom.Coord(); #if AMREX_SPACEDIM == 2 - auto problo = geomdata.ProbLoArray(); + auto problo = geom.ProbLoArray(); #endif amrex::ParallelFor(bx, @@ -1046,18 +1046,18 @@ extern "C" } void ca_derdivu(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto problo = geom.ProbLoArray(); - const int coord_type = geomdata.Coord(); + const int coord_type = geom.Coord(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -1112,7 +1112,7 @@ extern "C" } void ca_derstate(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1138,7 +1138,7 @@ extern "C" #ifdef MHD void ca_dermagcenx(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1154,7 +1154,7 @@ extern "C" } void ca_dermagceny(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1170,7 +1170,7 @@ extern "C" } void ca_dermagcenz(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1186,14 +1186,14 @@ extern "C" } void ca_derdivb(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept From 161c875ae969b327d0aef545145278f8a61ba1e0 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Wed, 13 Nov 2024 11:56:18 -0500 Subject: [PATCH 4/7] add a distance function (#2990) Adds a separate function to calculate distance, mainly because the old way gets the 2d spherical case wrong. --- Source/driver/Castro.cpp | 17 ++++++----- Source/driver/Castro_util.H | 22 ++++++++++++++ Source/driver/Derive.cpp | 49 ++++++++++++++----------------- Source/gravity/Gravity.cpp | 13 ++++---- Source/reactions/Castro_react.cpp | 10 ++++--- Source/sources/Castro_sponge.cpp | 3 +- 6 files changed, 69 insertions(+), 45 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 6149e7f988..2a47ac55dd 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -3641,19 +3641,17 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time const Real* probhi = geomdata.ProbHi(); const Real* dx = geomdata.CellSize(); - Real loc[3] = {0.0}; + GpuArray loc = {0.0}; - loc[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0]; + loc[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; #if AMREX_SPACEDIM >= 2 - loc[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1]; + loc[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; #endif #if AMREX_SPACEDIM == 3 - loc[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2]; + loc[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; #endif - Real r = std::sqrt((loc[0] - problem::center[0]) * (loc[0] - problem::center[0]) + - (loc[1] - problem::center[1]) * (loc[1] - problem::center[1]) + - (loc[2] - problem::center[2]) * (loc[2] - problem::center[2])); + Real r = distance(geomdata, loc); Real max_dist_lo = 0.0; Real max_dist_hi = 0.0; @@ -4357,9 +4355,12 @@ Castro::define_new_center(const MultiFab& S, Real time) // Now broadcast to everyone else. ParallelDescriptor::Bcast(&problem::center[0], AMREX_SPACEDIM, owner); - // Make sure if R-Z that center stays exactly on axis + // Make sure if R-Z and SPHERICAL that center stays exactly on axis if ( Geom().IsRZ() ) { problem::center[0] = 0; + } else if ( Geom().IsSPHERICAL() ) { + problem::center[0] = 0; + problem::center[1] = 0; } } diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index 016215aa32..c6162b0190 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -156,6 +156,28 @@ void position(int i, int j, int k, } + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +Real distance(GeometryData const& geomdata, GpuArray& loc) +{ + // Returns the distance from the center provided with loc, the position. + // loc is the form of [x, y, z,] in Cartesian, [r, z, phi] in cylindrical + // and [r, theta, phi] in spherical + + const auto coord = geomdata.Coord(); + + if (coord == CoordSys::cartesian) { + return std::sqrt(loc[0]*loc[0] + loc[1]*loc[1] + loc[2]*loc[2]); + } + + if (coord == CoordSys::RZ) { + return std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]); + } + + return std::abs(loc[0]); +} + + namespace geometry_util { diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index 63d31449d1..78a1dfcc35 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -579,21 +579,18 @@ extern "C" auto dx = geom.CellSizeArray(); auto problo = geom.ProbLoArray(); + auto geomdata = geom.data(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - - Real x = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; + GpuArray loc = {0.0}; + loc[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; #if AMREX_SPACEDIM >= 2 - Real y = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#else - Real y = 0.0_rt; + loc[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; #endif #if AMREX_SPACEDIM == 3 - Real z = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#else - Real z = 0.0_rt; + loc[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; #endif if (domain_is_plane_parallel) { @@ -607,15 +604,15 @@ extern "C" // where e_r and e_phi are the cylindrical unit vectors // we need the distance in the x-y plane from the origin - Real r = std::sqrt(x*x + y*y); - der(i,j,k,0) = (dat(i,j,k,1)*x + dat(i,j,k,2)*y) / (dat(i,j,k,0)*r); + Real r = std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]); + der(i,j,k,0) = (dat(i,j,k,1)*loc[0] + dat(i,j,k,2)*loc[1]) / (dat(i,j,k,0)*r); #endif } else { - Real r = std::sqrt(x*x + y*y + z*z); + Real r = distance(geomdata, loc); - der(i,j,k,0) = (dat(i,j,k,1)*x + - dat(i,j,k,2)*y + - dat(i,j,k,3)*z) / ( dat(i,j,k,0)*r ); + der(i,j,k,0) = (dat(i,j,k,1)*loc[0] + + dat(i,j,k,2)*loc[1] + + dat(i,j,k,3)*loc[2]) / ( dat(i,j,k,0)*r ); } }); @@ -634,21 +631,19 @@ extern "C" auto dx = geom.CellSizeArray(); auto problo = geom.ProbLoArray(); + auto geomdata = geom.data(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real x = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; + GpuArray loc = {0.0}; + loc[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; #if AMREX_SPACEDIM >= 2 - Real y = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#else - Real y = 0.0_rt; + loc[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; #endif #if AMREX_SPACEDIM == 3 - Real z = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#else - Real z = 0.0_rt; + loc[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; #endif if (domain_is_plane_parallel) { @@ -662,11 +657,11 @@ extern "C" // where e_r and e_phi are the cylindrical unit vectors // we need the distance in the x-y plane from the origin - Real r = std::sqrt(x*x + y*y); - der(i,j,k,0) = (-dat(i,j,k,1)*y + dat(i,j,k,2)*x) / (dat(i,j,k,0)*r); + Real r = std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]); + der(i,j,k,0) = (-dat(i,j,k,1)*loc[1] + dat(i,j,k,2)*loc[0]) / (dat(i,j,k,0)*r); #endif } else { - Real r = std::sqrt(x*x + y*y + z*z); + Real r = distance(geomdata, loc); // we really mean just the velocity component that is // perpendicular to radial, and in general 3-d (e.g. a @@ -676,9 +671,9 @@ extern "C" dat(i,j,k,2)*dat(i,j,k,2) + dat(i,j,k,3)*dat(i,j,k,3))/(dat(i,j,k,0)*dat(i,j,k,0)); - Real vr = (dat(i,j,k,1)*x + - dat(i,j,k,2)*y + - dat(i,j,k,3)*z) / ( dat(i,j,k,0)*r ); + Real vr = (dat(i,j,k,1)*loc[0] + + dat(i,j,k,2)*loc[1] + + dat(i,j,k,3)*loc[2]) / ( dat(i,j,k,0)*r ); der(i,j,k,0) = std::sqrt(amrex::max(vtot2 - vr*vr, 0.0_rt)); } diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 82003be86f..cfd40f65f0 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -1334,6 +1334,7 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab& const Real dr = dx[0] / static_cast(gravity::drdxfac); const auto problo = geom.ProbLoArray(); + const auto geomdata = geom.data(); #ifdef _OPENMP #pragma omp parallel @@ -1367,7 +1368,7 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab& loc[2] = 0.0_rt; #endif - Real r = std::sqrt(loc[0] * loc[0] + loc[1] * loc[1] + loc[2] * loc[2]); + Real r = distance(geomdata, loc); int index = static_cast(r / dr); @@ -1457,6 +1458,7 @@ Gravity::compute_radial_mass(const Box& bx, Real drinv = 1.0_rt / dr; const int coord_type = geom.Coord(); + const auto geomdata = geom.data(); AMREX_ALWAYS_ASSERT(coord_type >= 0 && coord_type <= 2); @@ -1497,16 +1499,17 @@ Gravity::compute_radial_mass(const Box& bx, amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xc = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; + GpuArray loc; + loc[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; Real lo_i = problo[0] + static_cast(i) * dx[0] - problem::center[0]; - Real yc = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; + loc[1]= problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; Real lo_j = problo[1] + static_cast(j) * dx[1] - problem::center[1]; - Real zc = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; + loc[2]= problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; Real lo_k = problo[2] + static_cast(k) * dx[2] - problem::center[2]; - Real r = std::sqrt(xc * xc + yc * yc + zc * zc); + Real r = distance(geomdata, loc); int index = static_cast(r * drinv); // We may be coming in here with a masked out zone (in a zone on a coarse diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index b7a422f1dd..9927bc4f2a 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -209,6 +209,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra const auto dx = geom.CellSizeArray(); #ifdef MODEL_PARSER const auto problo = geom.ProbLoArray(); + const auto geomdata = geom.data(); #endif #if defined(AMREX_USE_GPU) @@ -270,7 +271,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra #ifdef MODEL_PARSER if (drive_initial_convection) { - Real rr[3] = {0.0_rt}; + GpuArray rr = {0.0_rt}; rr[0] = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; #if AMREX_SPACEDIM >= 2 @@ -285,7 +286,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra if (domain_is_plane_parallel) { dist = rr[AMREX_SPACEDIM-1]; } else { - dist = std::sqrt(rr[0] * rr[0] + rr[1] * rr[1] + rr[2] * rr[2]); + dist = distance(geomdata, rr); } burn_state.T_fixed = interpolate(dist, model::itemp); @@ -565,6 +566,7 @@ Castro::react_state(Real time, Real dt) const auto dx = geom.CellSizeArray(); const auto problo = geom.ProbLoArray(); + const auto geomdata = geom.data(); #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -618,7 +620,7 @@ Castro::react_state(Real time, Real dt) #ifdef MODEL_PARSER if (drive_initial_convection) { - Real rr[3] = {0.0_rt}; + GpuArray rr = {0.0_rt}; rr[0] = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; #if AMREX_SPACEDIM >= 2 @@ -633,7 +635,7 @@ Castro::react_state(Real time, Real dt) if (domain_is_plane_parallel) { dist = rr[AMREX_SPACEDIM-1]; } else { - dist = std::sqrt(rr[0] * rr[0] + rr[1] * rr[1] + rr[2] * rr[2]); + dist = distance(geomdata, rr); } burn_state.T_fixed = interpolate(dist, model::itemp); diff --git a/Source/sources/Castro_sponge.cpp b/Source/sources/Castro_sponge.cpp index 2109d63374..06781f4403 100644 --- a/Source/sources/Castro_sponge.cpp +++ b/Source/sources/Castro_sponge.cpp @@ -81,6 +81,7 @@ Castro::apply_sponge(const Box& bx, auto dx = geom.CellSizeArray(); auto problo = geom.ProbLoArray(); + auto geomdata = geom.data(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -124,7 +125,7 @@ Castro::apply_sponge(const Box& bx, Real sponge_factor = 0.0_rt; if (sponge_lower_radius >= 0.0_rt && sponge_upper_radius > sponge_lower_radius) { - Real rad = std::sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]); + Real rad = distance(geomdata, r); if (rad < sponge_lower_radius) { sponge_factor = sponge_lower_factor; From aed0328358bb018bb16aaf5dd6a79cfa250bfc89 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Thu, 14 Nov 2024 08:04:34 -0500 Subject: [PATCH 5/7] xrb spherical problem setup (#2972) Initial problem setup for full-star xrb flame in spherical shell based on the flame_wave problem. --- Exec/science/xrb_spherical/GNUmakefile | 40 ++++ Exec/science/xrb_spherical/Make.package | 2 + Exec/science/xrb_spherical/README.md | 6 + Exec/science/xrb_spherical/_prob_params | 68 ++++++ .../xrb_spherical/analysis/r_profile.py | 56 +++++ Exec/science/xrb_spherical/analysis/slice.py | 94 ++++++++ Exec/science/xrb_spherical/initial_model.H | 1 + Exec/science/xrb_spherical/inputs.He.1000Hz | 147 ++++++++++++ .../xrb_spherical/problem_initialize.H | 222 ++++++++++++++++++ .../problem_initialize_state_data.H | 78 ++++++ Exec/science/xrb_spherical/problem_tagging.H | 56 +++++ 11 files changed, 770 insertions(+) create mode 100644 Exec/science/xrb_spherical/GNUmakefile create mode 100644 Exec/science/xrb_spherical/Make.package create mode 100644 Exec/science/xrb_spherical/README.md create mode 100644 Exec/science/xrb_spherical/_prob_params create mode 100755 Exec/science/xrb_spherical/analysis/r_profile.py create mode 100755 Exec/science/xrb_spherical/analysis/slice.py create mode 120000 Exec/science/xrb_spherical/initial_model.H create mode 100644 Exec/science/xrb_spherical/inputs.He.1000Hz create mode 100644 Exec/science/xrb_spherical/problem_initialize.H create mode 100644 Exec/science/xrb_spherical/problem_initialize_state_data.H create mode 100644 Exec/science/xrb_spherical/problem_tagging.H diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile new file mode 100644 index 0000000000..d53475be8d --- /dev/null +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -0,0 +1,40 @@ +PRECISION = DOUBLE +PROFILE = FALSE + +DEBUG = FALSE + +DIM = 2 + +COMP = gnu + +USE_MPI = TRUE + +USE_GRAV = TRUE +USE_REACT = FALSE + +USE_ROTATION = FALSE +USE_DIFFUSION = FALSE + +# define the location of the CASTRO top directory +CASTRO_HOME ?= ../../.. + +USE_JACOBIAN_CACHING = TRUE +USE_MODEL_PARSER = TRUE +NUM_MODELS := 2 + +# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos +EOS_DIR := helmholtz + +# This sets the network directory in $(MICROPHYSICS_HOME)/networks +NETWORK_DIR := subch_base + +INTEGRATOR_DIR := VODE + +CONDUCTIVITY_DIR := stellar + +PROBLEM_DIR ?= ./ + +Bpack := $(PROBLEM_DIR)/Make.package +Blocs := $(PROBLEM_DIR) + +include $(CASTRO_HOME)/Exec/Make.Castro diff --git a/Exec/science/xrb_spherical/Make.package b/Exec/science/xrb_spherical/Make.package new file mode 100644 index 0000000000..e5cc052427 --- /dev/null +++ b/Exec/science/xrb_spherical/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += initial_model.H + diff --git a/Exec/science/xrb_spherical/README.md b/Exec/science/xrb_spherical/README.md new file mode 100644 index 0000000000..325667c9c1 --- /dev/null +++ b/Exec/science/xrb_spherical/README.md @@ -0,0 +1,6 @@ +# xrb_spherical + +This is the full-star XRB flame setup based on flame_wave. +This setup uses a spherical 2D geometry to model XRB flame +on a spherical shell with initial temperature perturbation +on the north pole. \ No newline at end of file diff --git a/Exec/science/xrb_spherical/_prob_params b/Exec/science/xrb_spherical/_prob_params new file mode 100644 index 0000000000..89c1775615 --- /dev/null +++ b/Exec/science/xrb_spherical/_prob_params @@ -0,0 +1,68 @@ + +dtemp real 3.81e8_rt y + +theta_half_max real 1.745e-2_rt y + +theta_half_width real 4.9e-3_rt y + +# cutoff mass fraction of the first species for refinement +X_min real 1.e-4_rt y + +# do we dynamically refine based on density? or based on height? +tag_by_density integer 1 y + +# used for tagging if tag_by_density = 1 +cutoff_density real 500.e0_rt y + +# used if we are refining based on height rather than density +refine_height real 3600 y + +T_hi real 5.e8_rt y + +T_star real 1.e8_rt y + +T_lo real 5.e7_rt y + +dens_base real 2.e6_rt y + +H_star real 500.e0_rt y + +atm_delta real 25.e0_rt y + +fuel1_name string "helium-4" y + +fuel2_name string "" y + +fuel3_name string "" y + +fuel4_name string "" y + +ash1_name string "iron-56" y + +ash2_name string "" y + +ash3_name string "" y + +fuel1_frac real 1.0_rt y + +fuel2_frac real 0.0_rt y + +fuel3_frac real 0.0_rt y + +fuel4_frac real 0.0_rt y + +ash1_frac real 1.0_rt y + +ash2_frac real 0.0_rt y + +ash3_frac real 0.0_rt y + +low_density_cutoff real 1.e-4_rt y + +smallx real 1.e-10_rt y + +r_refine_distance real 1.e30_rt y + +max_hse_tagging_level integer 2 y + +max_base_tagging_level integer 1 y diff --git a/Exec/science/xrb_spherical/analysis/r_profile.py b/Exec/science/xrb_spherical/analysis/r_profile.py new file mode 100755 index 0000000000..468f2ba945 --- /dev/null +++ b/Exec/science/xrb_spherical/analysis/r_profile.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 + +# Spherical R profile at different theta + +import os +import sys +import yt +import matplotlib.pyplot as plt +import numpy as np +from functools import reduce +import itertools + +import matplotlib.ticker as ptick +from yt.frontends.boxlib.api import CastroDataset +from yt.units import cm + + +plotfile = sys.argv[1] +ds = CastroDataset(plotfile) + +rmin = ds.domain_left_edge[0] +rmax = rmin + 5000.0*cm +#rmax = ds.domain_right_edge[0] +print(ds.domain_left_edge[1]) +fig, _ax = plt.subplots(2,2) + +axes = list(itertools.chain(*_ax)) + +fig.set_size_inches(7.0, 8.0) + +fields = ["Temp", "density", "x_velocity", "y_velocity"] +nice_names = [r"$T$ (K)", r"$\rho$ (g/${cm}^3$)", r"$u$ (cm/s)", r"$v$ (cm/s)"] + +# 4 rays at different theta values +thetal = ds.domain_left_edge[1] +thetar = ds.domain_right_edge[1] +thetas = [thetal, 0.25*thetar, 0.5*thetar, 0.75*thetar] + +for i, f in enumerate(fields): + + for theta in thetas: + # simply go from (rmin, theta) -> (rmax, theta). Doesn't need to convert to physical R-Z + ray = ds.ray((rmin, theta, 0*cm), (rmax, theta, 0*cm)) + isrt = np.argsort(ray["t"]) + axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(float(theta))) + + axes[i].set_xlabel(r"$r$ (cm)") + axes[i].set_ylabel(nice_names[i]) + axes[i].set_yscale("symlog") + + if i == 0: + axes[0].legend(frameon=False, loc="lower left") + +#fig.set_size_inches(10.0, 9.0) +plt.tight_layout() +plt.savefig("{}_profiles.png".format(os.path.basename(plotfile))) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py new file mode 100755 index 0000000000..7e55d3848f --- /dev/null +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python3 + +import sys +import os +import yt +import numpy as np +import matplotlib.pyplot as plt +from yt.frontends.boxlib.api import CastroDataset + +from yt.units import cm + +""" +Given a plot file and field name, it gives slice plots at the top, +middle, and bottom of the domain (shell). +""" + +def slice(fname:str, field:str, + loc: str = "top", width_factor: float = 3.0) -> None: + """ + A slice plot of the dataset for Spherical2D geometry. + + Parameter + ======================= + fname: plot file name + field: field parameter + loc: location on the domain. {top, mid, bot} + """ + + ds = CastroDataset(fname) + currentTime = ds.current_time.in_units("s") + print(f"Current time of this plot file is {currentTime} s") + + # Some geometry properties + rr = ds.domain_right_edge[0].in_units("km") + rl = ds.domain_left_edge[0].in_units("km") + dr = rr - rl + r_center = 0.5 * (rr + rl) + + thetar = ds.domain_right_edge[1] + thetal = ds.domain_left_edge[1] + dtheta = thetar - thetal + theta_center = 0.5 * (thetar + thetal) + + # Domain width of the slice plot + width = width_factor * dr + box_widths = (width, width) + + loc = loc.lower() + loc_options = ["top", "mid", "bot"] + + if loc not in loc_options: + raise Exception("loc parameter must be top, mid or bot") + + # Centers for the Top, Mid and Bot panels + centers = {"top":(r_center*np.sin(thetal)+0.5*width, r_center*np.cos(thetal)), + "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)), + "bot":(r_center*np.sin(thetar)+0.5*width, r_center*np.cos(thetar))} + + # Note we can also set center during SlicePlot, however then we would enter in [r_center, theta_center, 0] + # rather than the physical R-Z coordinate if we do it via sp.set_center + + sp = yt.SlicePlot(ds, 'phi', field, width=box_widths) + sp.set_center(centers[loc]) + + sp.set_cmap(field, "viridis") + if field in ["x_velocity", "y_velocity", "z_velocity"]: + sp.set_cmap(field, "coolwarm") + elif field == "Temp": + sp.set_zlim(f, 5.e7, 2.5e9) + sp.set_cmap(f, "magma_r") + elif field == "enuc": + sp.set_zlim(f, 1.e18, 1.e20) + elif field == "density": + sp.set_zlim(f, 1.e-3, 5.e8) + + # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s") + sp.save(f"{ds}_{loc}") + +if __name__ == "__main__": + + if len(sys.argv) < 3: + raise Exception("Please enter parameters in order of: fname field_name width_factor[optional] loc[optional]") + + fname = sys.argv[1] + field = sys.argv[2] + loc = "top" + width_factor = 3.0 + + if len(sys.argv) == 4: + width_factor = float(sys.argv[3]) + elif len(sys.argv) > 4: + loc = sys.argv[4] + + slice(fname, field, loc=loc, width_factor=width_factor) diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H new file mode 120000 index 0000000000..9c923c3113 --- /dev/null +++ b/Exec/science/xrb_spherical/initial_model.H @@ -0,0 +1 @@ +../flame_wave/initial_model.H \ No newline at end of file diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz new file mode 100644 index 0000000000..17864ed3f3 --- /dev/null +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -0,0 +1,147 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 9900000 +stop_time = 3.0 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 +geometry.coord_sys = 2 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = 1.1e6 0 +geometry.prob_hi = 1.13072e6 3.1415926 +amr.n_cell = 768 2304 #192 1152 + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +castro.lo_bc = 3 3 # Inflow in lower R and Symmetry across Theta +castro.hi_bc = 2 3 # Outflow in upper boundaries + +# Allow non-square zones +castro.allow_non_unit_aspect_zones = 1 + +# Fill ambient states with outflow velocity in R-direction +castro.fill_ambient_bc = 1 +castro.ambient_fill_dir = 0 +castro.ambient_outflow_vel = 1 + +# WHICH PHYSICS +castro.do_hydro = 1 +castro.do_react = 1 +castro.do_rotation = 1 +castro.do_grav = 1 +castro.do_sponge = 1 + +castro.small_temp = 1.e6 +castro.small_dens = 1.e-5 + +castro.ppm_type = 1 +castro.grav_source_type = 2 +castro.use_pslope = 1 +castro.ppm_well_balanced = 1 +castro.pslope_cutoff_density = 1.e4 + +gravity.gravity_type = ConstantGrav + +# 1.4 Solar Mass NS with radius ~11 km +gravity.const_grav = -1.5e14 + +# 1000Hz Spinning Frequency +# Might want to use a more realistic spinning frequency like 500Hz +castro.rotational_period = 0.001 + +# Centrifugal is not important since NS would simply deform to accommodate for it +castro.rotation_include_centrifugal = 0 + +castro.diffuse_temp = 1 +castro.diffuse_cutoff_density_hi = 5.e4 +castro.diffuse_cutoff_density = 2.e4 + +castro.diffuse_cond_scale_fac = 1.0 + +castro.react_rho_min = 1.e2 +castro.react_rho_max = 1.5e7 + +castro.react_T_min = 6.e7 + +castro.sponge_upper_density = 1.e2 +castro.sponge_lower_density = 1.e0 +castro.sponge_timescale = 1.e-7 + +castro.abundance_failure_tolerance = 0.1 +castro.abundance_failure_rho_cutoff = 1.0 + +# TIME STEP CONTROL +castro.cfl = 0.8 # cfl number for hyperbolic system +castro.init_shrink = 0.1 # scale back initial timestep +castro.change_max = 1.1 # max time step growth + +castro.use_retry = 1 +castro.max_subcycles = 16 + +castro.retry_small_density_cutoff = 1.0 + +# DIAGNOSTICS & VERBOSITY +castro.sum_interval = 0 # timesteps between computing mass +castro.v = 1 # verbosity in Castro.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 #2 # maximum level number allowed +amr.ref_ratio = 4 2 2 2 # refinement ratio +amr.regrid_int = 2 2 2 2 # how often to regrid +amr.blocking_factor = 16 # block factor in grid generation +amr.max_grid_size = 128 +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.check_file = flame_wave_1000Hz_chk # root name of checkpoint file +amr.check_int = 250 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_file = flame_wave_1000Hz_plt # root name of plotfile +amr.plot_per = 5.e-3 # number of seconds between plotfiles +amr.derive_plot_vars = ALL + +amr.small_plot_file = flame_wave_1000Hz_smallplt # root name of plotfile +amr.small_plot_per = 1.e-6 #1.e-4 # number of seconds between plotfiles +amr.small_plot_vars = density Temp +amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc + +# problem initialization + +problem.dtemp = 1.2e9 +problem.theta_half_max = 1.745e-2 +problem.theta_half_width = 5.279e-3 + +problem.dens_base = 3.43e6 + +problem.T_star = 1.e8 +problem.T_hi = 2.e8 +problem.T_lo = 8.e6 + +problem.H_star = 2000.e0 +problem.atm_delta = 400.0 #50.0 + +problem.fuel1_name = "helium-4" +problem.fuel1_frac = 1.0e0 + +problem.ash1_name = "nickel-56" +problem.ash1_frac = 1.0e0 + +problem.low_density_cutoff = 1.e-4 + +problem.cutoff_density = 2.5e4 +problem.max_hse_tagging_level = 3 +problem.max_base_tagging_level = 2 + +problem.X_min = 1.e-2 + +problem.r_refine_distance = 9.216e4 + +# Microphysics + +integrator.rtol_spec = 1.e-6 +integrator.atol_spec = 1.e-6 + +network.use_tables = 1 diff --git a/Exec/science/xrb_spherical/problem_initialize.H b/Exec/science/xrb_spherical/problem_initialize.H new file mode 100644 index 0000000000..cfb11741ff --- /dev/null +++ b/Exec/science/xrb_spherical/problem_initialize.H @@ -0,0 +1,222 @@ +#ifndef problem_initialize_H +#define problem_initialize_H + +#include +#include +#include +#include +#include +#include +#include + +AMREX_INLINE +void problem_initialize () +{ + + const Geometry& dgeom = DefaultGeometry(); + + const Real* problo = dgeom.ProbLo(); + const Real* probhi = dgeom.ProbHi(); + + // check to make sure that small_dens is less than low_density_cutoff + // if not, funny things can happen above the atmosphere + + if (small_dens >= 0.99_rt * problem::low_density_cutoff) { + amrex::Error("ERROR: small_dens should be set lower than low_density_cutoff"); + } + + // make sure hse_fixed_temp is the same as T_star, if it's specified + + if (hse_fixed_temp > 0.0_rt && hse_fixed_temp != problem::T_star) { + amrex::Error("ERROR: hse_fixed_temp should be the same as T_star"); + } + + // get the species indices + + bool species_defined = true; + + int ifuel1 = network_spec_index(problem::fuel1_name); + if (ifuel1 < 0) { + species_defined = false; + } + + int ifuel2; + if (!problem::fuel2_name.empty()) { + ifuel2 = network_spec_index(problem::fuel2_name); + if (ifuel2 < 0) { + species_defined = false; + } + } + + int ifuel3; + if (!problem::fuel3_name.empty()) { + ifuel3 = network_spec_index(problem::fuel3_name); + if (ifuel3 < 0) { + species_defined = false; + } + } + + int ifuel4; + if (!problem::fuel4_name.empty()) { + ifuel4 = network_spec_index(problem::fuel4_name); + if (ifuel4 < 0) { + species_defined = false; + } + } + + int iash1 = network_spec_index(problem::ash1_name); + if (iash1 < 0) { + species_defined = false; + } + + int iash2; + if (!problem::ash2_name.empty()) { + iash2 = network_spec_index(problem::ash2_name); + if (iash2 < 0) { + species_defined = false; + } + } + + int iash3; + if (!problem::ash3_name.empty()) { + iash3 = network_spec_index(problem::ash3_name); + if (iash3 < 0) { + species_defined = false; + } + } + + if (! species_defined) { + std::cout << ifuel1 << " " << ifuel2 << " " << ifuel3 << " " << ifuel4 << std::endl; + std::cout << iash1 << " " << iash2 << " "<< iash3 << std::endl; + amrex::Error("ERROR: species not defined"); + } + + model_t model_params; + + // set the composition of the underlying star + + + for (Real &X : model_params.xn_star) { + X = problem::smallx; + } + model_params.xn_star[iash1] = problem::ash1_frac; + if (!problem::ash2_name.empty()) { + model_params.xn_star[iash2] = problem::ash2_frac; + } + if (!problem::ash3_name.empty()) { + model_params.xn_star[iash3] = problem::ash3_frac; + } + + // and the composition of the accreted layer + + for (Real &X : model_params.xn_base) { + X = problem::smallx; + } + model_params.xn_base[ifuel1] = problem::fuel1_frac; + if (!problem::fuel2_name.empty()) { + model_params.xn_base[ifuel2] = problem::fuel2_frac; + } + if (!problem::fuel3_name.empty()) { + model_params.xn_base[ifuel3] = problem::fuel3_frac; + } + if (!problem::fuel4_name.empty()) { + model_params.xn_base[ifuel4] = problem::fuel4_frac; + } + + // check if they sum to 1 + + Real sumX = 0.0_rt; + for (Real X : model_params.xn_star) { + sumX += X; + } + if (std::abs(sumX - 1.0_rt) > NumSpec * problem::smallx) { + amrex::Error("ERROR: ash mass fractions don't sum to 1"); + } + + sumX = 0.0_rt; + for (Real X : model_params.xn_base) { + sumX += X; + } + if (std::abs(sumX - 1.0_rt) > NumSpec * problem::smallx) { + amrex::Error("ERROR: fuel mass fractions don't sum to 1"); + } + + // we are going to generate an initial model from probl(1), i.e. r_min, to + // probhi(1), i.e., r_max, with nx_model zones. But to allow for a interpolated + // lower boundary, we'll add 4 ghostcells to this, so we need to + // compute dx + + // we use the fine grid dx for the model resolution + auto fine_geom = global::the_amr_ptr->Geom(global::the_amr_ptr->maxLevel()); + + auto dx = fine_geom.CellSizeArray(); + auto dx_model = dx[0]; + + int nx_model = static_cast((probhi[0] - problo[0]) / + dx_model); + + int ng = 4; + + // now generate the initial models + + model_params.dens_base = problem::dens_base; + model_params.T_star = problem::T_star; + model_params.T_hi = problem::T_hi; + model_params.T_lo = problem::T_lo; + + model_params.H_star = problem::H_star; + model_params.atm_delta = problem::atm_delta; + + model_params.low_density_cutoff = problem::low_density_cutoff; + + generate_initial_model(nx_model + ng, + problo[0] - ng * dx_model, + probhi[0], + model_params, 0); + + // now create a perturbed model -- we want the same base conditions + // a hotter temperature + + model_params.T_hi = model_params.T_hi + problem::dtemp; + + generate_initial_model(nx_model + ng, + problo[0] - ng * dx_model, + probhi[0], + model_params, 1); + + // set center + + for (int d = 0; d < AMREX_SPACEDIM; d++) { + // problem::center[d] = 0.5_rt * (problo[d] + probhi[d]); + problem::center[d] = 0.0_rt; + } + + // set the ambient state for the upper boundary condition + + ambient::ambient_state[URHO] = model::profile(0).state(model::npts-1, model::idens); + ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts-1, model::itemp); + for (int n = 0; n < NumSpec; n++) { + ambient::ambient_state[UFS+n] = + ambient::ambient_state[URHO] * model::profile(0).state(model::npts-1, model::ispec+n); + } + + ambient::ambient_state[UMX] = 0.0_rt; + ambient::ambient_state[UMY] = 0.0_rt; + ambient::ambient_state[UMZ] = 0.0_rt; + + // make the ambient state thermodynamically consistent + + eos_t eos_state; + eos_state.rho = ambient::ambient_state[URHO]; + eos_state.T = ambient::ambient_state[UTEMP]; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = ambient::ambient_state[UFS+n] / eos_state.rho; + } + + eos(eos_input_rt, eos_state); + + ambient::ambient_state[UEINT] = eos_state.rho * eos_state.e; + ambient::ambient_state[UEDEN] = eos_state.rho * eos_state.e; +} + +#endif diff --git a/Exec/science/xrb_spherical/problem_initialize_state_data.H b/Exec/science/xrb_spherical/problem_initialize_state_data.H new file mode 100644 index 0000000000..0431916a8e --- /dev/null +++ b/Exec/science/xrb_spherical/problem_initialize_state_data.H @@ -0,0 +1,78 @@ +#ifndef problem_initialize_state_data_H +#define problem_initialize_state_data_H + +#include +#include +#include +#include +#include +#include + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_initialize_state_data (int i, int j, int k, + Array4 const& state, + const GeometryData& geomdata) +{ + + const Real* dx = geomdata.CellSize(); + const Real* problo = geomdata.ProbLo(); + + Real r = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; + Real theta = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; + + // blending factor + + Real f; + + if (theta < problem::theta_half_max) { + f = 1.0_rt; + + } else if (theta > problem::theta_half_max + problem::theta_half_width) { + f = 0.0_rt; + + } else { + f = -(theta - problem::theta_half_max) / problem::theta_half_width + 1.0_rt; + } + + state(i,j,k,URHO) = f * interpolate(r, model::idens, 1) + + (1.0_rt - f) * interpolate(r, model::idens, 0); + + state(i,j,k,UTEMP) = f * interpolate(r, model::itemp, 1) + + (1.0_rt - f) * interpolate(r, model::itemp, 0); + + Real temppres = f * interpolate(r, model::ipres, 1) + + (1.0_rt - f) * interpolate(r, model::ipres, 0); + + for (int n = 0; n < NumSpec; n++) { + state(i,j,k,UFS+n) = f * interpolate(r, model::ispec+n, 1) + + (1.0_rt - f) * interpolate(r, model::ispec+n, 0); + } + + eos_t eos_state; + eos_state.rho = state(i,j,k,URHO); + eos_state.T = state(i,j,k,UTEMP); + eos_state.p = temppres; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = state(i,j,k,UFS+n); + } + + eos(eos_input_rp, eos_state); + + state(i,j,k,UTEMP) = eos_state.T; + state(i,j,k,UEINT) = eos_state.rho * eos_state.e; + state(i,j,k,UEDEN) = state(i,j,k,UEINT); + + // Initial velocities = 0 + + state(i,j,k,UMX) = 0.e0_rt; + state(i,j,k,UMY) = 0.e0_rt; + state(i,j,k,UMZ) = 0.e0_rt; + + // convert to partial densities + + for (int n = 0; n < NumSpec; n++) { + state(i,j,k,UFS+n) = state(i,j,k,URHO) * state(i,j,k,UFS+n); + } +} + +#endif diff --git a/Exec/science/xrb_spherical/problem_tagging.H b/Exec/science/xrb_spherical/problem_tagging.H new file mode 100644 index 0000000000..48b2f254c6 --- /dev/null +++ b/Exec/science/xrb_spherical/problem_tagging.H @@ -0,0 +1,56 @@ +#ifndef problem_tagging_H +#define problem_tagging_H +#include +#include +#include + +/// +/// Define problem-specific tagging criteria +/// +/// @param i x-index +/// @param j y-index +/// @param k z-index +/// @param tag tag array (TagBox) +/// @param state simulation state (Fab) +/// @param level AMR level +/// @param geomdata geometry data +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_tagging(int i, int j, int k, + Array4 const& tag, + Array4 const& state, + int level, const GeometryData& geomdata) +{ + + GpuArray loc; + position(i, j, k, geomdata, loc); + + if (problem::tag_by_density) { + if (state(i,j,k,URHO) > problem::cutoff_density && + state(i,j,k,UFS) / state(i,j,k,URHO) > problem::X_min) { + + Real dist = std::abs(loc[0]); + + if (level < problem::max_hse_tagging_level && dist < geomdata.ProbLo(0) + problem::r_refine_distance) { + tag(i,j,k) = TagBox::SET; + } + } + + if (state(i,j,k,URHO) > problem::cutoff_density) { + if (level < problem::max_base_tagging_level) { + tag(i,j,k) = TagBox::SET; + } + } + + } else { + + // tag everything below a certain height + if (loc[0] < geomdata.ProbLo(0) + problem::refine_height) { + if (level < problem::max_base_tagging_level) { + tag(i,j,k) = TagBox::SET; + } + } + } + +} +#endif From be7843093e870734e219dcf73df2094ab4abece7 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Sun, 17 Nov 2024 14:40:02 -0500 Subject: [PATCH 6/7] remove extra pressure scale in scale_flux for 1d cartesian (#2991) addresses #2803 . In #2468, we modified mom_flux_has_p to include pressure in momentum flux with 1d cartesian, so we don't need to scale pressure again scale_flux. --- .../ci-benchmarks/Rad2TShock-1d.out | 52 +++++++++---------- Source/hydro/Castro_ctu_hydro.cpp | 29 ++++------- Source/hydro/Castro_hydro.H | 3 -- Source/hydro/Castro_mol_hydro.cpp | 31 ++++------- Source/hydro/advection_util.cpp | 14 ----- 5 files changed, 46 insertions(+), 83 deletions(-) diff --git a/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out b/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out index c4b21ef48f..0e1151359f 100644 --- a/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out +++ b/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out @@ -1,40 +1,40 @@ plotfile = plt_00025 - time = 0.00026083835470476599 + time = 0.00026084061812361401 variables minimum value maximum value - density 5.4596904436e-13 1.2720098871e-12 - xmom 1.2327065163e-07 1.4587103381e-07 + density 5.4596904436e-13 1.2717694091e-12 + xmom 1.2328047524e-07 1.4591226172e-07 ymom 0 0 zmom 0 0 - rho_E 0.021940622286 0.039199658391 - rho_e 0.0068091599527 0.032333529333 - Temp 100.00003116 217.3564598 - rho_X 5.4596904436e-13 1.2720098871e-12 - rad 7.5662679486e-07 1.4083425328e-05 - pressure 0.0045394399687 0.021555686223 - kineng 0.0064467244856 0.015131462333 - soundspeed 117717.6283 173551.23637 + rho_E 0.021940622285 0.039194728059 + rho_e 0.0068091599514 0.032333527949 + Temp 100.00003115 217.29660609 + rho_X 5.4596904436e-13 1.2717694091e-12 + rad 7.5662675907e-07 1.4083418723e-05 + pressure 0.0045394399678 0.0215556853 + kineng 0.0064473434191 0.015131462334 + soundspeed 117717.62829 173527.33922 Gamma_1 1.6666666667 1.6666666667 - MachNumber 0.60689740202 1.9999997205 - uplusc 269483.97988 362284.15407 - uminusc -66810.704217 117717.5954 - entropy 2458725989.8 2507243428 + MachNumber 0.6068973935 1.9999997207 + uplusc 269518.17809 362247.33487 + uminusc -66804.446805 117717.59541 + entropy 2458725989.8 2507184972.2 magvort 0 0 - divu -20286.061564 245.50876324 - eint_E 12471696011 27108028479 - eint_e 12471696011 27108028479 - logden -12.26283198 -11.895509513 - StateErr_0 5.4596904436e-13 1.2720098871e-12 - StateErr_1 100.00003116 217.3564598 + divu -20283.836079 245.94831837 + eint_E 12471696009 27100563708 + eint_e 12471696009 27100563708 + logden -12.26283198 -11.895591626 + StateErr_0 5.4596904436e-13 1.2717694091e-12 + StateErr_1 100.00003115 217.29660609 StateErr_2 1 1 X(X) 1 1 abar 1 1 - x_velocity 102258.95554 235435.2237 + x_velocity 102299.57948 235435.22371 y_velocity 0 0 z_velocity 0 0 - magvel 102258.95554 235435.2237 - radvel -235435.2237 114038.34336 - circvel 0 0.005524271728 - magmom 1.2327065163e-07 1.4587103381e-07 + magvel 102299.57948 235435.22371 + radvel -235435.22371 114042.93363 + circvel 0 0.0047841596539 + magmom 1.2328047524e-07 1.4591226172e-07 angular_momentum_x 0 0 angular_momentum_y 0 0 angular_momentum_z -0 -0 diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index d443bbe060..e982187f92 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -1259,35 +1259,26 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co Array4 const flux_arr = (flux[idir]).array(); Array4 const area_arr = (area[idir]).array(mfi); - scale_flux(nbx, -#if AMREX_SPACEDIM == 1 - qex_arr, -#endif - flux_arr, area_arr, dt); + scale_flux(nbx, flux_arr, area_arr, dt); #ifdef RADIATION Array4 const rad_flux_arr = (rad_flux[idir]).array(); scale_rad_flux(nbx, rad_flux_arr, area_arr, dt); #endif - if (idir == 0) { #if AMREX_SPACEDIM <= 2 - Array4 pradial_fab = pradial.array(); -#endif + // get the scaled radial pressure -- we need to treat this specially - // get the scaled radial pressure -- we need to treat this specially -#if AMREX_SPACEDIM <= 2 - if (!mom_flux_has_p(0, 0, coord)) { - amrex::ParallelFor(nbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; - }); - } + if (idir == 0 && !mom_flux_has_p(0, 0, coord)) { + Array4 pradial_fab = pradial.array(); -#endif + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; + }); } - +#endif // Store the fluxes from this advance. For simplified SDC integration we // only need to do this on the last iteration. diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index aa8c9582c5..cc7cd154cb 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -784,9 +784,6 @@ void scale_flux(const amrex::Box& bx, -#if AMREX_SPACEDIM == 1 - amrex::Array4 const& qint, -#endif amrex::Array4 const& flux, amrex::Array4 const& area, const amrex::Real dt); diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 75899d67a1..8b00cd958b 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -657,9 +657,6 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) Array4 pradial_fab = pradial.array(); #endif -#if AMREX_SPACEDIM == 1 - Array4 const qex_arr = qe[0].array(); -#endif for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) { @@ -668,27 +665,19 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) Array4 const flux_arr = (flux[idir]).array(); Array4 const area_arr = (area[idir]).array(mfi); - scale_flux(nbx, -#if AMREX_SPACEDIM == 1 - qex_arr, -#endif - flux_arr, area_arr, dt); - + scale_flux(nbx, flux_arr, area_arr, dt); - if (idir == 0) { - // get the scaled radial pressure -- we need to treat this specially - Array4 const qex_fab = qe[idir].array(); - const int prescomp = GDPRES; +#if AMREX_SPACEDIM <= 2 + // get the scaled radial pressure -- we need to treat this specially + if (idir == 0 && !mom_flux_has_p(0, 0, coord)) { + Array4 const qex_arr = qe[idir].array(); -#if AMREX_SPACEDIM <= 2 - if (!mom_flux_has_p(0, 0, coord)) { - amrex::ParallelFor(nbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pradial_fab(i,j,k) = qex_fab(i,j,k,prescomp) * dt; - }); - } + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; + }); #endif } } diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 515778c4d1..05ce02e699 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -503,28 +503,14 @@ Castro::normalize_species_fluxes(const Box& bx, void // NOLINTNEXTLINE(readability-convert-member-functions-to-static) Castro::scale_flux(const Box& bx, -#if AMREX_SPACEDIM == 1 - Array4 const& qint, -#endif Array4 const& flux, Array4 const& area_arr, const Real dt) { -#if AMREX_SPACEDIM == 1 - const int coord_type = geom.Coord(); -#endif - amrex::ParallelFor(bx, NUM_STATE, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - flux(i,j,k,n) = dt * flux(i,j,k,n) * area_arr(i,j,k); -#if AMREX_SPACEDIM == 1 - // Correct the momentum flux with the grad p part. - if (coord_type == 0 && n == UMX) { - flux(i,j,k,n) += dt * area_arr(i,j,k) * qint(i,j,k,GDPRES); - } -#endif }); } From 2fb7228782adc65abef7236f6dfa72bab0f58b45 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Sun, 17 Nov 2024 14:43:21 -0500 Subject: [PATCH 7/7] pass in the correct dloga when doing tracing (#2994) computes dloga in y-direction and pass in the correct dloga for the corresponding dir. This was missed in pr #2964 . --- Source/driver/Castro.H | 5 +++-- Source/driver/Castro.cpp | 13 +++++++++++-- Source/hydro/Castro_ctu.cpp | 27 ++++++++++++++++++--------- Source/hydro/Castro_ctu_hydro.cpp | 22 +++++++++++++++++----- Source/hydro/Castro_hydro.H | 24 ++++++++++++++++++------ 5 files changed, 67 insertions(+), 24 deletions(-) diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index a3c6983c4a..0a77218ee6 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -1381,9 +1381,10 @@ protected: /// amrex::MultiFab volume; amrex::Array area; - amrex::MultiFab dLogArea[1]; amrex::Vector< amrex::Vector > radius; - +#if AMREX_SPACEDIM <= 2 + amrex::Array dLogArea; +#endif /// /// For initialization of the C++ values of the runtime parameters diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 2a47ac55dd..0133c68bfe 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -845,9 +845,18 @@ Castro::buildMetrics () area[dir].setVal(0.0); } - dLogArea[0].clear(); #if (AMREX_SPACEDIM <= 2) - geom.GetDLogA(dLogArea[0],grids,dmap,0,NUM_GROW); + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) + { + dLogArea[dir].clear(); + geom.GetDLogA(dLogArea[dir],grids,dmap,dir,NUM_GROW); + } + for (int dir = AMREX_SPACEDIM; dir < 2; dir++) + { + dLogArea[dir].clear(); + dLogArea[dir].define(grids, dmap, 1, 0); + dLogArea[dir].setVal(0.0); + } #endif wall_time_start = 0.0; diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 6c222cdddf..18756b0b71 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -101,7 +101,10 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, Array4 const& qzp, #endif #if AMREX_SPACEDIM < 3 - Array4 const& dloga, + Array4 const& dlogaX, +#endif +#if AMREX_SPACEDIM == 2 + Array4 const& dlogaY, #endif const Real dt) { @@ -118,7 +121,7 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ, qxm, qxp, #if AMREX_SPACEDIM <= 2 - dloga, + dlogaX, #endif vbx, dt); @@ -131,7 +134,7 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ, qym, qyp, #if AMREX_SPACEDIM <= 2 - dloga, + dlogaY, #endif vbx, dt); @@ -173,7 +176,10 @@ Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx, Array4 const& qzp, #endif #if AMREX_SPACEDIM < 3 - Array4 const& dloga, + Array4 const& dlogaX, +#endif +#if AMREX_SPACEDIM == 2 + Array4 const& dlogaY, #endif const Real dt) { @@ -187,7 +193,7 @@ Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx, U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ, qxm, qxp, #if AMREX_SPACEDIM <= 2 - dloga, + dlogaX, #endif vbx, dt); @@ -200,7 +206,7 @@ Castro::ctu_ppm_rad_states(const Box& bx, const Box& vbx, U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ, qym, qyp, #if AMREX_SPACEDIM <= 2 - dloga, + dlogaY, #endif vbx, dt); @@ -243,7 +249,10 @@ Castro::ctu_plm_states(const Box& bx, const Box& vbx, Array4 const& qzp, #endif #if AMREX_SPACEDIM < 3 - Array4 const& dloga, + Array4 const& dlogaX, +#endif +#if AMREX_SPACEDIM == 2 + Array4 const& dlogaY, #endif const Real dt) { @@ -273,7 +282,7 @@ Castro::ctu_plm_states(const Box& bx, const Box& vbx, U_arr, rho_inv_arr, q_arr, qaux_arr, qxm, qxp, #if AMREX_SPACEDIM < 3 - dloga, + dlogaX, #endif srcQ, vbx, dt); @@ -285,7 +294,7 @@ Castro::ctu_plm_states(const Box& bx, const Box& vbx, U_arr, rho_inv_arr, q_arr, qaux_arr, qym, qyp, #if AMREX_SPACEDIM < 3 - dloga, + dlogaY, #endif srcQ, vbx, dt); diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index e982187f92..9791280851 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -242,7 +242,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co #endif #if AMREX_SPACEDIM < 3 - Array4 const dLogArea_arr = (dLogArea[0]).array(mfi); + Array4 const dLogAreaX_arr = (dLogArea[0]).array(mfi); +#endif +#if AMREX_SPACEDIM == 2 + Array4 const dLogAreaY_arr = (dLogArea[1]).array(mfi); #endif const Box& xbx = amrex::surroundingNodes(bx, 0); @@ -342,8 +345,11 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co #if AMREX_SPACEDIM == 3 qzm_arr, qzp_arr, #endif -#if (AMREX_SPACEDIM < 3) - dLogArea_arr, +#if AMREX_SPACEDIM < 3 + dLogAreaX_arr, +#endif +#if AMREX_SPACEDIM == 2 + dLogAreaY_arr, #endif dt); @@ -361,7 +367,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co qzm_arr, qzp_arr, #endif #if AMREX_SPACEDIM < 3 - dLogArea_arr, + dLogAreaX_arr, +#endif +#if AMREX_SPACEDIM == 2 + dLogAreaY_arr, #endif dt); #else @@ -377,7 +386,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co qzm_arr, qzp_arr, #endif #if AMREX_SPACEDIM < 3 - dLogArea_arr, + dLogAreaX_arr, +#endif +#if AMREX_SPACEDIM == 2 + dLogAreaY_arr, #endif dt); #endif diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index cc7cd154cb..8b46a76c9e 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -218,7 +218,8 @@ /// @param qyp right interface state in y, q_{i,j-1/2,k,R} /// @param qzm left interface state in z, q_{i,j,k-1/2,L} /// @param qzp right interface state in z, q_{i,j,k-1/2,R} -/// @param dloga the geometric factor d(log area) +/// @param dlogaX the geometric factor d(log area) in x (r) dir +/// @param dlogaY the geometric factor d(log area) in y (theta) dir /// @param dt timestep /// void ctu_ppm_states(const amrex::Box& bx, const amrex::Box& vbx, @@ -238,7 +239,10 @@ amrex::Array4 const& qzp, #endif #if AMREX_SPACEDIM < 3 - amrex::Array4 const& dloga, + amrex::Array4 const& dlogaX, +#endif +#if AMREX_SPACEDIM == 2 + amrex::Array4 const& dlogaY, #endif const amrex::Real dt); @@ -259,7 +263,8 @@ /// @param qyp right interface state in y, q_{i,j-1/2,k,R} /// @param qzm left interface state in z, q_{i,j,k-1/2,L} /// @param qzp right interface state in z, q_{i,j,k-1/2,R} -/// @param dloga the geometric factor d(log area) +/// @param dlogaX the geometric factor d(log area) in x (r) dir +/// @param dlogaY the geometric factor d(log area) in y (theta) dir /// @param dt timestep /// void ctu_plm_states(const amrex::Box& bx, const amrex::Box& vbx, @@ -279,7 +284,10 @@ amrex::Array4 const& qzp, #endif #if AMREX_SPACEDIM < 3 - amrex::Array4 const& dloga, + amrex::Array4 const& dlogaX, +#endif +#if AMREX_SPACEDIM == 2 + amrex::Array4 const& dlogaY, #endif const amrex::Real dt); @@ -301,7 +309,8 @@ /// @param qyp right interface state in y, q_{i,j-1/2,k,R} /// @param qzm left interface state in z, q_{i,j,k-1/2,L} /// @param qzp right interface state in z, q_{i,j,k-1/2,R} -/// @param dloga the geometric factor d(log area) +/// @param dlogaX the geometric factor d(log area) in x (r) dir +/// @param dlogaY the geometric factor d(log area) in y (theta) dir /// @param dt timestep /// void ctu_ppm_rad_states(const amrex::Box& bx, const amrex::Box& vbx, @@ -321,7 +330,10 @@ amrex::Array4 const& qzp, #endif #if AMREX_SPACEDIM < 3 - amrex::Array4 const& dloga, + amrex::Array4 const& dlogaX, +#endif +#if AMREX_SPACEDIM == 2 + amrex::Array4 const& dlogaY, #endif const amrex::Real dt); #endif