diff --git a/CHANGES.md b/CHANGES.md index 7b7e78a0da..665eff135e 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,20 @@ +# 24.09 + + * Code clean-ups / clang-tidy (#2942, #2949) + + * update the `hse_convergence` readme to reflect current convergence + (#2946) + + * update the `bubble_convergence` plotting script (#2947) + + * new Frontier scaling numbers (#2948) + + * more GPU error printing (@3944) + + * science problem updates: `flame_wave` (#2943) + + * documentation updates (#2939) + # 24.08 * lazy QueueReduction has been enabled for the timing diagnostics diff --git a/Exec/gravity_tests/hse_convergence/README.md b/Exec/gravity_tests/hse_convergence/README.md index fa639f4708..fab88d7bd7 100644 --- a/Exec/gravity_tests/hse_convergence/README.md +++ b/Exec/gravity_tests/hse_convergence/README.md @@ -7,29 +7,34 @@ in the plotfiles. To run this problem, use one of the convergence scripts: - * ``convergence_plm.sh`` : + * `convergence_plm.sh` : - this runs CTU + PLM using the default HSE BCs and default - use_pslope, then with reflect BCs, then without use_pslope, and - finally runs with reflect instead of HSE BCs. + this runs CTU + PLM using: + 1. the default HSE BCs and `use_pslope` + 2. the HSE BCs with reflection and `use_pslope` + 3. reflect BCs instead of HSE BCs without `use_pslope` + 4. reflect BCs with `use_pslope` - These tests show that the best results come from HSE BCs + reflect vel + These tests show that the best results (by far) come from + `use_pslope=1` and reflecting BCs * convergence_ppm.sh : this runs CTU + PPM in a similar set of configurations as PLM above - (with one additional one: grav_source_type = 4) + 1. the default HSE BCs + 2. HSE BCs with reflection + 3. reflecting BCs + 4. reflecting BCs with `use_pslope` - These tests show that the best results come from HSE BCs + reflect vel + These tests show that the best results (by far) come from + reflecting BCs with `use_pslope=1`, just like the PLM case. * convergence_sdc.sh : - this uses the TRUE_SDC integration, first with SDC-2 + PLM and reflecting BCs, - the SDC-2 + PPM and reflecting BCs, then the same but HSE BCs, and finally - SDC-4 + reflect + this uses the TRUE_SDC integration, first with SDC-2 + PLM and + reflecting BCs, the SDC-2 + PPM and reflecting BCs, then the same + but HSE BCs, and finally SDC-4 + reflect These tests show that the PLM + reflect (which uses the well-balanced use_pslope) and the SDC-4 + reflect give the lowest - errors and expected (or better) convergence: - - + errors and expected (or better) convergence. diff --git a/Exec/gravity_tests/hse_convergence/convergence_plm.sh b/Exec/gravity_tests/hse_convergence/convergence_plm.sh index 285cfed67d..3d12427ad5 100755 --- a/Exec/gravity_tests/hse_convergence/convergence_plm.sh +++ b/Exec/gravity_tests/hse_convergence/convergence_plm.sh @@ -58,43 +58,15 @@ pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} -## plm + hse reflect + no pslope - -ofile=plm-hsereflect-nopslope.converge.out - -RUNPARAMS=" -castro.ppm_type=0 -castro.use_pslope=0 -castro.hse_interp_temp=1 -castro.hse_reflect_vels=1 -""" - -${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out -pfile=`ls -t | grep -i hse_64_plt | head -1` -fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel > ${ofile} - -${EXEC} inputs.ppm.128 ${RUNPARAMS} >& 128.out -pfile=`ls -t | grep -i hse_128_plt | head -1` -fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} - -${EXEC} inputs.ppm.256 ${RUNPARAMS} >& 256.out -pfile=`ls -t | grep -i hse_256_plt | head -1` -fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} - -${EXEC} inputs.ppm.512 ${RUNPARAMS} >& 512.out -pfile=`ls -t | grep -i hse_512_plt | head -1` -fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} - - -## plm + reflect +## plm + reflect + nopslope -ofile=plm-reflect.converge.out +ofile=plm-reflect-nopslope.converge.out RUNPARAMS=" castro.ppm_type=0 -castro.use_pslope=1 castro.lo_bc=3 castro.hi_bc=3 +castro.use_pslope=0 """ ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out @@ -114,16 +86,15 @@ pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} +## plm + reflect + pslope -## plm + reflect + nopslope - -ofile=plm-reflect-nopslope.converge.out +ofile=plm-reflect-pslope.converge.out RUNPARAMS=" castro.ppm_type=0 castro.lo_bc=3 castro.hi_bc=3 -castro.use_pslope=0 +castro.use_pslope=1 """ ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out diff --git a/Exec/gravity_tests/hse_convergence/convergence_ppm.sh b/Exec/gravity_tests/hse_convergence/convergence_ppm.sh index ff6e2c7620..27b4bb86ef 100755 --- a/Exec/gravity_tests/hse_convergence/convergence_ppm.sh +++ b/Exec/gravity_tests/hse_convergence/convergence_ppm.sh @@ -50,12 +50,13 @@ pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} -## ppm + grav_source_type = 4 +## ppm + reflect -ofile=ppm-grav4.converge.out +ofile=ppm-reflect.converge.out RUNPARAMS=" -castro.grav_source_type=4 +castro.lo_bc=3 +castro.hi_bc=3 """ ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out @@ -75,13 +76,14 @@ pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} -## ppm + reflect +## ppm + reflect + pslope -ofile=ppm-reflect.converge.out +ofile=ppm-reflect-pslope.converge.out RUNPARAMS=" castro.lo_bc=3 castro.hi_bc=3 +castro.use_pslope=1 """ ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out diff --git a/Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H b/Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H index adc5a4c050..f87fb240af 100644 --- a/Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H +++ b/Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H @@ -73,7 +73,7 @@ void problem_initialize_state_data (int i, int j, int k, } u_tot += u_phi; - reint += p/(gamma_const - 1.0_rt); + reint += p/(eos_rp::eos_gamma - 1.0_rt); } } } diff --git a/Exec/radiation_tests/Rad2Tshock/problem_initialize.H b/Exec/radiation_tests/Rad2Tshock/problem_initialize.H index 1405a5d6c2..4cffa04c3a 100644 --- a/Exec/radiation_tests/Rad2Tshock/problem_initialize.H +++ b/Exec/radiation_tests/Rad2Tshock/problem_initialize.H @@ -14,8 +14,8 @@ void problem_initialize () eos_state.rho = problem::rho0; eos_state.T = problem::T0; - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = 0.0_rt; + for (auto & X : eos_state.xn) { + X = 0.0_rt; } eos_state.xn[0] = 1.0_rt; diff --git a/Exec/radiation_tests/Rad2Tshock/problem_initialize_rad_data.H b/Exec/radiation_tests/Rad2Tshock/problem_initialize_rad_data.H index 860d2dcd60..87b4721254 100644 --- a/Exec/radiation_tests/Rad2Tshock/problem_initialize_rad_data.H +++ b/Exec/radiation_tests/Rad2Tshock/problem_initialize_rad_data.H @@ -12,6 +12,9 @@ void problem_initialize_rad_data (int i, int j, int k, const GeometryData& geomdata) { + amrex::ignore_unused(nugroup); + amrex::ignore_unused(dnugroup); + const Real* dx = geomdata.CellSize(); const Real* problo = geomdata.ProbLo(); diff --git a/Exec/radiation_tests/Rad2Tshock/problem_initialize_state_data.H b/Exec/radiation_tests/Rad2Tshock/problem_initialize_state_data.H index 96fe4057e0..c5e67920a5 100644 --- a/Exec/radiation_tests/Rad2Tshock/problem_initialize_state_data.H +++ b/Exec/radiation_tests/Rad2Tshock/problem_initialize_state_data.H @@ -16,7 +16,7 @@ void problem_initialize_state_data (int i, int j, int k, // Provides the simulation to be run in the x,y,or z direction // where length direction is the length side in a square prism - Real length_cell; + Real length_cell{}; if (problem::idir == 1) { length_cell = problo[0] + dx[0] * (static_cast(i) + 0.5_rt); } else if (problem::idir == 2) { diff --git a/Exec/reacting_tests/bubble_convergence/analysis/slice_multi.py b/Exec/reacting_tests/bubble_convergence/analysis/slice_multi.py index d05bd7fba9..2f8c1e16b8 100755 --- a/Exec/reacting_tests/bubble_convergence/analysis/slice_multi.py +++ b/Exec/reacting_tests/bubble_convergence/analysis/slice_multi.py @@ -1,5 +1,8 @@ #!/usr/bin/env python3 +import matplotlib +matplotlib.use('agg') + import os import sys import yt @@ -26,35 +29,37 @@ fig = plt.figure() fig.set_size_inches(12.0, 9.0) -grid = ImageGrid(fig, 111, nrows_ncols=(2, 2), axes_pad=0.75, cbar_pad="2%", +grid = ImageGrid(fig, 111, nrows_ncols=(2, 2), + axes_pad=0.75, cbar_pad="2%", label_mode="L", cbar_mode="each") -fields = ["Temp", "magvel", "X(C12)", "rho_enuc"] +fields = ["Temp", "magvel", "X(C12)", "enuc"] for i, f in enumerate(fields): - sp = yt.SlicePlot(ds, "z", f, center=[xctr, yctr, 0.0], width=[L_x, L_y, 0.0], fontsize="12") + sp = yt.SlicePlot(ds, "z", f, center=[xctr, yctr, 0.0*cm], + width=[L_x, L_y, 0.0*cm], fontsize="12") sp.set_buff_size((2000,2000)) if f == "X(C12)": sp.set_log(f, True) - sp.set_cmap(f, "plasma") - sp.set_zlim(f, 1.e-8, 2.e-4) + sp.set_cmap(f, "magma") + sp.set_zlim(f, 1.e-8, 1.e-4) elif f == "magvel": sp.set_log(f, False) #sp.set_zlim(f, 1.e-3, 2.5e-2) - sp.set_cmap(f, "magma") + sp.set_cmap(f, "cividis") elif f == "Temp": - sp.set_log(f, False) - #sp.set_zlim(f, 1.e-3, 2.5e-2) - - elif f == "rho_enuc": sp.set_log(f, True) + sp.set_zlim(f, 5.e7, 2.e8) + + elif f == "enuc": + sp.set_log(f, True, linthresh=1.e11) + sp.set_zlim(f, 1.e11, 1.e14) sp.set_cmap(f, "plasma") - #sp.set_zlim(f, 1.e-3, 2.5e-2) sp.set_axes_unit("cm") @@ -71,5 +76,4 @@ fig.set_size_inches(8.0, 8.0) plt.tight_layout() -plt.savefig("{}_slice.pdf".format(os.path.basename(plotfile))) - +plt.savefig("{}_slice.png".format(os.path.basename(plotfile))) diff --git a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt index 7d3acf86ec..ff734d8ae0 100644 --- a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt +++ b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt @@ -94,7 +94,7 @@ castro.sdc_quadrature = 0 castro.sdc_extra = 0 castro.sdc_solver = 1 - castro.use_axisymmetric_geom_source = 1 + castro.use_geom_source = 1 castro.add_sdc_react_source_to_advection = 1 castro.hydro_memory_footprint_ratio = -1 castro.fixed_dt = -1 diff --git a/Exec/science/flame_wave/scaling/frontier/frontier-scaling-2024-08-21.txt b/Exec/science/flame_wave/scaling/frontier/frontier-scaling-2024-08-21.txt new file mode 100644 index 0000000000..660b7af7b9 --- /dev/null +++ b/Exec/science/flame_wave/scaling/frontier/frontier-scaling-2024-08-21.txt @@ -0,0 +1,39 @@ +# new scaling numbers for the 3D XRB +# using the same inputs.He.25cm.static.1000Hz as previously +# modules: +# +# module load PrgEnv-gnu +# module load cray-mpich/8.1.28 +# module load craype-accel-amd-gfx90a +# module load amd-mixed/6.1.3 +# module unload darshan-runtime +# +# build info: +# +# EOS: /ccs/home/zingale/Microphysics/EOS/helmholtz +# NETWORK: /ccs/home/zingale/Microphysics/networks/iso7 +# CONDUCTIVITY: /ccs/home/zingale/Microphysics/conductivity/stellar +# INTEGRATOR: VODE +# SCREENING: screen5 +# +# Castro git describe: 24.08-3-g15327db6b +# AMReX git describe: 24.08-25-g6dcaa1223 +# Microphysics git describe: 24.08-2-g8ce3375a + +# nodes rocm mag_grid_size avg time / std dev +# step + 48 6.1.3 128 59.0711 0.2525 + 64 6.1.3 128 42.6938 0.285659 + 128 6.1.3 128 24.5353 1.36496 + 256 6.1.3 128 13.3647 0.108731 + 512 6.1.3 128 7.88166 0.0856889 +1024 6.1.3 128 5.54221 0.0979851 +2048 6.1.3 128 4.55679 0.0528629 + + +# note that the 2048 run uses a blocking factor of 16) + +# in contrast to the previous run, we've disabled all inlining with +# ROCm to get around some compiler bugs, so that might explain some +# slight slowdown here. + diff --git a/Exec/science/subch_planar/Problem_Derive.cpp b/Exec/science/subch_planar/Problem_Derive.cpp index 3275166238..9cb25c92e7 100644 --- a/Exec/science/subch_planar/Problem_Derive.cpp +++ b/Exec/science/subch_planar/Problem_Derive.cpp @@ -23,9 +23,8 @@ void ca_dergradpoverp(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco #if AMREX_SPACEDIM == 3 amrex::Error("3D not supported"); -#endif -#if AMREX_SPACEDIM == 1 - amrex::Error("1D not supported"); +#elif AMREX_SPACEDIM == 1 + return; // Skip for 1D #endif Real dxinv = 1.0_rt / dx[0]; @@ -226,9 +225,8 @@ void ca_dergradpoverp1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nc #if AMREX_SPACEDIM == 3 amrex::Error("3D not supported"); -#endif -#if AMREX_SPACEDIM == 1 - amrex::Error("1D not supported"); +#elif AMREX_SPACEDIM == 1 + return; // Skip for 1D #endif Real dxinv = 1.0_rt / dx[0]; @@ -434,9 +432,8 @@ void ca_dergradpx(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/ #if AMREX_SPACEDIM == 3 amrex::Error("3D not supported"); -#endif -#if AMREX_SPACEDIM == 1 - amrex::Error("1D not supported"); +#elif AMREX_SPACEDIM == 1 + return; // Skip for 1D #endif Real dxinv = 1.0_rt / dx[0]; @@ -545,9 +542,8 @@ void ca_dergradpy(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/ #if AMREX_SPACEDIM == 3 amrex::Error("3D not supported"); -#endif -#if AMREX_SPACEDIM == 1 - amrex::Error("1D not supported"); +#elif AMREX_SPACEDIM == 1 + return; // Skip for 1D #endif Real dxinv = 1.0_rt / dx[0]; @@ -658,9 +654,8 @@ void ca_dergradrhooverrho(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int / #if AMREX_SPACEDIM == 3 amrex::Error("3D not supported"); -#endif -#if AMREX_SPACEDIM == 1 - amrex::Error("1D not supported"); +#elif AMREX_SPACEDIM == 1 + return; // Skip for 1D #endif Real dxinv = 1.0_rt / dx[0]; @@ -719,4 +714,4 @@ void ca_dergradrhooverrho(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int / }); -} +} \ No newline at end of file diff --git a/Exec/science/subch_planar/problem_initialize_state_data.H b/Exec/science/subch_planar/problem_initialize_state_data.H index f002f1339c..48f0aceac3 100644 --- a/Exec/science/subch_planar/problem_initialize_state_data.H +++ b/Exec/science/subch_planar/problem_initialize_state_data.H @@ -25,7 +25,9 @@ void problem_initialize_state_data (int i, int j, int k, Real z = problo[2] + dx[2] * (static_cast(k) + 0.5_rt); #endif -#if AMREX_SPACEDIM == 2 +#if AMREX_SPACEDIM == 1 + Real height = x; +#elif AMREX_SPACEDIM == 2 Real height = y; #else Real height = z; diff --git a/Exec/science/wdmerger/Prob.cpp b/Exec/science/wdmerger/Prob.cpp index ac6a41eceb..ba28a87a18 100644 --- a/Exec/science/wdmerger/Prob.cpp +++ b/Exec/science/wdmerger/Prob.cpp @@ -604,7 +604,7 @@ Castro::update_relaxation(Real time, Real dt) { const Real ldt = new_time - old_time; - force[lev].reset(new MultiFab(getLevel(lev).grids, getLevel(lev).dmap, NUM_STATE, 0)); + force[lev] = std::make_unique(getLevel(lev).grids, getLevel(lev).dmap, NUM_STATE, 0); force[lev]->setVal(0.0); MultiFab& S_new = getLevel(lev).get_new_data(State_Type); diff --git a/Exec/science/wdmerger/Problem.H b/Exec/science/wdmerger/Problem.H index db48b224e7..8696516ab6 100644 --- a/Exec/science/wdmerger/Problem.H +++ b/Exec/science/wdmerger/Problem.H @@ -27,11 +27,13 @@ void volInBoundary (amrex::Real time, amrex::Real& vol_p, amrex::Real& vol_s, am // Computes standard dot product of two three-vectors. +static amrex::Real dot_product(const amrex::Real a[], const amrex::Real b[]); // Computes norm of a three vector. +static amrex::Real norm(const amrex::Real a[]); // Problem post-initialization routine. @@ -48,6 +50,7 @@ void problem_post_timestep(); // Write out the git hashes for the various parts of the code. +static void writeGitHashes(std::ostream& log); // Update relaxation process. diff --git a/Exec/science/wdmerger/problem_checkpoint.H b/Exec/science/wdmerger/problem_checkpoint.H index 1ddee5e049..9ee3b16846 100644 --- a/Exec/science/wdmerger/problem_checkpoint.H +++ b/Exec/science/wdmerger/problem_checkpoint.H @@ -7,7 +7,7 @@ #include AMREX_INLINE -void problem_checkpoint (std::string checkpoint_dir) +void problem_checkpoint (const std::string& checkpoint_dir) { std::ofstream com; com.open(checkpoint_dir + "/COM", std::ios::out); diff --git a/Exec/science/wdmerger/problem_restart.H b/Exec/science/wdmerger/problem_restart.H index 2246732b64..3acbe99165 100644 --- a/Exec/science/wdmerger/problem_restart.H +++ b/Exec/science/wdmerger/problem_restart.H @@ -5,7 +5,7 @@ #include AMREX_INLINE -void problem_restart (std::string checkpoint_dir) +void problem_restart (const std::string& checkpoint_dir) { std::ifstream com; com.open(checkpoint_dir + "/COM", std::ios::in); diff --git a/Exec/science/wdmerger/wdmerger_util.cpp b/Exec/science/wdmerger/wdmerger_util.cpp index 65b09f89e5..6cb481b5e4 100644 --- a/Exec/science/wdmerger/wdmerger_util.cpp +++ b/Exec/science/wdmerger/wdmerger_util.cpp @@ -658,7 +658,7 @@ void binary_setup () Real v_P_r, v_S_r, v_P_phi, v_S_phi; - kepler_third_law(problem::radius_P, problem::mass_P, problem::radius_S, problem::mass_S, + kepler_third_law(problem::radius_P, problem::mass_P, problem::radius_S, problem::mass_S, // NOLINT(readability-suspicious-call-argument) castro::rotational_period, problem::orbital_eccentricity, problem::orbital_angle, problem::a, problem::r_P_initial, problem::r_S_initial, v_P_r, v_S_r, v_P_phi, v_S_phi); diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 0600a16d1c..d3fe48b979 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -737,7 +737,7 @@ Castro::Castro (Amr& papa, } radiation->regrid(level, grids, dmap); - rad_solver.reset(new RadSolve(parent, level, grids, dmap)); + rad_solver = std::make_unique(parent, level, grids, dmap); } #endif @@ -751,7 +751,7 @@ Castro::Castro (Amr& papa, Castro::~Castro () // NOLINT(modernize-use-equals-default) { #ifdef RADIATION - if (radiation != 0) { + if (radiation != nullptr) { //radiation->cleanup(level); radiation->close(level); } diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index 4202132493..016215aa32 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -71,9 +71,15 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord) if (coord != CoordSys::cartesian) { return false; - } else { - return true; } + return true; + + } else if (mom_dir == 1 && flux_dir == 1) { + + if (coord == CoordSys::SPHERICAL) { + return false; + } + return true; } else { return (mom_dir == flux_dir); @@ -157,10 +163,12 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE Real volume(const int& i, const int& j, const int& k, const GeometryData& geomdata) { + // // Given 3D cell-centered indices (i,j,k), return the volume of the zone. - // Note that Castro has no support for angular coordinates, so - // this function only provides Cartesian in 1D/2D/3D, Cylindrical (R-Z) - // in 2D, and Spherical in 1D. + // this function only provides Cartesian in 1D/2D/3D, + // Cylindrical (R-Z) in 2D, + // and Spherical in 1D and 2D (R-THETA). + // amrex::ignore_unused(i); amrex::ignore_unused(j); @@ -210,8 +218,7 @@ Real volume(const int& i, const int& j, const int& k, vol = dx[0] * dx[1]; - } - else { + } else if (coord == 1) { // Cylindrical @@ -220,6 +227,20 @@ Real volume(const int& i, const int& j, const int& k, vol = M_PI * (r_l + r_r) * dx[0] * dx[1]; + } else { + + // Spherical + + constexpr Real twoThirdsPi = 2.0_rt * M_PI / 3.0_rt; + + Real r_l = geomdata.ProbLo()[0] + static_cast(i) * dx[0]; + Real r_r = geomdata.ProbLo()[0] + static_cast(i+1) * dx[0]; + Real theta_l = geomdata.ProbLo()[1] + static_cast(j) * dx[1]; + Real theta_r = geomdata.ProbLo()[1] + static_cast(j+1) * dx[1]; + + vol = twoThirdsPi * (std::cos(theta_l) - std::cos(theta_r)) * dx[0] * + (r_r * r_r + r_r * r_l + r_l * r_l); + } #else @@ -290,8 +311,7 @@ Real area(const int& i, const int& j, const int& k, a = dx[0]; } - } - else { + } else if (coord == 1) { // Cylindrical @@ -304,6 +324,24 @@ Real area(const int& i, const int& j, const int& k, a = 2.0_rt * M_PI * r * dx[0]; } + } else { + + // Spherical + + if (idir == 0) { + Real r = geomdata.ProbLo()[0] + static_cast(i) * dx[0]; + Real theta_l = geomdata.ProbLo()[1] + static_cast(j) * dx[1]; + Real theta_r = geomdata.ProbLo()[1] + static_cast(j+1) * dx[1]; + + a = 2.0_rt * M_PI * r * r * (std::cos(theta_l) - std::cos(theta_r)); + } + else { + Real r = geomdata.ProbLo()[0] + (static_cast(i) + 0.5_rt) * dx[0]; + Real theta = geomdata.ProbLo()[1] + static_cast(j) * dx[1]; + + a = 2.0_rt * M_PI * std::sin(theta) * r * dx[0]; + } + } #else diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 44f78e6ebf..f11b20719a 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -268,8 +268,9 @@ sdc_extra int 0 # which SDC nonlinear solver to use? 1 = Newton, 2 = VODE, 3 = VODE for first iter sdc_solver int 1 -# for 2-d axisymmetry, do we include the geometry source terms from Bernand-Champmartin? -use_axisymmetric_geom_source bool 1 +# Do we include geometry source terms due to local unit vectors in non-Cartesian Coord? +# We currently support R-Z cylinderical 2D (Bernand-Champmartin) and R-THETA spherical 2D +use_geom_source bool 1 # for simplified-SDC, do we add the reactive source prediction to the interface states # used in the advective source construction? diff --git a/Source/driver/math.H b/Source/driver/math.H index db639a105d..30005d81e7 100644 --- a/Source/driver/math.H +++ b/Source/driver/math.H @@ -3,6 +3,11 @@ #include #include +#include +#include + +using namespace amrex::literals; + AMREX_GPU_HOST_DEVICE AMREX_INLINE void @@ -16,4 +21,11 @@ cross_product(amrex::GpuArray const& a, } + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real cot(amrex::Real x) { + + AMREX_ASSERT(x != 0.0_rt || x != M_PI); + return std::cos(x) / std::sin(x); +} #endif diff --git a/Source/gravity/Gravity.H b/Source/gravity/Gravity.H index 0f6aba88ac..feaf0e23ca 100644 --- a/Source/gravity/Gravity.H +++ b/Source/gravity/Gravity.H @@ -567,12 +567,6 @@ public: const amrex::Vector& res, amrex::Real time); - static inline amrex::Real - get_const_grav() { - return gravity::const_grav; - } - - }; /// diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 18eebefbfc..edeb7359d2 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -881,11 +881,17 @@ Gravity::get_old_grav_vector(int level, MultiFab& grav_vector, Real time) MultiFab grav(grids[level], dmap[level], AMREX_SPACEDIM, ng); grav.setVal(0.0,ng); - if (gravity::gravity_type == "ConstantGrav") { + const Geometry& geom = parent->Geom(level); - // Set to constant value in the AMREX_SPACEDIM direction and zero in all others. + if (gravity::gravity_type == "ConstantGrav") { - grav.setVal(gravity::const_grav,AMREX_SPACEDIM-1,1,ng); + if (AMREX_SPACEDIM == 2 && geom.Coord() == 2) { + // 2D spherical r-theta, we want g in the radial direction + grav.setVal(gravity::const_grav, 0, 1, ng); + } else { + // Set to constant value in the AMREX_SPACEDIM direction and zero in all others. + grav.setVal(gravity::const_grav, AMREX_SPACEDIM-1, 1, ng); + } } else if (gravity::gravity_type == "MonopoleGrav") { @@ -895,7 +901,6 @@ Gravity::get_old_grav_vector(int level, MultiFab& grav_vector, Real time) } else if (gravity::gravity_type == "PoissonGrav") { - const Geometry& geom = parent->Geom(level); amrex::average_face_to_cellcenter(grav, amrex::GetVecOfConstPtrs(grad_phi_prev[level]), geom); grav.mult(-1.0, ng); // g = - grad(phi) @@ -952,11 +957,17 @@ Gravity::get_new_grav_vector(int level, MultiFab& grav_vector, Real time) MultiFab grav(grids[level],dmap[level],AMREX_SPACEDIM,ng); grav.setVal(0.0,ng); + const Geometry& geom = parent->Geom(level); if (gravity::gravity_type == "ConstantGrav") { - // Set to constant value in the AMREX_SPACEDIM direction - grav.setVal(gravity::const_grav,AMREX_SPACEDIM-1,1,ng); + if (AMREX_SPACEDIM == 2 && geom.Coord() == 2) { + // 2D spherical r-theta, we want g in the radial direction + grav.setVal(gravity::const_grav, 0, 1, ng); + } else { + // Set to constant value in the AMREX_SPACEDIM direction + grav.setVal(gravity::const_grav, AMREX_SPACEDIM-1, 1, ng); + } } else if (gravity::gravity_type == "MonopoleGrav") { @@ -967,7 +978,6 @@ Gravity::get_new_grav_vector(int level, MultiFab& grav_vector, Real time) } else if (gravity::gravity_type == "PoissonGrav") { - const Geometry& geom = parent->Geom(level); amrex::average_face_to_cellcenter(grav, amrex::GetVecOfConstPtrs(grad_phi_curr[level]), geom); grav.mult(-1.0, ng); // g = - grad(phi) @@ -2777,12 +2787,12 @@ Gravity::fill_direct_sum_BCs(int crse_level, int fine_level, const Vector::max()); BL_ASSERT(nPtsYZ <= std::numeric_limits::max()); - ParallelDescriptor::ReduceRealSum(bcXYLo.dataPtr(), nPtsXY); - ParallelDescriptor::ReduceRealSum(bcXYHi.dataPtr(), nPtsXY); - ParallelDescriptor::ReduceRealSum(bcXZLo.dataPtr(), nPtsXZ); - ParallelDescriptor::ReduceRealSum(bcXZHi.dataPtr(), nPtsXZ); - ParallelDescriptor::ReduceRealSum(bcYZLo.dataPtr(), nPtsYZ); - ParallelDescriptor::ReduceRealSum(bcYZHi.dataPtr(), nPtsYZ); + ParallelDescriptor::ReduceRealSum(bcXYLo.dataPtr(), static_cast(nPtsXY)); + ParallelDescriptor::ReduceRealSum(bcXYHi.dataPtr(), static_cast(nPtsXY)); + ParallelDescriptor::ReduceRealSum(bcXZLo.dataPtr(), static_cast(nPtsXZ)); + ParallelDescriptor::ReduceRealSum(bcXZHi.dataPtr(), static_cast(nPtsXZ)); + ParallelDescriptor::ReduceRealSum(bcYZLo.dataPtr(), static_cast(nPtsYZ)); + ParallelDescriptor::ReduceRealSum(bcYZHi.dataPtr(), static_cast(nPtsYZ)); #ifdef _OPENMP #pragma omp parallel diff --git a/Source/gravity/binary.H b/Source/gravity/binary.H index 684d875139..1d2438b7d4 100644 --- a/Source/gravity/binary.H +++ b/Source/gravity/binary.H @@ -67,7 +67,7 @@ void lagrange_iterate (Real& r, Real mass_1, Real mass_2, Real r1, Real r2, Real const Real tolerance = 1.0e-8_rt; const int max_iters = 200; - Real rm, rp; + Real rm{}, rp{}; if (r_min == 0.0_rt && r_max == 0.0_rt) { amrex::Abort("Lagrange point iteration must have at least one non-zero bound provided."); diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 7998b7a107..6c222cdddf 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -42,7 +42,7 @@ Castro::consup_hydro(const Box& bx, #endif ) * volinv; - // Add the p div(u) source term to (rho e). + // Add the p div(u) source term to (rho e). if (n == UEINT) { Real pdu = (qx(i+1,j,k,GDPRES) + qx(i,j,k,GDPRES)) * @@ -65,13 +65,19 @@ Castro::consup_hydro(const Box& bx, U_new(i,j,k,n) = U_new(i,j,k,n) - dt * pdu; - } else if (n == UMX) { - // Add gradp term to momentum equation -- only for axisymmetric - // coords (and only for the radial flux). + } else if (n == UMX && !mom_flux_has_p(0, 0, geomdata.Coord())) { + // Add gradp term to momentum equation -- only for axisymmetric + // coords (and only for the radial flux). - if (!mom_flux_has_p(0, 0, geomdata.Coord())) { - U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize()[0]; - } + U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize(0); + +#if AMREX_SPACEDIM >= 2 + } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) { + // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry + + Real r = geomdata.ProbLo(0) + (static_cast(i) + 0.5_rt) * geomdata.CellSize(0); + U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize(1)); +#endif } }); } diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index c70fcd64eb..aa8c9582c5 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -387,6 +387,7 @@ /// @param area1 area of y faces /// @param area2 area of z faces /// @param q0 Godunuv state on x faces +/// @param q1 Godunuv state on y faces /// @param vol cell volume /// void mol_consup(const amrex::Box& bx, @@ -412,6 +413,7 @@ #endif #if AMREX_SPACEDIM <= 2 amrex::Array4 const& q0, + amrex::Array4 const& q1, #endif amrex::Array4 const& vol); diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index fe2004a60c..688b298a70 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -226,6 +226,7 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function #endif #if AMREX_SPACEDIM <= 2 Array4 const& q0, + Array4 const& q1, #endif Array4 const& vol) { @@ -238,6 +239,7 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function #if AMREX_SPACEDIM <= 2 const auto dx = geom.CellSizeArray(); auto coord = geom.Coord(); + auto prob_lo = geom.ProbLoArray(); #endif amrex::ParallelFor(bx, NUM_STATE, @@ -258,12 +260,18 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function #endif #if AMREX_SPACEDIM <= 2 - if (n == UMX && do_hydro == 1) { - // Add gradp term to momentum equation -- only for axisymmetric - // coords (and only for the radial flux). + if (do_hydro == 1) { + if (n == UMX && !mom_flux_has_p(0, 0, coord)) { + // Add gradp term to radial momentum equation -- only for axisymmetric + // coords. - if (!mom_flux_has_p(0, 0, coord)) { update(i,j,k,UMX) -= (q0(i+1,j,k,GDPRES) - q0(i,j,k,GDPRES)) / dx[0]; + + } else if (n == UMY && !mom_flux_has_p(1, 1, coord)) { + // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry + + Real r = prob_lo[0] + (static_cast(i) + 0.5_rt) * dx[0]; + update(i,j,k,UMY) -= (q1(i,j+1,k,GDPRES) - q1(i,j,k,GDPRES)) / (r * dx[1]); } } #endif diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 9fee21334d..75899d67a1 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -620,6 +620,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) #endif #if AMREX_SPACEDIM <= 2 qe[0].array(), + qe[1].array(), #endif volume.array(mfi)); diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 095192931b..869ec477b8 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -205,30 +205,30 @@ Castro::divu(const Box& bx, #if AMREX_SPACEDIM == 1 if (coord_type == 0) { - div(i,j,k) = (q_arr(i,j,k,QU) - q_arr(i-1,j,k,QU)) * dxinv; + div(i,j,k) = (q_arr(i,j,k,QU) - q_arr(i-1,j,k,QU)) * dxinv; } else if (coord_type == 1) { - // axisymmetric - if (i == 0) { - div(i,j,k) = 0.0_rt; - } else { - Real rl = (i - 0.5_rt) * dx[0] + problo[0]; - Real rr = (i + 0.5_rt) * dx[0] + problo[0]; - Real rc = (i) * dx[0] + problo[0]; - - div(i,j,k) = (rr * q_arr(i,j,k,QU) - rl * q_arr(i-1,j,k,QU)) * dxinv / rc; - } + // axisymmetric + if (i == 0) { + div(i,j,k) = 0.0_rt; + } else { + Real rl = (i - 0.5_rt) * dx[0] + problo[0]; + Real rr = (i + 0.5_rt) * dx[0] + problo[0]; + Real rc = (i) * dx[0] + problo[0]; + + div(i,j,k) = (rr * q_arr(i,j,k,QU) - rl * q_arr(i-1,j,k,QU)) * dxinv / rc; + } } else { - // spherical - if (i == 0) { - div(i,j,k) = 0.0_rt; - } else { - Real rl = (i - 0.5_rt) * dx[0] + problo[0]; - Real rr = (i + 0.5_rt) * dx[0] + problo[0]; - Real rc = (i) * dx[0] + problo[0]; - - div(i,j,k) = (rr * rr * q_arr(i,j,k,QU) - rl * rl * q_arr(i-1,j,k,QU)) * dxinv / (rc * rc); - } + // spherical + if (i == 0) { + div(i,j,k) = 0.0_rt; + } else { + Real rl = (i - 0.5_rt) * dx[0] + problo[0]; + Real rr = (i + 0.5_rt) * dx[0] + problo[0]; + Real rc = (i) * dx[0] + problo[0]; + + div(i,j,k) = (rr * rr * q_arr(i,j,k,QU) - rl * rl * q_arr(i-1,j,k,QU)) * dxinv / (rc * rc); + } } #endif @@ -237,31 +237,64 @@ Castro::divu(const Box& bx, Real vy = 0.0_rt; if (coord_type == 0) { - ux = 0.5_rt * (q_arr(i,j,k,QU) - q_arr(i-1,j,k,QU) + q_arr(i,j-1,k,QU) - q_arr(i-1,j-1,k,QU)) * dxinv; - vy = 0.5_rt * (q_arr(i,j,k,QV) - q_arr(i,j-1,k,QV) + q_arr(i-1,j,k,QV) - q_arr(i-1,j-1,k,QV)) * dyinv; + + // Cartesian + + ux = 0.5_rt * (q_arr(i,j,k,QU) - q_arr(i-1,j,k,QU) + q_arr(i,j-1,k,QU) - q_arr(i-1,j-1,k,QU)) * dxinv; + vy = 0.5_rt * (q_arr(i,j,k,QV) - q_arr(i,j-1,k,QV) + q_arr(i-1,j,k,QV) - q_arr(i-1,j-1,k,QV)) * dyinv; + + } else if (coord_type == 1) { + + // Cylindrical R-Z + + if (i == 0) { + ux = 0.0_rt; + vy = 0.0_rt; // is this part correct? + } else { + Real rl = (i - 0.5_rt) * dx[0] + problo[0]; + Real rr = (i + 0.5_rt) * dx[0] + problo[0]; + Real rc = (i) * dx[0] + problo[0]; + + // These are transverse averages in the y-direction + Real ul = 0.5_rt * (q_arr(i-1,j,k,QU) + q_arr(i-1,j-1,k,QU)); + Real ur = 0.5_rt * (q_arr(i,j,k,QU) + q_arr(i,j-1,k,QU)); + + // Take 1/r d/dr(r*u) + ux = (rr * ur - rl * ul) * dxinv / rc; + + // These are transverse averages in the x-direction + Real vb = 0.5_rt * (q_arr(i,j-1,k,QV) + q_arr(i-1,j-1,k,QV)); + Real vt = 0.5_rt * (q_arr(i,j,k,QV) + q_arr(i-1,j,k,QV)); + + vy = (vt - vb) * dyinv; + } } else { - if (i == 0) { - ux = 0.0_rt; - vy = 0.0_rt; // is this part correct? - } else { + + // Spherical R-Theta + Real rl = (i - 0.5_rt) * dx[0] + problo[0]; Real rr = (i + 0.5_rt) * dx[0] + problo[0]; Real rc = (i) * dx[0] + problo[0]; + // cell-centered sin(theta) of top, bot cell and face-centered + Real sint = std::sin((j + 0.5_rt) * dx[1] + problo[1]); + Real sinb = std::sin((j - 0.5_rt) * dx[1] + problo[1]); + Real sinc = std::sin(j * dx[1] + problo[1]); + // These are transverse averages in the y-direction Real ul = 0.5_rt * (q_arr(i-1,j,k,QU) + q_arr(i-1,j-1,k,QU)); Real ur = 0.5_rt * (q_arr(i,j,k,QU) + q_arr(i,j-1,k,QU)); - // Take 1/r d/dr(r*u) - ux = (rr * ur - rl * ul) * dxinv / rc; + // Finite difference to get divergence. ux = 1/r^2 d/dr(r^2 * u) + ux = (ur * rr * rr - ul * rl * rl) * dxinv / (rc * rc); // These are transverse averages in the x-direction Real vb = 0.5_rt * (q_arr(i,j-1,k,QV) + q_arr(i-1,j-1,k,QV)); Real vt = 0.5_rt * (q_arr(i,j,k,QV) + q_arr(i-1,j,k,QV)); - vy = (vt - vb) * dyinv; - } + // Finite difference to get divergence. vy = 1/(r sin) * d/dtheta(v sin) + vy = (vt * sint - vb * sinb) * dyinv / (rc * sinc); } div(i,j,k) = ux + vy; diff --git a/Source/radiation/HypreABec.H b/Source/radiation/HypreABec.H index 10cc5e96f6..b41f641427 100644 --- a/Source/radiation/HypreABec.H +++ b/Source/radiation/HypreABec.H @@ -33,6 +33,10 @@ class HypreABec { ~HypreABec(); + HypreABec(const HypreABec& a) = delete; + HypreABec(const HypreABec&& a) = delete; + HypreABec& operator= (const HypreABec& a) = delete; + HypreABec& operator= (const HypreABec&& a) = delete; /// /// @param v @@ -48,10 +52,10 @@ class HypreABec { /// void setScalars(amrex::Real alpha, amrex::Real beta); - amrex::Real getAlpha() const { + [[nodiscard]] amrex::Real getAlpha() const { return alpha; } - amrex::Real getBeta() const { + [[nodiscard]] amrex::Real getBeta() const { return beta; } diff --git a/Source/radiation/HypreExtMultiABec.H b/Source/radiation/HypreExtMultiABec.H index 911b2ca657..ee217fc1e1 100644 --- a/Source/radiation/HypreExtMultiABec.H +++ b/Source/radiation/HypreExtMultiABec.H @@ -23,12 +23,16 @@ class HypreExtMultiABec : public HypreMultiABec { a2coefs(fine_level+1), ccoefs( fine_level+1), d1coefs(fine_level+1), - d2coefs(fine_level+1), - alpha2(0.0), gamma(0.0), delta1(0.0), delta2(0.0) - { - } + d2coefs(fine_level+1) + {} + + ~HypreExtMultiABec() override; + + HypreExtMultiABec(const HypreExtMultiABec& src) = delete; + HypreExtMultiABec(const HypreExtMultiABec&& src) = delete; - ~HypreExtMultiABec(); + HypreExtMultiABec& operator= (const HypreExtMultiABec& src) = delete; + HypreExtMultiABec& operator= (const HypreExtMultiABec&& src) = delete; amrex::Real& a2Multiplier() { return alpha2; @@ -109,12 +113,12 @@ class HypreExtMultiABec : public HypreMultiABec { amrex::MultiFab& dest, int icomp, amrex::MultiFab& rhs, ///< will not be altered - BC_Mode inhom); + BC_Mode inhom) override; void loadLevelVectorB(int level, amrex::MultiFab& rhs, ///< will not be altered - BC_Mode inhom); + BC_Mode inhom) override; - void loadMatrix(); ///< once all level coeffs and scalars have been set + void loadMatrix() override; ///< once all level coeffs and scalars have been set /// @@ -134,7 +138,7 @@ class HypreExtMultiABec : public HypreMultiABec { amrex::Vector > > ccoefs; ///< face-based, 2 component amrex::Vector > > d1coefs; ///< cell-based but directional amrex::Vector > > d2coefs; ///< face-based - amrex::Real alpha2, gamma, delta1, delta2; ///< multipliers for the above + amrex::Real alpha2{}, gamma{}, delta1{}, delta2{}; ///< multipliers for the above }; #endif diff --git a/Source/radiation/HypreMultiABec.H b/Source/radiation/HypreMultiABec.H index 0a39aec47f..4aec6e89b9 100644 --- a/Source/radiation/HypreMultiABec.H +++ b/Source/radiation/HypreMultiABec.H @@ -21,38 +21,34 @@ class AuxVar { class Connex { public: - Connex() { - other = NULL; - } + Connex() : + other(nullptr) + {} /// /// @param p /// @param r /// - Connex(AuxVar* p, amrex::Real r) { - val = r; - other = p; - } + Connex(AuxVar* p, amrex::Real r) : + val(r), other(p) + {} /// /// @param lev /// @param iv /// @param r /// - Connex(int lev, const amrex::IntVect& iv, amrex::Real r) { - val = r; - index = iv; - level = lev; - other = NULL; - } + Connex(int lev, const amrex::IntVect& iv, amrex::Real r) : + val(r), index(iv), level(lev), other(nullptr) + {} /// /// @param c /// - bool same_target(const Connex& c) { - return ((other != NULL) + [[nodiscard]] bool same_target(const Connex& c) const { + return ((other != nullptr) ? (other == c.other) - : (c.other == NULL && level == c.level && index == c.index)); + : (c.other == nullptr && level == c.level && index == c.index)); } amrex::Real val; @@ -63,8 +59,7 @@ class AuxVar { public: - AuxVar() : secondary_flag(0) { - } + AuxVar() = default; /// @@ -72,7 +67,7 @@ class AuxVar { /// @param r /// void push(AuxVar* p, amrex::Real r) { - a.push_back(Connex(p,r)); + a.emplace_back(p, r); } @@ -82,7 +77,7 @@ class AuxVar { /// @param r /// void push(int lev, const amrex::IntVect& iv, amrex::Real r) { - a.push_back(Connex(lev,iv,r)); + a.emplace_back(lev, iv, r); } @@ -95,7 +90,7 @@ class AuxVar { /// @param p->secondary_flag /// if (p->secondary_flag == 0) { // don't count the same secondary twice - a.push_back(Connex(p,1.0)); + a.emplace_back(p, 1.0); p->secondary_flag = 1; } } @@ -104,7 +99,7 @@ class AuxVar { return a.empty(); } - bool secondary() { + [[nodiscard]] bool secondary() const { return secondary_flag; } @@ -128,7 +123,7 @@ class AuxVar { protected: std::list a; - int secondary_flag; + int secondary_flag{}; }; /// @@ -149,7 +144,7 @@ class AuxVarBox { /// @param bx /// AuxVarBox(const amrex::Box& bx) : domain(bx) { - int numpts = domain.numPts(); + auto numpts = domain.numPts(); dptr = new AuxVar[numpts]; } @@ -157,13 +152,19 @@ class AuxVarBox { delete[] dptr; } + AuxVarBox(const AuxVarBox& src) = delete; + AuxVarBox(const AuxVarBox&& src) = delete; + + AuxVarBox& operator= (const AuxVarBox& src) = delete; + AuxVarBox& operator= (const AuxVarBox&& src) = delete; + AuxVar& operator()(const amrex::IntVect& p) { BL_ASSERT(!(dptr == 0)); BL_ASSERT(domain.contains(p)); return dptr[domain.index(p)]; } - const amrex::Box& box() const { + [[nodiscard]] const amrex::Box& box() const { return domain; } @@ -346,7 +347,9 @@ class CrseBndryAuxVar : public BndryAuxVarBase { /// @param ori /// @param i /// - int size (const amrex::Orientation ori, int i) const { return aux[ori][i].size(); } + int size (const amrex::Orientation ori, int i) const { + return static_cast(aux[ori][i].size()); + } AuxVarBox& operator() (const amrex::Orientation ori, int i, int j) { return *aux[ori][i][j]; @@ -408,6 +411,11 @@ class HypreMultiABec { virtual ~HypreMultiABec(); + HypreMultiABec(const HypreMultiABec& src) = delete; + HypreMultiABec(const HypreMultiABec&& src) = delete; + + HypreMultiABec operator=(const HypreMultiABec& src) = delete; + HypreMultiABec operator=(const HypreMultiABec&& src) = delete; /// /// @param level @@ -422,10 +430,10 @@ class HypreMultiABec { const amrex::DistributionMapping& _dmap, amrex::IntVect _crse_ratio); - int crseLevel() { + [[nodiscard]] int crseLevel() const { return crse_level; } - int fineLevel() { + [[nodiscard]] int fineLevel() const { return fine_level; } @@ -501,10 +509,10 @@ class HypreMultiABec { /// void setScalars(amrex::Real alpha, amrex::Real beta); - amrex::Real getAlpha() const { + [[nodiscard]] amrex::Real getAlpha() const { return alpha; } - amrex::Real getBeta() const { + [[nodiscard]] amrex::Real getBeta() const { return beta; } diff --git a/Source/radiation/MGRadBndry.H b/Source/radiation/MGRadBndry.H index 8b922323d4..a565e2eb43 100644 --- a/Source/radiation/MGRadBndry.H +++ b/Source/radiation/MGRadBndry.H @@ -22,31 +22,28 @@ public: const int _ngroups, const amrex::Geometry& _geom); - ~MGRadBndry(); - - /// /// @param phys_bc /// @param geom /// @param ratio /// - virtual void setBndryConds(const amrex::BCRec& phys_bc, - const amrex::Geometry& geom, amrex::IntVect& ratio); + void setBndryConds(const amrex::BCRec& phys_bc, + const amrex::Geometry& geom, amrex::IntVect& ratio) override; /// /// @param bc /// @param phys_bc_mode /// - virtual void setBndryFluxConds(const amrex::BCRec& bc, - const BC_Mode phys_bc_mode = Inhomogeneous_BC); + void setBndryFluxConds(const amrex::BCRec& bc, + const BC_Mode phys_bc_mode = Inhomogeneous_BC) override; /// /// @param _face /// - virtual int mixedBndry(const amrex::Orientation& _face) const { - return (bcflag[_face] > 1) ? 1 : 0; + int mixedBndry(const amrex::Orientation& _face) const override { + return (bcflag[_face] > 1) ? 1 : 0; } @@ -99,7 +96,7 @@ public: NGBndry* operator()(const amrex::BoxArray& _grids, const amrex::DistributionMapping& _dmap, int _ngroups, - const amrex::Geometry& _geom) const { + const amrex::Geometry& _geom) const override { /// /// @param _grids diff --git a/Source/radiation/MGRadBndry.cpp b/Source/radiation/MGRadBndry.cpp index 73c4abf41a..32ab3824c9 100644 --- a/Source/radiation/MGRadBndry.cpp +++ b/Source/radiation/MGRadBndry.cpp @@ -58,10 +58,6 @@ MGRadBndry::MGRadBndry(const BoxArray& _grids, } } -MGRadBndry::~MGRadBndry() -{ -} - void MGRadBndry::init(const int _ngroups) { // obsolete implementation of the Marshak boundary condition requires diff --git a/Source/radiation/NGBndry.H b/Source/radiation/NGBndry.H index 0fccdaa3c3..e90c1faef2 100644 --- a/Source/radiation/NGBndry.H +++ b/Source/radiation/NGBndry.H @@ -22,16 +22,10 @@ class NGBndry : public RadInterpBndryData public: NGBndry(const amrex::BoxArray& _grids, const amrex::DistributionMapping& _dmap, int _ncomp, const amrex::Geometry& _geom) : - -/// -/// @param _grids -/// @param _dmap -/// @param _ncomp -/// @param _geom -/// RadInterpBndryData(_grids,_dmap,_ncomp,_geom) { } + /// /// @param bc /// @param phys_bc_mode @@ -51,7 +45,7 @@ public: /// /// @param _face /// - virtual int mixedBndry(const amrex::Orientation& _face) const { + virtual int mixedBndry(const amrex::Orientation& /* _face */) const { return 0; } @@ -62,13 +56,6 @@ protected: /// amrex::Vector< std::unique_ptr > > bctypearray[2*AMREX_SPACEDIM]; - -/// -/// @param src -/// -private: - NGBndry(const NGBndry& src); - NGBndry& operator=(const NGBndry& src); }; /// diff --git a/Source/radiation/RadBndry.H b/Source/radiation/RadBndry.H index 2f137f4e72..4781ee6c11 100644 --- a/Source/radiation/RadBndry.H +++ b/Source/radiation/RadBndry.H @@ -28,30 +28,28 @@ public: RadBndry(const amrex::BoxArray& _grids, const amrex::DistributionMapping& _dmap, const amrex::Geometry& _geom, amrex::Real bv); - ~RadBndry(); - /// /// @param phys_bc /// @param geom /// @param ratio /// - virtual void setBndryConds(const amrex::BCRec& phys_bc, - const amrex::Geometry& geom, amrex::IntVect& ratio); + void setBndryConds(const amrex::BCRec& phys_bc, + const amrex::Geometry& geom, amrex::IntVect& ratio) override; /// /// @param bc /// @param phys_bc_mode /// - virtual void setBndryFluxConds(const amrex::BCRec& bc, - const BC_Mode phys_bc_mode = Inhomogeneous_BC); + void setBndryFluxConds(const amrex::BCRec& bc, + const BC_Mode phys_bc_mode = Inhomogeneous_BC) override; /// /// @param _face /// - virtual int mixedBndry(const amrex::Orientation& _face) const { + int mixedBndry(const amrex::Orientation& _face) const override { return (bcflag[_face] > 1) ? 1 : 0; } @@ -110,8 +108,8 @@ public: /// NGBndry* operator()(const amrex::BoxArray& _grids, const amrex::DistributionMapping& _dmap, - int _ncomp, - const amrex::Geometry& _geom) const { + int /* _ncomp */, + const amrex::Geometry& _geom) const override { return new RadBndry(_grids, _dmap, _geom); } }; diff --git a/Source/radiation/RadBndry.cpp b/Source/radiation/RadBndry.cpp index 3635f8a512..aebb5b39fc 100644 --- a/Source/radiation/RadBndry.cpp +++ b/Source/radiation/RadBndry.cpp @@ -65,10 +65,6 @@ RadBndry::RadBndry(const BoxArray& _grids, const DistributionMapping& _dmap, setBndryValues(bv); } -RadBndry::~RadBndry() -{ -} - void RadBndry::init() { // obsolete implementation of the Marshak boundary condition requires diff --git a/Source/radiation/RadSolve.H b/Source/radiation/RadSolve.H index d85c6dee43..2fb694ab03 100644 --- a/Source/radiation/RadSolve.H +++ b/Source/radiation/RadSolve.H @@ -29,7 +29,6 @@ class RadSolve { RadSolve (amrex::Amr* Parent, int level, const amrex::BoxArray& grids, const amrex::DistributionMapping& dmap); - ~RadSolve () {} /// /// query runtime parameters diff --git a/Source/radiation/Radiation.H b/Source/radiation/Radiation.H index 58be11061a..7aff7e5bc5 100644 --- a/Source/radiation/Radiation.H +++ b/Source/radiation/Radiation.H @@ -111,8 +111,6 @@ public: /// @param restart /// Radiation(amrex::Amr* Parent, class Castro* castro, int restart = 0); - ~Radiation() { } - /// /// @param level @@ -193,11 +191,11 @@ public: return delta_T_rat_level[lev]; } - amrex::Real deltaEnergyTol() { + [[nodiscard]] amrex::Real deltaEnergyTol() const { return delta_e_rat_dt_tol; } - amrex::Real deltaTTol() { + [[nodiscard]] amrex::Real deltaTTol() const { return delta_T_rat_dt_tol; } diff --git a/Source/radiation/_interpbndry/RadBndryData.H b/Source/radiation/_interpbndry/RadBndryData.H index 3df1041532..c58eab3d9b 100644 --- a/Source/radiation/_interpbndry/RadBndryData.H +++ b/Source/radiation/_interpbndry/RadBndryData.H @@ -94,11 +94,25 @@ public: //@ManMemo: administrative functions //@ManDoc: default constructor RadBndryData() : amrex::BndryRegister() {}; + + ~RadBndryData() = default; + + // + // Disabled! + // +//@ManDoc: copy constructor + RadBndryData(const RadBndryData& src) = delete; +//@ManDoc: copy operator + RadBndryData& operator = (const RadBndryData& src) = delete; +//@ManDoc: move constructor + RadBndryData(const RadBndryData&& src) = delete; +//@ManDoc: move operator + RadBndryData& operator = (const RadBndryData&& src) = delete; + //@ManDoc: constructor specifying number of components and box of physical domain (cell-centered) RadBndryData(const amrex::BoxArray& _grids, const amrex::DistributionMapping& _dmap, int _ncomp, const ProxyGeometry& geom); -//@ManDoc: destructor - virtual ~RadBndryData(); + //@ManDoc: allocate bndry fabs along given face void define(const amrex::BoxArray& _grids, const amrex::DistributionMapping& _dmap, int _ncomp, const ProxyGeometry& geom); @@ -168,15 +182,6 @@ public: } -private: - // - // Disabled! - // -//@ManDoc: copy constructor - RadBndryData(const RadBndryData& src); -//@ManDoc: copy operator - RadBndryData& operator = (const RadBndryData& src); }; #endif - diff --git a/Source/radiation/_interpbndry/RadBndryData.cpp b/Source/radiation/_interpbndry/RadBndryData.cpp index 850a36e797..0439f707dc 100644 --- a/Source/radiation/_interpbndry/RadBndryData.cpp +++ b/Source/radiation/_interpbndry/RadBndryData.cpp @@ -32,10 +32,6 @@ RadBndryData::RadBndryData(const BoxArray& _grids, const DistributionMapping& _d // (*this) = src; // } -RadBndryData::~RadBndryData() -{ -} - std::ostream& operator << (std::ostream& os, const RadBndryData &mgb) { const BoxArray& grds = mgb.boxes(); @@ -129,4 +125,3 @@ RadBndryData::define(const BoxArray& _grids, const DistributionMapping& _dmap, } } } - diff --git a/Source/sources/Castro_geom.cpp b/Source/sources/Castro_geom.cpp index e3c92732e0..334672edbc 100644 --- a/Source/sources/Castro_geom.cpp +++ b/Source/sources/Castro_geom.cpp @@ -1,31 +1,43 @@ -#include "Castro.H" +#include +#include using namespace amrex; /// -/// this adds the geometric source term for axisymmetric -/// coordinates as described in Bernand-Champmartin. This only -/// applies to axisymmetric geometry. +/// this adds the geometric source terms for non-Cartesian Coordinates +/// This includes 2D Cylindrical (R-Z) coordinate as described in Bernand-Champmartin +/// and 2D Spherical (R-THETA) coordinate. /// + void Castro::construct_old_geom_source(MultiFab& source, MultiFab& state_in, Real time, Real dt) { - if (geom.Coord() != 1) { + amrex::ignore_unused(source); + amrex::ignore_unused(state_in); + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + + if (geom.Coord() == 0) { return; } - if (use_axisymmetric_geom_source == 0) { + if (use_geom_source == 0) { return; } +#if AMREX_SPACEDIM == 2 const Real strt_time = ParallelDescriptor::second(); MultiFab geom_src(grids, dmap, source.nComp(), 0); geom_src.setVal(0.0); - fill_geom_source(time, dt, state_in, geom_src); + if (geom.Coord() == 1) { + fill_RZ_geom_source(time, dt, state_in, geom_src); + } else { + fill_RTheta_geom_source(time, dt, state_in, geom_src); + } Real mult_factor = 1.0; @@ -47,22 +59,31 @@ Castro::construct_old_geom_source(MultiFab& source, MultiFab& state_in, Real tim #endif } +#endif } void -Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFab& state_new, Real time, Real dt) +Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, + MultiFab& state_new, Real time, Real dt) { - if (geom.Coord() != 1) { + amrex::ignore_unused(source); + amrex::ignore_unused(state_old); + amrex::ignore_unused(state_new); + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + + if (geom.Coord() == 0) { return; } - if (use_axisymmetric_geom_source == 0) { + if (use_geom_source == 0) { return; } +#if AMREX_SPACEDIM == 2 const Real strt_time = ParallelDescriptor::second(); MultiFab geom_src(grids, dmap, source.nComp(), 0); @@ -72,7 +93,11 @@ Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFa // Subtract off the old-time value first Real old_time = time - dt; - fill_geom_source(old_time, dt, state_old, geom_src); + if (geom.Coord() == 1) { + fill_RZ_geom_source(old_time, dt, state_old, geom_src); + } else { + fill_RTheta_geom_source(old_time, dt, state_old, geom_src); + } Real mult_factor = -0.5; @@ -84,7 +109,11 @@ Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFa mult_factor = 0.5; - fill_geom_source(time, dt, state_new, geom_src); + if (geom.Coord() == 1) { + fill_RZ_geom_source(time, dt, state_new, geom_src); + } else { + fill_RTheta_geom_source(time, dt, state_new, geom_src); + } MultiFab::Saxpy(source, mult_factor, geom_src, 0, 0, source.nComp(), 0); // NOLINT(readability-suspicious-call-argument) @@ -104,20 +133,22 @@ Castro::construct_new_geom_source(MultiFab& source, MultiFab& state_old, MultiFa #endif } - +#endif } void -Castro::fill_geom_source ([[maybe_unused]] Real time, [[maybe_unused]] Real dt, - MultiFab& cons_state, MultiFab& geom_src) +Castro::fill_RZ_geom_source (Real time, Real dt, MultiFab& cons_state, MultiFab& geom_src) { - // Compute the geometric source for axisymmetric coordinates + // Compute the geometric source for axisymmetric coordinates (R-Z) // resulting from taking the divergence of (rho U U) in cylindrical // coordinates. See the paper by Bernard-Champmartin + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + auto dx = geom.CellSizeArray(); auto prob_lo = geom.ProbLoArray(); @@ -148,3 +179,53 @@ Castro::fill_geom_source ([[maybe_unused]] Real time, [[maybe_unused]] Real dt, }); } } + + + +void +Castro::fill_RTheta_geom_source (Real time, Real dt, MultiFab& cons_state, MultiFab& geom_src) +{ + + // Compute the geometric source resulting from taking the divergence of (rho U U) + // in Spherical 2D (R-Theta) coordinate. + + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + + auto dx = geom.CellSizeArray(); + auto prob_lo = geom.ProbLoArray(); + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(geom_src, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box& bx = mfi.tilebox(); + + Array4 const U_arr = cons_state.array(mfi); + + Array4 const src = geom_src.array(mfi); + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + + // Cell-centered Spherical Radius and Theta + Real r = prob_lo[0] + (static_cast(i) + 0.5_rt)*dx[0]; + Real theta = prob_lo[1] + (static_cast(j) + 0.5_rt)*dx[1]; + + // radial momentum: F = rho (v_theta**2 + v_phi**2) / r + src(i,j,k,UMX) = (U_arr(i,j,k,UMY) * U_arr(i,j,k,UMY) + + U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ)) / (U_arr(i,j,k,URHO) * r); + + // Theta momentum F = rho v_phi**2 cot(theta) / r - rho v_r v_theta / r + src(i,j,k,UMY) = (U_arr(i,j,k,UMZ) * U_arr(i,j,k,UMZ) * cot(theta) - + U_arr(i,j,k,UMX) * U_arr(i,j,k,UMY)) / (U_arr(i,j,k,URHO) * r); + + // Phi momentum: F = - rho v_r v_phi / r - rho v_theta v_phi cot(theta) / r + src(i,j,k,UMZ) = (- U_arr(i,j,k,UMY) * U_arr(i,j,k,UMZ) * cot(theta) - + U_arr(i,j,k,UMX) * U_arr(i,j,k,UMZ)) / (U_arr(i,j,k,URHO) * r); + + }); + } +} diff --git a/Source/sources/Castro_sources.H b/Source/sources/Castro_sources.H index 2381d4e82d..3fab4399c0 100644 --- a/Source/sources/Castro_sources.H +++ b/Source/sources/Castro_sources.H @@ -330,15 +330,27 @@ /// -/// Fill ``ext_src`` with axisymmetric geometry sources +/// Fill ``ext_src`` with 2D cylindrical R-Z geometry sources /// /// @param time current time /// @param dt timestep /// @param S state /// @param ext_src MultiFab to fill with sources /// - void fill_geom_source(amrex::Real time, amrex::Real dt, - amrex::MultiFab& cons_state, amrex::MultiFab& geom_src); + void fill_RZ_geom_source(amrex::Real time, amrex::Real dt, + amrex::MultiFab& cons_state, amrex::MultiFab& geom_src); + + +/// +/// Fill ``ext_src`` with 2D spherical R-Theta geometry sources +/// +/// @param time current time +/// @param dt timestep +/// @param S state +/// @param ext_src MultiFab to fill with sources +/// + void fill_RTheta_geom_source(amrex::Real time, amrex::Real dt, + amrex::MultiFab& cons_state, amrex::MultiFab& geom_src); /// diff --git a/external/Microphysics b/external/Microphysics index 14b8b0e317..3961b439ba 160000 --- a/external/Microphysics +++ b/external/Microphysics @@ -1 +1 @@ -Subproject commit 14b8b0e3173041968943d4bbac2c4803a33abceb +Subproject commit 3961b439ba5c8193975018300e4c23a09b533206 diff --git a/external/amrex b/external/amrex index ac5dde35b6..74127d6d8f 160000 --- a/external/amrex +++ b/external/amrex @@ -1 +1 @@ -Subproject commit ac5dde35b6c10f5d91e289edeff218bde84878a4 +Subproject commit 74127d6d8fa83f922069a25e7ef9f153aa73f68c