From 0f4a868944ecf430bb2a17167f7d77da2c022556 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 27 Jul 2024 21:09:13 -0400 Subject: [PATCH 1/6] clean up the artifical viscosity routines we just need to compute divu on the face once this is prep for limiting the artifical viscosity coefficient for species conservation --- Source/hydro/advection_util.cpp | 118 +++++++++++++++++--------------- 1 file changed, 62 insertions(+), 56 deletions(-) diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index ec7505a438..08a4112e41 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -302,49 +302,52 @@ Castro::apply_av(const Box& bx, Real diff_coeff = difmag; - amrex::ParallelFor(bx, NUM_STATE, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - if (n == UTEMP) { - return; - } + Real div1; + if (idir == 0) { + div1 = 0.25_rt * (div(i,j,k) + div(i,j+dg1,k) + + div(i,j,k+dg2) + div(i,j+dg1,k+dg2)); + } else if (idir == 1) { + div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + + div(i,j,k+dg2) + div(i+1,j,k+dg2)); + } else { + div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + + div(i,j+dg1,k) + div(i+1,j+dg1,k)); + } + + div1 = diff_coeff * std::min(0.0_rt, div1); + + for (int n = 0; n < NUM_STATE, ++n) { + + if (n == UTEMP) { + continue; + } #ifdef SHOCK_VAR - if (n == USHK) { - return; - } + if (n == USHK) { + continue; + } #endif - #ifdef NSE_NET - if (n == UMUP || n == UMUN) { - return; - } + if (n == UMUP || n == UMUN) { + continue; + } #endif - Real div1; - if (idir == 0) { - - div1 = 0.25_rt * (div(i,j,k) + div(i,j+dg1,k) + - div(i,j,k+dg2) + div(i,j+dg1,k+dg2)); - div1 = diff_coeff * amrex::min(0.0_rt, div1); - div1 = div1 * (uin(i,j,k,n) - uin(i-1,j,k,n)); - } else if (idir == 1) { + Real div_var{}; - div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + - div(i,j,k+dg2) + div(i+1,j,k+dg2)); - div1 = diff_coeff * amrex::min(0.0_rt, div1); - div1 = div1 * (uin(i,j,k,n) - uin(i,j-dg1,k,n)); - - } else { + if (idir == 0) { + div_var = div1 * (uin(i,j,k,n) - uin(i-1,j,k,n)); + } else if (idir == 1) { + div_var = div1 * (uin(i,j,k,n) - uin(i,j-dg1,k,n)); + } else { + div_var = div1 * (uin(i,j,k,n) - uin(i,j,k-dg2,n)); + } - div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + - div(i,j+dg1,k) + div(i+1,j+dg1,k)); - div1 = diff_coeff * amrex::min(0.0_rt, div1); - div1 = div1 * (uin(i,j,k,n) - uin(i,j,k-dg2,n)); - - } - - flux(i,j,k,n) += dx[idir] * div1; + flux(i,j,k,n) += dx[idir] * div_var; + } }); } @@ -361,35 +364,38 @@ Castro::apply_av_rad(const Box& bx, Real diff_coeff = difmag; - amrex::ParallelFor(bx, Radiation::nGroups, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real div1; - if (idir == 0) { - - div1 = 0.25_rt * (div(i,j,k) + div(i,j+dg1,k) + - div(i,j,k+dg2) + div(i,j+dg1,k+dg2)); - div1 = diff_coeff * amrex::min(0.0_rt, div1); - div1 = div1 * (Erin(i,j,k,n) - Erin(i-1,j,k,n)); + Real div1; + if (idir == 0) { + div1 = 0.25_rt * (div(i,j,k) + div(i,j+dg1,k) + + div(i,j,k+dg2) + div(i,j+dg1,k+dg2)); + } else if (idir == 1) { + div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + + div(i,j,k+dg2) + div(i+1,j,k+dg2)); + } else { + div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + + div(i,j+dg1,k) + div(i+1,j+dg1,k)); + } - } else if (idir == 1) { + div1 = diff_coeff * std::min(0.0_rt, div1); - div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + - div(i,j,k+dg2) + div(i+1,j,k+dg2)); - div1 = diff_coeff * amrex::min(0.0_rt, div1); - div1 = div1 * (Erin(i,j,k,n) - Erin(i,j-dg1,k,n)); + for (int n = 0; n < Radiation::nGroups; ++n) { - } else { + Real div_var{}; - div1 = 0.25_rt * (div(i,j,k) + div(i+1,j,k) + - div(i,j+dg1,k) + div(i+1,j+dg1,k)); - div1 = diff_coeff * amrex::min(0.0_rt, div1); - div1 = div1 * (Erin(i,j,k,n) - Erin(i,j,k-dg2,n)); + if (idir == 0) { + div_var = div1 * (Erin(i,j,k,n) - Erin(i-1,j,k,n)); + } else if (idir == 1) { + div_var = div1 * (Erin(i,j,k,n) - Erin(i,j-dg1,k,n)); + } else { + div_var = div1 * (Erin(i,j,k,n) - Erin(i,j,k-dg2,n)); + } - } - - radflux(i,j,k,n) += dx[idir] * div1; + radflux(i,j,k,n) += dx[idir] * div_var; + } }); } #endif From 660f4e97f27784cc3d3cde6223f4de1994fad8b0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 27 Jul 2024 21:21:56 -0400 Subject: [PATCH 2/6] fix compilation --- Source/hydro/advection_util.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 08a4112e41..52eb3ba749 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -320,7 +320,7 @@ Castro::apply_av(const Box& bx, div1 = diff_coeff * std::min(0.0_rt, div1); - for (int n = 0; n < NUM_STATE, ++n) { + for (int n = 0; n < NUM_STATE; ++n) { if (n == UTEMP) { continue; From 9253308885f847a8d8fc05120de18c55d96e9143 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 28 Jul 2024 14:58:08 -0400 Subject: [PATCH 3/6] limit artificial viscosity if it breaks species conservation --- Source/driver/_cpp_parameters | 5 ++++ Source/hydro/advection_util.cpp | 41 +++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index f5cdc29bd3..b4a0d87a10 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -39,6 +39,11 @@ allow_non_unit_aspect_zones int 0 # the coefficient of the artificial viscosity difmag Real 0.1 +# when applying the artifical viscosity, do we limit the magnitude of +# the diffusion coefficient to ensure that the species fluxes all keep +# the same sign? +limit_avisc_on_species bool 0 + # the small density cutoff. Densities below this value will be reset small_dens Real -1.e200 diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 52eb3ba749..a6cf5e682b 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -320,6 +320,43 @@ Castro::apply_av(const Box& bx, div1 = diff_coeff * std::min(0.0_rt, div1); + if (div1 == 0.0_rt) { + return; + } + + // we want to prevent the artificial viscosity from breaking + // species conservation (sum X_k = 1). This can happen if the + // artificial viscosity flips the sign of some of the fluxes. + // Here we scale div1 to ensure that this is not an issue. + + if (limit_avisc_on_species) { + + amrex::Real limit_div1{div1}; + + for (int n = 0; n < NumSpec; ++n) { + + Real du{}; + + if (idir == 0) { + du = uin(i,j,k,UFS+n) - uin(i-1,j,k,UFS+n); + } else if (idir == 1) { + du = uin(i,j,k,UFS+n) - uin(i,j-dg1,k,UFS+n); + } else { + du = uin(i,j,k,UFS+n) - uin(i,j,k-dg2,UFS+n); + } + + Real fnew = flux(i,j,k,UFS+n) + dx[idir] * div1 * du; + if (fnew * flux(i,j,k,UFS+n) < 0) { + limit_div1 = std::min(limit_div1, + std::abs(flux(i,j,k,UFS+n) / (dx[idir] * du))); + } + } + + div1 = limit_div1; + + } + + for (int n = 0; n < NUM_STATE; ++n) { if (n == UTEMP) { @@ -382,6 +419,10 @@ Castro::apply_av_rad(const Box& bx, div1 = diff_coeff * std::min(0.0_rt, div1); + if (div1 == 0.0_rt) { + return; + } + for (int n = 0; n < Radiation::nGroups; ++n) { Real div_var{}; From 76ed84696cd96bf2fb181f4ac4dc2e2b7406810f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 28 Jul 2024 15:02:18 -0400 Subject: [PATCH 4/6] codespell --- Source/driver/_cpp_parameters | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index b4a0d87a10..edd7245c03 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -39,7 +39,7 @@ allow_non_unit_aspect_zones int 0 # the coefficient of the artificial viscosity difmag Real 0.1 -# when applying the artifical viscosity, do we limit the magnitude of +# when applying the artificial viscosity, do we limit the magnitude of # the diffusion coefficient to ensure that the species fluxes all keep # the same sign? limit_avisc_on_species bool 0 From 8949302c6adcaaa23f56715098a2782cf7f0981c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 28 Jul 2024 16:06:29 -0400 Subject: [PATCH 5/6] turn it on --- Source/driver/_cpp_parameters | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index edd7245c03..5279f85f6a 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -42,7 +42,7 @@ difmag Real 0.1 # when applying the artificial viscosity, do we limit the magnitude of # the diffusion coefficient to ensure that the species fluxes all keep # the same sign? -limit_avisc_on_species bool 0 +limit_avisc_on_species bool 1 # the small density cutoff. Densities below this value will be reset small_dens Real -1.e200 From 4410afdad0bb2b0b25f31049e89dbbe53e24b605 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 28 Jul 2024 16:21:22 -0400 Subject: [PATCH 6/6] fix capture --- Source/hydro/advection_util.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 52eb3ba749..095192931b 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -362,7 +362,9 @@ Castro::apply_av_rad(const Box& bx, const auto dx = geom.CellSizeArray(); - Real diff_coeff = difmag; + amrex::Real diff_coeff = difmag; + + int ngroups = Radiation::nGroups; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -382,7 +384,7 @@ Castro::apply_av_rad(const Box& bx, div1 = diff_coeff * std::min(0.0_rt, div1); - for (int n = 0; n < Radiation::nGroups; ++n) { + for (int n = 0; n < ngroups; ++n) { Real div_var{};