diff --git a/Docs/source/mpi_plus_x.rst b/Docs/source/mpi_plus_x.rst index fe583bc8fb..7d88e7be77 100644 --- a/Docs/source/mpi_plus_x.rst +++ b/Docs/source/mpi_plus_x.rst @@ -106,6 +106,19 @@ To enable this, compile with:: to the ``make`` line or ``GNUmakefile``. +.. note:: + + CUDA 11.2 and later can do link time optimization. This can + increase performance by 10-30% (depending on the application), but + may greatly increase the compilation time. This is disabled by + default. To enable link time optimization, add: + + .. code:: + + CUDA_LTO=TRUE + + to the ``make`` line of ``GNUmakefile``. + AMD GPUs -------- @@ -123,4 +136,3 @@ Working at Supercomputing Centers Our best practices for running any of the AMReX Astrophysics codes at different supercomputing centers is produced in our workflow documentation: https://amrex-astro.github.io/workflow/ - diff --git a/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp b/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp index de815abb96..5c82caa1a8 100644 --- a/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp +++ b/Exec/hydro_tests/gamma_law_bubble/Problem_Derive.cpp @@ -52,11 +52,11 @@ void ca_derpi(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/, 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) / dat(i,j,k,URHO); + 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) / dat(i,j,k,URHO); + eos_state.aux[n] = dat(i,j,k,UFX+n) * rhoInv; } #endif 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..6a01ac74dc 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 @@ -33,6 +33,8 @@ castro.grav_source_type = 2 castro.use_pslope = 1 castro.pslope_cutoff_density = 1.e4 +castro.react_rho_min = 100.0 + gravity.gravity_type = ConstantGrav gravity.const_grav = -1.466e9 @@ -71,7 +73,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 diff --git a/Exec/science/xrb_layered/analysis/vol-zvel.py b/Exec/science/xrb_layered/analysis/vol-zvel.py index 013a772f75..f780d24391 100644 --- a/Exec/science/xrb_layered/analysis/vol-zvel.py +++ b/Exec/science/xrb_layered/analysis/vol-zvel.py @@ -62,7 +62,7 @@ def doit(plotfile): cam.position = [2.5*ds.domain_right_edge[0], 2.5*ds.domain_right_edge[1], - center[2]+0.25*ds.domain_right_edge[2]] + center[2]+0.5*ds.domain_right_edge[2]] # look toward the center -- we are dealing with an octant normal = (center - cam.position) @@ -71,6 +71,7 @@ def doit(plotfile): cam.switch_orientation(normal_vector=normal, north_vector=[0., 0., 1.]) cam.set_width(ds.domain_width) + cam.zoom(1.3) sc.camera = cam @@ -81,10 +82,11 @@ def doit(plotfile): sc.save(f"{plotfile}_zvel_axes.png", sigma_clip=3.0) sc.save_annotated(f"{plotfile}_zvel_annotated.png", - sigma_clip=4.0, label_fmt="%.2f", + sigma_clip=4.0, #label_fmt="%6.2g", + label_fontsize="16", text_annotate=[[(0.05, 0.05), - f"t = {ds.current_time.d:7.5f} s", - dict(horizontalalignment="left")]]) + f"t = {ds.current_time.d:6.4f} s", + dict(horizontalalignment="left", fontsize=16)]]) if __name__ == "__main__":