From 50b362169db508859a596d4a099f5c9f53c6ff36 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 21 Jun 2024 14:01:45 -0400 Subject: [PATCH] add a derived variable for gradP_over_P (#2877) this will help understand the shock flag --- Exec/science/subch_planar/Problem_Derive.H | 10 + Exec/science/subch_planar/Problem_Derive.cpp | 413 ++++++++++++++++++ Exec/science/subch_planar/Problem_Derives.H | 8 + Exec/science/subch_planar/inputs_2d | 4 +- .../science/subch_planar/problem_initialize.H | 5 +- .../problem_initialize_state_data.H | 5 +- 6 files changed, 436 insertions(+), 9 deletions(-) create mode 100644 Exec/science/subch_planar/Problem_Derive.H create mode 100644 Exec/science/subch_planar/Problem_Derive.cpp create mode 100644 Exec/science/subch_planar/Problem_Derives.H diff --git a/Exec/science/subch_planar/Problem_Derive.H b/Exec/science/subch_planar/Problem_Derive.H new file mode 100644 index 0000000000..2110a535a1 --- /dev/null +++ b/Exec/science/subch_planar/Problem_Derive.H @@ -0,0 +1,10 @@ + +void ca_dergradpoverp + (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, + const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); + +void ca_dergradpoverp1 + (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, + const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); diff --git a/Exec/science/subch_planar/Problem_Derive.cpp b/Exec/science/subch_planar/Problem_Derive.cpp new file mode 100644 index 0000000000..6a67c0532f --- /dev/null +++ b/Exec/science/subch_planar/Problem_Derive.cpp @@ -0,0 +1,413 @@ +#include + +#include +#include + +using namespace amrex; + +using RealVector = amrex::Gpu::ManagedVector; + +void ca_dergradpoverp(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, + const FArrayBox& datfab, const Geometry& geomdata, + Real /*time*/, const int* /*bcrec*/, int /*level*/) +{ + + // Compute grad p . U / (p |U|) -- this is what we do in the shock + // detection algorithm + + const auto dx = geomdata.CellSizeArray(); + const int coord_type = geomdata.Coord(); + + auto const dat = datfab.array(); + auto const der = derfab.array(); + +#if AMREX_SPACEDIM == 3 + amrex::Error("3D not supported"); +#endif +#if AMREX_SPACEDIM == 1 + amrex::Error("1D not supported"); +#endif + + Real dxinv = 1.0_rt / dx[0]; + Real dyinv = 1.0_rt / dx[1]; + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + Real div_u = 0.0_rt; + + // get the pressure + eos_t eos_state; + + Real rhoInv = 1.0 / dat(i,j,k,URHO); + + eos_state.rho = dat(i,j,k,URHO); + eos_state.T = dat(i,j,k,UTEMP); + eos_state.e = dat(i,j,k,UEINT) * rhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + Real p_zone = eos_state.p; + + Real up = dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO); + Real um = dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO); + Real u0 = dat(i,j,k,UMX) / dat(i,j,k,URHO); + + Real vp = dat(i,j+1,k,UMY) / dat(i,j+1,k,URHO); + Real vm = dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO); + Real v0 = dat(i,j,k,UMY) / dat(i,j,k,URHO); + + // construct div{U} + if (coord_type == 0) { + + // Cartesian + div_u += 0.5_rt * (up - um) * dxinv; + div_u += 0.5_rt * (vp - vm) * dyinv; + + } else if (coord_type == 1) { + + // r-z + Real rc = (i + 0.5_rt) * dx[0]; + Real rm = (i - 1 + 0.5_rt) * dx[0]; + Real rp = (i + 1 + 0.5_rt) * dx[0]; + + div_u += 0.5_rt * (rp * up - rm * um) / (rc * dx[0]) + + 0.5_rt * (vp - vm) * dyinv; + +#ifndef AMREX_USE_GPU + } else { + amrex::Error("ERROR: invalid coord_type in shock"); +#endif + } + + + // we need to compute p in the full stencil + + Real p_ip1{}; + { + Real lrhoInv = 1.0 / dat(i+1,j,k,URHO); + + eos_state.rho = dat(i+1,j,k,URHO); + eos_state.T = dat(i+1,j,k,UTEMP); + eos_state.e = dat(i+1,j,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i+1,j,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i+1,j,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_ip1 = eos_state.p; + } + + Real p_im1{}; + { + Real lrhoInv = 1.0 / dat(i-1,j,k,URHO); + + eos_state.rho = dat(i-1,j,k,URHO); + eos_state.T = dat(i-1,j,k,UTEMP); + eos_state.e = dat(i-1,j,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i-1,j,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i-1,j,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_im1 = eos_state.p; + } + + Real dP_x = 0.5_rt * (p_ip1 - p_im1); + + Real p_jp1{}; + { + Real lrhoInv = 1.0 / dat(i,j+1,k,URHO); + + eos_state.rho = dat(i,j+1,k,URHO); + eos_state.T = dat(i,j+1,k,UTEMP); + eos_state.e = dat(i,j+1,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i,j+1,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i,j+1,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_jp1 = eos_state.p; + } + + Real p_jm1{}; + { + Real lrhoInv = 1.0 / dat(i,j-1,k,URHO); + + eos_state.rho = dat(i,j-1,k,URHO); + eos_state.T = dat(i,j-1,k,UTEMP); + eos_state.e = dat(i,j-1,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i,j-1,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i,j-1,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_jm1 = eos_state.p; + } + + Real dP_y = 0.5_rt * (p_jp1 - p_jm1); + + //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / dat(i,j,k,QPRES); + Real vel = std::sqrt(u0 * u0 + v0 * v0); + + Real gradPdx_over_P{0.0_rt}; + if (vel != 0.0) { + gradPdx_over_P = std::abs(dP_x * u0 + dP_y * v0) / vel; + } + gradPdx_over_P /= p_zone; + + der(i,j,k,0) = gradPdx_over_P; + + }); + +} + +void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, + const FArrayBox& datfab, const Geometry& geomdata, + Real /*time*/, const int* /*bcrec*/, int /*level*/) +{ + + // compute grad p projected onto div{U} as an alternative to the shock detection + + const auto dx = geomdata.CellSizeArray(); + const int coord_type = geomdata.Coord(); + + auto const dat = datfab.array(); + auto const der = derfab.array(); + +#if AMREX_SPACEDIM == 3 + amrex::Error("3D not supported"); +#endif +#if AMREX_SPACEDIM == 1 + amrex::Error("1D not supported"); +#endif + + Real dxinv = 1.0_rt / dx[0]; + Real dyinv = 1.0_rt / dx[1]; + + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + Real div_u = 0.0_rt; + + // get the pressure + eos_t eos_state; + + Real rhoInv = 1.0 / dat(i,j,k,URHO); + + eos_state.rho = dat(i,j,k,URHO); + eos_state.T = dat(i,j,k,UTEMP); + eos_state.e = dat(i,j,k,UEINT) * rhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i,j,k,UFS+n) * rhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + Real p_zone = eos_state.p; + + Real up = dat(i+1,j,k,UMX) / dat(i+1,j,k,URHO); + Real um = dat(i-1,j,k,UMX) / dat(i-1,j,k,URHO); + Real u0 = dat(i,j,k,UMX) / dat(i,j,k,URHO); + + Real vp = dat(i,j+1,k,UMY) / dat(i,j+1,k,URHO); + Real vm = dat(i,j-1,k,UMY) / dat(i,j-1,k,URHO); + Real v0 = dat(i,j,k,UMY) / dat(i,j,k,URHO); + + // construct div{U} + if (coord_type == 0) { + + // Cartesian + div_u += 0.5_rt * (up - um) * dxinv; + div_u += 0.5_rt * (vp - vm) * dyinv; + + } else if (coord_type == 1) { + + // r-z + Real rc = (i + 0.5_rt) * dx[0]; + Real rm = (i - 1 + 0.5_rt) * dx[0]; + Real rp = (i + 1 + 0.5_rt) * dx[0]; + + div_u += 0.5_rt * (rp * up - rm * um) / (rc * dx[0]) + + 0.5_rt * (vp - vm) * dyinv; + +#ifndef AMREX_USE_GPU + } else { + amrex::Error("ERROR: invalid coord_type in shock"); +#endif + } + + // we need to compute p in the full stencil + + Real p_ip1{}; + { + Real lrhoInv = 1.0 / dat(i+1,j,k,URHO); + + eos_state.rho = dat(i+1,j,k,URHO); + eos_state.T = dat(i+1,j,k,UTEMP); + eos_state.e = dat(i+1,j,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i+1,j,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i+1,j,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_ip1 = eos_state.p; + } + + Real p_im1{}; + { + Real lrhoInv = 1.0 / dat(i-1,j,k,URHO); + + eos_state.rho = dat(i-1,j,k,URHO); + eos_state.T = dat(i-1,j,k,UTEMP); + eos_state.e = dat(i-1,j,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i-1,j,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i-1,j,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_im1 = eos_state.p; + } + + Real dP_x = 0.5_rt * (p_ip1 - p_im1); + + Real p_jp1{}; + { + Real lrhoInv = 1.0 / dat(i,j+1,k,URHO); + + eos_state.rho = dat(i,j+1,k,URHO); + eos_state.T = dat(i,j+1,k,UTEMP); + eos_state.e = dat(i,j+1,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i,j+1,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i,j+1,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_jp1 = eos_state.p; + } + + Real p_jm1{}; + { + Real lrhoInv = 1.0 / dat(i,j-1,k,URHO); + + eos_state.rho = dat(i,j-1,k,URHO); + eos_state.T = dat(i,j-1,k,UTEMP); + eos_state.e = dat(i,j-1,k,UEINT) * lrhoInv; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = dat(i,j-1,k,UFS+n) * lrhoInv; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = dat(i,j-1,k,UFX+n) * lrhoInv; + } +#endif + + if (eos_state.e <= 0.0_rt) { + eos(eos_input_rt, eos_state); + } else { + eos(eos_input_re, eos_state); + } + p_jm1 = eos_state.p; + } + + Real dP_y = 0.5_rt * (p_jp1 - p_jm1); + + //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / dat(i,j,k,QPRES); + Real du_x = std::min(up - um, 0.0); + Real dv_y = std::min(vp - vm, 0.0); + + Real divu_mag = std::sqrt(du_x * du_x + dv_y * dv_y + 1.e-30); + + Real gradPdx_over_P = std::abs(dP_x * du_x + dP_y * dv_y) / divu_mag; + gradPdx_over_P /= p_zone; + + der(i,j,k,0) = gradPdx_over_P; + + }); + +} diff --git a/Exec/science/subch_planar/Problem_Derives.H b/Exec/science/subch_planar/Problem_Derives.H new file mode 100644 index 0000000000..ae147a39c1 --- /dev/null +++ b/Exec/science/subch_planar/Problem_Derives.H @@ -0,0 +1,8 @@ + // + // gradpoverp + // + derive_lst.add("gradp_over_p", IndexType::TheCellType(), 1, ca_dergradpoverp, grow_box_by_one); + derive_lst.addComponent("gradp_over_p", desc_lst, State_Type, URHO, NUM_STATE); + + derive_lst.add("gradp_over_p1", IndexType::TheCellType(), 1, ca_dergradpoverp1, grow_box_by_one); + derive_lst.addComponent("gradp_over_p1", desc_lst, State_Type, URHO, NUM_STATE); diff --git a/Exec/science/subch_planar/inputs_2d b/Exec/science/subch_planar/inputs_2d index 56d6d9854f..d1eb341121 100644 --- a/Exec/science/subch_planar/inputs_2d +++ b/Exec/science/subch_planar/inputs_2d @@ -7,7 +7,7 @@ geometry.is_periodic = 0 0 geometry.coord_sys = 1 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 0.0 0.0 geometry.prob_hi = 4.e8 4.e8 -amr.n_cell = 200 200 +amr.n_cell = 400 400 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< # 0 = Interior 3 = Symmetry @@ -71,7 +71,7 @@ amr.plot_per = 0.002 #simulation time (s) between plotfiles amr.derive_plot_vars = ALL # PROBLEM PARAMETERS -problem.model_name = "toy_subch.hse.tanh.delta_50.00km.dx_20.00km" +problem.model_name = "toy_subch.hse.tanh.delta_50.00km.dx_10.00km" problem.pert_temp_factor = 20.0 problem.pert_rad_factor = 0.5 diff --git a/Exec/science/subch_planar/problem_initialize.H b/Exec/science/subch_planar/problem_initialize.H index 01bced54c2..7a8dc5e7b3 100644 --- a/Exec/science/subch_planar/problem_initialize.H +++ b/Exec/science/subch_planar/problem_initialize.H @@ -12,9 +12,6 @@ void problem_initialize () const Geometry& dgeom = DefaultGeometry(); - const Real* problo = dgeom.ProbLo(); - const Real* probhi = dgeom.ProbHi(); - problem::ihe4 = network_spec_index("helium-4"); if (problem::ihe4 < 0) { @@ -70,4 +67,4 @@ void problem_initialize () amrex::Print() << "base of He layer found at " << problem::R_He_base << std::endl; } -#endif \ No newline at end of file +#endif diff --git a/Exec/science/subch_planar/problem_initialize_state_data.H b/Exec/science/subch_planar/problem_initialize_state_data.H index 0f29bd07e1..f002f1339c 100644 --- a/Exec/science/subch_planar/problem_initialize_state_data.H +++ b/Exec/science/subch_planar/problem_initialize_state_data.H @@ -21,9 +21,8 @@ void problem_initialize_state_data (int i, int j, int k, y = problo[1] + dx[1] * (static_cast(j) + 0.5_rt); #endif - Real z = 0.0; #if AMREX_SPACEDIM == 3 - z = problo[2] + dx[2] * (static_cast(k) + 0.5_rt); + Real z = problo[2] + dx[2] * (static_cast(k) + 0.5_rt); #endif #if AMREX_SPACEDIM == 2 @@ -123,4 +122,4 @@ void problem_initialize_state_data (int i, int j, int k, state(i,j,k,UEINT) = burn_state.e * state(i,j,k,URHO); state(i,j,k,UEDEN) = burn_state.e * state(i,j,k,URHO); } -#endif \ No newline at end of file +#endif