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/science/flame_wave/ci-benchmarks/job_info_params.txt b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt index c5c1a3f2af..f3b20994e7 100644 --- a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt +++ b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt @@ -95,7 +95,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/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/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..e7b8751d12 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -157,10 +157,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 +212,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 +221,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 +305,7 @@ Real area(const int& i, const int& j, const int& k, a = dx[0]; } - } - else { + } else if (coord == 1) { // Cylindrical @@ -304,6 +318,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 a7bed8b17e..4983d32b86 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -272,8 +272,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 c5f7bac162..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) 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); ///