From dbd972b75fb437a20e6d593fb4e0dc4de28265f1 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 15:03:29 -0500 Subject: [PATCH 01/10] subtract off center in position and create separate distance function --- Source/driver/Castro.cpp | 19 ++--- Source/driver/Castro_util.H | 25 +++++- Source/driver/Derive.cpp | 116 ++++++--------------------- Source/driver/sum_utils.cpp | 4 - Source/gravity/Castro_gravity.cpp | 8 +- Source/gravity/Gravity.cpp | 28 ++----- Source/hydro/Castro_ctu_hydro.cpp | 3 - Source/hydro/Castro_hybrid.cpp | 23 ++---- Source/hydro/Castro_mol_hydro.cpp | 3 - Source/hydro/advection_util.cpp | 4 - Source/hydro/hybrid.H | 6 +- Source/reactions/Castro_react.cpp | 27 ++----- Source/rotation/Rotation.H | 8 -- Source/rotation/Rotation.cpp | 13 +-- Source/rotation/rotation_sources.cpp | 13 --- Source/scf/scf_relax.cpp | 38 ++------- Source/sources/Castro_sponge.cpp | 18 +---- 17 files changed, 86 insertions(+), 270 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 6149e7f988..c66ff0f6e1 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -3641,19 +3641,9 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time const Real* probhi = geomdata.ProbHi(); const Real* dx = geomdata.CellSize(); - Real loc[3] = {0.0}; + GpuArray loc; - loc[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0]; -#if AMREX_SPACEDIM >= 2 - loc[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1]; -#endif -#if AMREX_SPACEDIM == 3 - loc[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2]; -#endif - - Real r = std::sqrt((loc[0] - problem::center[0]) * (loc[0] - problem::center[0]) + - (loc[1] - problem::center[1]) * (loc[1] - problem::center[1]) + - (loc[2] - problem::center[2]) * (loc[2] - problem::center[2])); + Real r = distance(geomdata, loc); Real max_dist_lo = 0.0; Real max_dist_hi = 0.0; @@ -4357,9 +4347,12 @@ Castro::define_new_center(const MultiFab& S, Real time) // Now broadcast to everyone else. ParallelDescriptor::Bcast(&problem::center[0], AMREX_SPACEDIM, owner); - // Make sure if R-Z that center stays exactly on axis + // Make sure if R-Z and SPHERICAL that center stays exactly on axis if ( Geom().IsRZ() ) { problem::center[0] = 0; + } else if ( Geom().IsSPHERICAL() ) { + problem::center[0] = 0; + problem::center[1] = 0; } } diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index 016215aa32..a9f3a87fae 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -94,7 +94,6 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord) } - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void position(int i, int j, int k, GeometryData const& geomdata, GpuArray& loc, @@ -147,7 +146,7 @@ void position(int i, int j, int k, } for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] = offset[dir] + idx[dir] * dx[dir]; + loc[dir] = offset[dir] + idx[dir] * dx[dir] - problem::center[dir]; } for (int dir = AMREX_SPACEDIM; dir < 3; ++dir) { @@ -156,6 +155,28 @@ void position(int i, int j, int k, } + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +Real distance(GeometryData const& geomdata, GpuArray& loc) +{ + // Returns the distance from the center provided with loc, the position. + // loc is the form of [x,y,z] in Cartesian, [r, z, phi] in cylindrical + // and [r, theta, phi] in spherical + + const auto coord = geomdata.Coord(); + + if (coord == CoordSys::cartesian) { + return std::sqrt(loc[0]*loc[0] + loc[1]*loc[1] + loc[2]*loc[2]); + } + + if (coord == CoordSys::RZ) { + return std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]); + } + + return loc[0]; +} + + namespace geometry_util { diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index 3064fc3dd2..4564db79df 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -584,17 +584,9 @@ extern "C" [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real x = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - Real y = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#else - Real y = 0.0_rt; -#endif -#if AMREX_SPACEDIM == 3 - Real z = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#else - Real z = 0.0_rt; -#endif + GpuArray loc; + + position(i, j, k, geomdata, loc); if (domain_is_plane_parallel) { #if AMREX_SPACEDIM == 2 @@ -607,15 +599,15 @@ extern "C" // where e_r and e_phi are the cylindrical unit vectors // we need the distance in the x-y plane from the origin - Real r = std::sqrt(x*x + y*y); - der(i,j,k,0) = (dat(i,j,k,1)*x + dat(i,j,k,2)*y) / (dat(i,j,k,0)*r); + Real r = std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]); + der(i,j,k,0) = (dat(i,j,k,1)*loc[0] + dat(i,j,k,2)*loc[1]) / (dat(i,j,k,0)*r); #endif } else { - Real r = std::sqrt(x*x + y*y + z*z); + Real r = distance(geomdata, loc); - der(i,j,k,0) = (dat(i,j,k,1)*x + - dat(i,j,k,2)*y + - dat(i,j,k,3)*z) / ( dat(i,j,k,0)*r ); + der(i,j,k,0) = (dat(i,j,k,1)*loc[0] + + dat(i,j,k,2)*loc[1] + + dat(i,j,k,3)*loc[2]) / ( dat(i,j,k,0)*r ); } }); @@ -639,17 +631,9 @@ extern "C" [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real x = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - Real y = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#else - Real y = 0.0_rt; -#endif -#if AMREX_SPACEDIM == 3 - Real z = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#else - Real z = 0.0_rt; -#endif + GpuArray loc; + + position(i, j, k, geomdata, loc); if (domain_is_plane_parallel) { #if AMREX_SPACEDIM == 2 @@ -662,11 +646,11 @@ extern "C" // where e_r and e_phi are the cylindrical unit vectors // we need the distance in the x-y plane from the origin - Real r = std::sqrt(x*x + y*y); - der(i,j,k,0) = (-dat(i,j,k,1)*y + dat(i,j,k,2)*x) / (dat(i,j,k,0)*r); + Real r = std::sqrt(loc[0]*loc[0] + loc[1]*loc[1]); + der(i,j,k,0) = (-dat(i,j,k,1)*loc[1] + dat(i,j,k,2)*loc[0]) / (dat(i,j,k,0)*r); #endif } else { - Real r = std::sqrt(x*x + y*y + z*z); + Real r = distance(geomdata, loc); // we really mean just the velocity component that is // perpendicular to radial, and in general 3-d (e.g. a @@ -676,9 +660,9 @@ extern "C" dat(i,j,k,2)*dat(i,j,k,2) + dat(i,j,k,3)*dat(i,j,k,3))/(dat(i,j,k,0)*dat(i,j,k,0)); - Real vr = (dat(i,j,k,1)*x + - dat(i,j,k,2)*y + - dat(i,j,k,3)*z) / ( dat(i,j,k,0)*r ); + Real vr = (dat(i,j,k,1)*loc[0] + + dat(i,j,k,2)*loc[1] + + dat(i,j,k,3)*loc[2]) / ( dat(i,j,k,0)*r ); der(i,j,k,0) = std::sqrt(amrex::max(vtot2 - vr*vr, 0.0_rt)); } @@ -721,27 +705,8 @@ extern "C" amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real loc[3]; - - //loc calculated like sum_utils.cpp - //This might be equivalent and more modular: position(i, j, k, geomdata, loc); - loc[0] = problo[0] + (0.5_rt + i) * dx[0]; - -#if AMREX_SPACEDIM >= 2 - loc[1] = problo[1] + (0.5_rt + j) * dx[1]; -#else - loc[1] = 0.0_rt; -#endif - -#if AMREX_SPACEDIM == 3 - loc[2] = problo[2] + (0.5_rt + k) * dx[2]; -#else - loc[2] = 0.0_rt; -#endif - - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } + GpuArray loc; + position(i, j, k, geomdata, loc); // Explicitly computing only the required cross-product as in inertial_to_rotational_velocity if (idir == 0) { // cross_product(loc, mom): ang_mom(1)->x) @@ -773,24 +738,8 @@ extern "C" amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real loc[3]; - - loc[0] = problo[0] + (0.5_rt + i) * dx[0]; - -#if AMREX_SPACEDIM >= 2 - loc[1] = problo[1] + (0.5_rt + j) * dx[1]; -#else - loc[1] = 0.0_rt; -#endif - -#if AMREX_SPACEDIM == 3 - loc[2] = problo[2] + (0.5_rt + k) * dx[2]; -#else - loc[2] = 0.0_rt; -#endif - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } + GpuArray loc; + position(i, j, k, geomdata, loc); if (idir == 0) { // cross_product(loc, mom): ang_mom(1)->x) L(i,j,k,0) = loc[1] * dat(i,j,k,3) - loc[2] * dat(i,j,k,2); @@ -821,25 +770,8 @@ extern "C" amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real loc[3]; - - loc[0] = problo[0] + (0.5_rt + i) * dx[0]; - -#if AMREX_SPACEDIM >= 2 - loc[1] = problo[1] + (0.5_rt + j) * dx[1]; -#else - loc[1] = 0.0_rt; -#endif - -#if AMREX_SPACEDIM == 3 - loc[2] = problo[2] + (0.5_rt + k) * dx[2]; -#else - loc[2] = 0.0_rt; -#endif - - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } + GpuArray loc; + position(i, j, k, geomdata, loc); if (idir == 0) { // cross_product(loc, mom): ang_mom(1)->x) L(i,j,k,0) = loc[1] * dat(i,j,k,3) - loc[2] * dat(i,j,k,2); diff --git a/Source/driver/sum_utils.cpp b/Source/driver/sum_utils.cpp index 6423b50b70..82d59bcf9f 100644 --- a/Source/driver/sum_utils.cpp +++ b/Source/driver/sum_utils.cpp @@ -384,10 +384,6 @@ Castro::gwstrain (Real time, GpuArray r; position(i, j, k, geomdata, r); - for (int n = 0; n < 3; ++n) { - r[n] -= problem::center[n]; - } - Real rhoInv; if (rho(i,j,k) * maskFactor > 0.0_rt) { rhoInv = 1.0_rt / rho(i,j,k); diff --git a/Source/gravity/Castro_gravity.cpp b/Source/gravity/Castro_gravity.cpp index 395e0e675e..cbada962bb 100644 --- a/Source/gravity/Castro_gravity.cpp +++ b/Source/gravity/Castro_gravity.cpp @@ -375,10 +375,7 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in, #ifdef HYBRID_MOMENTUM GpuArray loc; - for (int n = 0; n < 3; ++n) { - position(i, j, k, geomdata, loc); - loc[n] -= problem::center[n]; - } + position(i, j, k, geomdata, loc); GpuArray hybrid_src; @@ -574,9 +571,6 @@ void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old, #ifdef HYBRID_MOMENTUM GpuArray loc; position(i, j, k, geomdata, loc); - for (int n = 0; n < 3; ++n) { - loc[n] -= problem::center[n]; - } GpuArray hybrid_src; diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 82003be86f..82bc7e699e 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -1352,22 +1352,8 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab& [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { GpuArray loc; - - loc[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; - -#if AMREX_SPACEDIM >= 2 - loc[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#else - loc[1] = 0.0_rt; -#endif - -#if AMREX_SPACEDIM == 3 - loc[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#else - loc[2] = 0.0_rt; -#endif - - Real r = std::sqrt(loc[0] * loc[0] + loc[1] * loc[1] + loc[2] * loc[2]); + position(i, j, k, geomdata, loc); + Real r = distance(geomdata, loc); int index = static_cast(r / dr); @@ -1497,16 +1483,14 @@ Gravity::compute_radial_mass(const Box& bx, amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xc = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; Real lo_i = problo[0] + static_cast(i) * dx[0] - problem::center[0]; - - Real yc = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; Real lo_j = problo[1] + static_cast(j) * dx[1] - problem::center[1]; - - Real zc = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; Real lo_k = problo[2] + static_cast(k) * dx[2] - problem::center[2]; - Real r = std::sqrt(xc * xc + yc * yc + zc * zc); + GpuArray loc; + position(i, j, k, geomdata, loc); + Real r = distance(geomdata, loc); + int index = static_cast(r * drinv); // We may be coming in here with a masked out zone (in a zone on a coarse diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index d443bbe060..70b2d74008 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -1212,9 +1212,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) - loc[dir] -= problem::center[dir]; - Real R = amrex::max(std::sqrt(loc[0] * loc[0] + loc[1] * loc[1]), std::numeric_limits::min()); Real RInv = 1.0_rt / R; diff --git a/Source/hydro/Castro_hybrid.cpp b/Source/hydro/Castro_hybrid.cpp index 0191587b6f..8dfea244c9 100644 --- a/Source/hydro/Castro_hybrid.cpp +++ b/Source/hydro/Castro_hybrid.cpp @@ -107,9 +107,6 @@ Castro::fill_hybrid_hydro_source(MultiFab& sources, const MultiFab& state_in, Re position(i, j, k, geomdata, loc); - loc[0] -= problem::center[0]; - loc[1] -= problem::center[1]; - Real R = amrex::max(std::sqrt(loc[0] * loc[0] + loc[1] * loc[1]), std::numeric_limits::min()); @@ -150,21 +147,19 @@ Castro::linear_to_hybrid_momentum(MultiFab& state_in, int ng) position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) - loc[dir] -= problem::center[dir]; - GpuArray linear_mom; - for (int dir = 0; dir < 3; ++dir) + for (int dir = 0; dir < 3; ++dir) { linear_mom[dir] = u(i,j,k,UMX+dir); + } GpuArray hybrid_mom; linear_to_hybrid(loc, linear_mom, hybrid_mom); - for (int dir = 0; dir < 3; ++dir) + for (int dir = 0; dir < 3; ++dir) { u(i,j,k,UMR+dir) = hybrid_mom[dir]; - + } }); } } @@ -196,21 +191,19 @@ Castro::hybrid_to_linear_momentum(MultiFab& state_in, int ng) position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) - loc[dir] -= problem::center[dir]; - GpuArray hybrid_mom; - for (int dir = 0; dir < 3; ++dir) + for (int dir = 0; dir < 3; ++dir) { hybrid_mom[dir] = u(i,j,k,UMR+dir); + } GpuArray linear_mom; hybrid_to_linear(loc, hybrid_mom, linear_mom); - for (int dir = 0; dir < 3; ++dir) + for (int dir = 0; dir < 3; ++dir) { u(i,j,k,UMX+dir) = linear_mom[dir]; - + } }); } } diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 75899d67a1..d4f8747a07 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -636,9 +636,6 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) - loc[dir] -= problem::center[dir]; - Real R = amrex::max(std::sqrt(loc[0] * loc[0] + loc[1] * loc[1]), std::numeric_limits::min()); Real RInv = 1.0_rt / R; diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 515778c4d1..f87386d6dd 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -726,10 +726,6 @@ Castro::do_enforce_minimum_density(const Box& bx, position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } - GpuArray linear_mom; for (int dir = 0; dir < 3; ++dir) { diff --git a/Source/hydro/hybrid.H b/Source/hydro/hybrid.H index 0188a1536e..632947f61f 100644 --- a/Source/hydro/hybrid.H +++ b/Source/hydro/hybrid.H @@ -123,13 +123,11 @@ void compute_hybrid_flux(const GpuArray& state, const GeometryData& position(i, j, k, geomdata, loc, ccx, ccy, ccz); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) - loc[dir] -= problem::center[dir]; - GpuArray linear_mom; - for (int dir = 0; dir < 3; ++dir) + for (int dir = 0; dir < 3; ++dir) { linear_mom[dir] = state[GDRHO] * state[GDU + dir]; + } GpuArray hybrid_mom; diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index b7a422f1dd..c04891d531 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -210,6 +210,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra #ifdef MODEL_PARSER const auto problo = geom.ProbLoArray(); #endif + const auto geomdata = goem.data(); #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -270,22 +271,15 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra #ifdef MODEL_PARSER if (drive_initial_convection) { - Real rr[3] = {0.0_rt}; - - rr[0] = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - rr[1] = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; -#endif -#if AMREX_SPACEDIM == 3 - rr[2] = problo[2] + dx[2] * (static_cast(k) + 0.5_rt) - problem::center[2]; -#endif + GpuArray rr; + position(i, j, k, geomdata, rr); Real dist; if (domain_is_plane_parallel) { dist = rr[AMREX_SPACEDIM-1]; } else { - dist = std::sqrt(rr[0] * rr[0] + rr[1] * rr[1] + rr[2] * rr[2]); + dist = distance(geomdata, rr); } burn_state.T_fixed = interpolate(dist, model::itemp); @@ -618,22 +612,15 @@ Castro::react_state(Real time, Real dt) #ifdef MODEL_PARSER if (drive_initial_convection) { - Real rr[3] = {0.0_rt}; - - rr[0] = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - rr[1] = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; -#endif -#if AMREX_SPACEDIM == 3 - rr[2] = problo[2] + dx[2] * (static_cast(k) + 0.5_rt) - problem::center[2]; -#endif + GpuArray rr; + position(i, j, k, geomdata, rr); Real dist; if (domain_is_plane_parallel) { dist = rr[AMREX_SPACEDIM-1]; } else { - dist = std::sqrt(rr[0] * rr[0] + rr[1] * rr[1] + rr[2] * rr[2]); + dist = distance(geomdata, rr); } burn_state.T_fixed = interpolate(dist, model::itemp); diff --git a/Source/rotation/Rotation.H b/Source/rotation/Rotation.H index 878ae4a6b0..ef6c3e4129 100644 --- a/Source/rotation/Rotation.H +++ b/Source/rotation/Rotation.H @@ -129,10 +129,6 @@ inertial_to_rotational_velocity (const int i, const int j, const int k, position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } - auto omega = get_omega(); // do the cross product Omega x loc @@ -159,10 +155,6 @@ rotational_to_inertial_velocity (const int i, const int j, const int k, position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } - auto omega = get_omega(); // do the cross product Omega x loc diff --git a/Source/rotation/Rotation.cpp b/Source/rotation/Rotation.cpp index 9a5f858fec..082d7beff8 100644 --- a/Source/rotation/Rotation.cpp +++ b/Source/rotation/Rotation.cpp @@ -30,18 +30,7 @@ Castro::fill_rotational_psi(const Box& bx, { GpuArray r; - - r[0] = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - r[1] = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; -#else - r[1] = 0.0_rt; -#endif -#if AMREX_SPACEDIM == 3 - r[2] = problo[2] + dx[2] * (static_cast(k) + 0.5_rt) - problem::center[2]; -#else - r[2] = 0.0_rt; -#endif + position(i, j, k, geomdata, r); if (denom != 0.0) { psi(i,j,k) = rotational_potential(r) / denom; diff --git a/Source/rotation/rotation_sources.cpp b/Source/rotation/rotation_sources.cpp index 4c034810d8..e09554d903 100644 --- a/Source/rotation/rotation_sources.cpp +++ b/Source/rotation/rotation_sources.cpp @@ -27,10 +27,6 @@ Castro::rsrc(const Box& bx, GpuArray loc; position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } - Real rho = uold(i,j,k,URHO); Real rhoInv = 1.0_rt / rho; @@ -241,10 +237,6 @@ Castro::corrrsrc(const Box& bx, GpuArray loc; position(i, j, k, geomdata, loc); - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - loc[dir] -= problem::center[dir]; - } - Real rhoo = uold(i,j,k,URHO); Real rhooinv = 1.0_rt / uold(i,j,k,URHO); @@ -448,10 +440,6 @@ Castro::corrrsrc(const Box& bx, position(ie, je, ke, geomdata, loc, ccx, ccy, ccz); - for (int n = 0; n < AMREX_SPACEDIM; ++n) { - loc[n] -= problem::center[n]; - } - GpuArray temp_vel{}; Real temp_Sr[3]; @@ -487,4 +475,3 @@ Castro::corrrsrc(const Box& bx, }); } - diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp index 22935bb4e9..7b1899ee13 100644 --- a/Source/scf/scf_relax.cpp +++ b/Source/scf/scf_relax.cpp @@ -247,15 +247,8 @@ Castro::do_hscf_solve() // The below assumes we are rotating on the z-axis. - Real r[3] = {0.0}; - - r[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - r[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#endif -#if AMREX_SPACEDIM == 3 - r[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#endif + GpuArrayr; + position(i, j, k, geomdata, r); // Do a trilinear interpolation to find the contribution from // this grid point. Limit so that only the nearest zone centers @@ -381,14 +374,7 @@ Castro::do_hscf_solve() // The below assumes we are rotating on the z-axis. GpuArray r = {0.0}; - - r[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - r[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#endif -#if AMREX_SPACEDIM == 3 - r[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#endif + position(i, j, k, geomdata, r); // Do a trilinear interpolation to find the contribution from // this grid point. Limit so that only the nearest zone centers @@ -459,14 +445,7 @@ Castro::do_hscf_solve() const auto *problo = geomdata.ProbLo(); GpuArray r = {0.0}; - - r[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - r[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#endif -#if AMREX_SPACEDIM == 3 - r[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#endif + position(i, j, k, geomdata, r); enthalpy_arr(i,j,k) = bernoulli - phi_arr(i,j,k) - rotational_potential(r); }); @@ -642,14 +621,7 @@ Castro::do_hscf_solve() const auto* problo = geomdata.ProbLo(); GpuArray r = {0.0}; - - r[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; -#if AMREX_SPACEDIM >= 2 - r[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#endif -#if AMREX_SPACEDIM == 3 - r[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#endif + position(i, j, k, geomdata, r); if (state_arr(i,j,k,URHO) > 0.0) { diff --git a/Source/sources/Castro_sponge.cpp b/Source/sources/Castro_sponge.cpp index 2109d63374..ebf0fef625 100644 --- a/Source/sources/Castro_sponge.cpp +++ b/Source/sources/Castro_sponge.cpp @@ -81,6 +81,7 @@ Castro::apply_sponge(const Box& bx, auto dx = geom.CellSizeArray(); auto problo = geom.ProbLoArray(); + const auto geomdata = geom.data(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -89,20 +90,7 @@ Castro::apply_sponge(const Box& bx, Real src[NSRC] = {0.0}; GpuArray r; - - r[0] = problo[0] + (static_cast(i) + 0.5_rt) * dx[0] - problem::center[0]; - -#if AMREX_SPACEDIM >= 2 - r[1] = problo[1] + (static_cast(j) + 0.5_rt) * dx[1] - problem::center[1]; -#else - r[1] = 0.0_rt; -#endif - -#if AMREX_SPACEDIM == 3 - r[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; -#else - r[2] = 0.0_rt; -#endif + position(i, j, k, geomdata, r); Real rho = state_in(i,j,k,URHO); Real rhoInv = 1.0_rt / rho; @@ -124,7 +112,7 @@ Castro::apply_sponge(const Box& bx, Real sponge_factor = 0.0_rt; if (sponge_lower_radius >= 0.0_rt && sponge_upper_radius > sponge_lower_radius) { - Real rad = std::sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]); + Real rad = distance(geomdata, r); if (rad < sponge_lower_radius) { sponge_factor = sponge_lower_factor; From ec20cb7bc5d2200c859af218f45ff82e140510da Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 15:57:43 -0500 Subject: [PATCH 02/10] fix compile --- Source/driver/Castro.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index c66ff0f6e1..d9900c479d 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -3639,10 +3639,9 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time { const Real* problo = geomdata.ProbLo(); const Real* probhi = geomdata.ProbHi(); - const Real* dx = geomdata.CellSize(); GpuArray loc; - + position(i, j, k, geomdata, loc); Real r = distance(geomdata, loc); Real max_dist_lo = 0.0; From ac9568bbc1815c600bbe2cf342c7d19ea1bcbcd2 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 16:18:01 -0500 Subject: [PATCH 03/10] name variable to geom when its amrex::Geometry to avoid confusion --- Source/driver/Derive.H | 70 ++++++++++++------------ Source/driver/Derive.cpp | 114 +++++++++++++++++++-------------------- 2 files changed, 92 insertions(+), 92 deletions(-) diff --git a/Source/driver/Derive.H b/Source/driver/Derive.H index 7dd302da16..0533c2d824 100644 --- a/Source/driver/Derive.H +++ b/Source/driver/Derive.H @@ -12,125 +12,125 @@ extern "C" void ca_derpres (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dereint1 (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dereint2 (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derlogden (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_deruplusc (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_deruminusc (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dersoundspeed (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dergamma1 (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermachnumber (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derentropy (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #ifdef DIFFUSION void ca_dercond (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdiffcoeff (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdiffterm (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #endif void ca_derenuc (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derenuctimescale (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dervel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagvel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermaggrav (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derradialvel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dercircvel (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagmom (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derangmomx (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_derangmomy (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_derangmomz (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_derkineng (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real time, const int* bcrec, int level); void ca_dernull @@ -144,53 +144,53 @@ extern "C" void ca_derspec (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derabar (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derye (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagvort (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdivu (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derstate (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #ifdef MHD void ca_dermagcenx (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagceny (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_dermagcenz (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); void ca_derdivb (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/, - const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata, + const amrex::FArrayBox& datfab, const amrex::Geometry& geom, amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/); #endif diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index 3064fc3dd2..63d31449d1 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -21,7 +21,7 @@ extern "C" // need to explicitly synchronize after GPU kernels. void ca_derpres(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -54,7 +54,7 @@ extern "C" } void ca_dereint1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -76,7 +76,7 @@ extern "C" } void ca_dereint2(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -92,7 +92,7 @@ extern "C" } void ca_derlogden(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -107,7 +107,7 @@ extern "C" } void ca_deruplusc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -142,7 +142,7 @@ extern "C" } void ca_deruminusc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -177,7 +177,7 @@ extern "C" } void ca_dersoundspeed(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -213,7 +213,7 @@ extern "C" void ca_dergamma1(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -248,7 +248,7 @@ extern "C" } void ca_dermachnumber(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -286,7 +286,7 @@ extern "C" } void ca_derentropy(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -321,7 +321,7 @@ extern "C" #ifdef DIFFUSION void ca_dercond(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -333,7 +333,7 @@ extern "C" } void ca_derdiffcoeff(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -345,7 +345,7 @@ extern "C" } void ca_derdiffterm(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -360,9 +360,9 @@ extern "C" fill_temp_cond(obx, dat, coeff_arr); - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); - const int coord_type = geomdata.Coord(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); + const int coord_type = geom.Coord(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -432,14 +432,14 @@ extern "C" #ifdef REACTIONS void ca_derenuctimescale(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); Real dd = 0.0_rt; #if AMREX_SPACEDIM == 1 @@ -494,7 +494,7 @@ extern "C" } void ca_derenuc(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); @@ -512,7 +512,7 @@ extern "C" #endif void ca_dervel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -528,7 +528,7 @@ extern "C" void ca_dermagvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -549,7 +549,7 @@ extern "C" void ca_dermaggrav(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -568,7 +568,7 @@ extern "C" } void ca_derradialvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -577,8 +577,8 @@ extern "C" auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -623,7 +623,7 @@ extern "C" void ca_dercircvel(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -632,8 +632,8 @@ extern "C" auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -688,7 +688,7 @@ extern "C" void ca_dermagmom(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -707,13 +707,13 @@ extern "C" } void ca_derangmomx (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { int idir = 0; - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); auto const dat = datfab.array(); auto const L = derfab.array(); @@ -724,7 +724,7 @@ extern "C" Real loc[3]; //loc calculated like sum_utils.cpp - //This might be equivalent and more modular: position(i, j, k, geomdata, loc); + //This might be equivalent and more modular: position(i, j, k, geom, loc); loc[0] = problo[0] + (0.5_rt + i) * dx[0]; #if AMREX_SPACEDIM >= 2 @@ -759,13 +759,13 @@ extern "C" } void ca_derangmomy (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { int idir = 1; - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); auto const dat = datfab.array(); auto const L = derfab.array(); @@ -807,13 +807,13 @@ extern "C" } void ca_derangmomz (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { int idir = 2; - auto dx = geomdata.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto dx = geom.CellSizeArray(); + auto problo = geom.ProbLoArray(); auto const dat = datfab.array(); auto const L = derfab.array(); @@ -856,7 +856,7 @@ extern "C" } void ca_derkineng (const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -887,7 +887,7 @@ extern "C" } void ca_derspec(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -903,7 +903,7 @@ extern "C" void ca_derye(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -925,7 +925,7 @@ extern "C" } void ca_derabar(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -947,19 +947,19 @@ extern "C" } void ca_dermagvort(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); - const int coord_type = geomdata.Coord(); + const int coord_type = geom.Coord(); #if AMREX_SPACEDIM == 2 - auto problo = geomdata.ProbLoArray(); + auto problo = geom.ProbLoArray(); #endif amrex::ParallelFor(bx, @@ -1046,18 +1046,18 @@ extern "C" } void ca_derdivu(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); - auto problo = geomdata.ProbLoArray(); + auto problo = geom.ProbLoArray(); - const int coord_type = geomdata.Coord(); + const int coord_type = geom.Coord(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -1112,7 +1112,7 @@ extern "C" } void ca_derstate(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1138,7 +1138,7 @@ extern "C" #ifdef MHD void ca_dermagcenx(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1154,7 +1154,7 @@ extern "C" } void ca_dermagceny(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1170,7 +1170,7 @@ extern "C" } void ca_dermagcenz(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& /*geomdata*/, + const FArrayBox& datfab, const Geometry& /*geom*/, Real /*time*/, const int* /*bcrec*/, int /*level*/) { @@ -1186,14 +1186,14 @@ extern "C" } void ca_derdivb(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/, - const FArrayBox& datfab, const Geometry& geomdata, + const FArrayBox& datfab, const Geometry& geom, Real /*time*/, const int* /*bcrec*/, int /*level*/) { auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geomdata.CellSizeArray(); + auto dx = geom.CellSizeArray(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept From b9c7c1f96f860b31d6135115f302e45af3f47e8b Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 16:41:53 -0500 Subject: [PATCH 04/10] fix compiler warning and compilation error --- Source/driver/Derive.cpp | 10 ---------- Source/gravity/Gravity.cpp | 2 ++ 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index fcca2c219e..f1f5f70f2c 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -577,8 +577,6 @@ extern "C" auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geom.CellSizeArray(); - auto problo = geom.ProbLoArray(); auto geomdata = geom.data(); amrex::ParallelFor(bx, @@ -623,8 +621,6 @@ extern "C" auto const dat = datfab.array(); auto const der = derfab.array(); - auto dx = geom.CellSizeArray(); - auto problo = geom.ProbLoArray(); auto geomdata = geom.data(); amrex::ParallelFor(bx, @@ -694,8 +690,6 @@ extern "C" { int idir = 0; - auto dx = geom.CellSizeArray(); - auto problo = geom.ProbLoArray(); auto geomdata = geom.data(); auto const dat = datfab.array(); @@ -728,8 +722,6 @@ extern "C" { int idir = 1; - auto dx = geom.CellSizeArray(); - auto problo = geom.ProbLoArray(); auto geomdata = geom.data(); auto const dat = datfab.array(); @@ -761,8 +753,6 @@ extern "C" { int idir = 2; - auto dx = geom.CellSizeArray(); - auto problo = geom.ProbLoArray(); auto geomdata = geom.data(); auto const dat = datfab.array(); diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 82bc7e699e..34b9a74354 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -1334,6 +1334,7 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab& const Real dr = dx[0] / static_cast(gravity::drdxfac); const auto problo = geom.ProbLoArray(); + const auto geomdata = geom.data(); #ifdef _OPENMP #pragma omp parallel @@ -1443,6 +1444,7 @@ Gravity::compute_radial_mass(const Box& bx, Real drinv = 1.0_rt / dr; const int coord_type = geom.Coord(); + const auto geomdata = geom.data(); AMREX_ALWAYS_ASSERT(coord_type >= 0 && coord_type <= 2); From bb9877c81e6540aa66b908cdb675e09fc25f54d3 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 16:47:58 -0500 Subject: [PATCH 05/10] fix typo --- Source/reactions/Castro_react.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index c04891d531..7d005e726c 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -210,7 +210,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra #ifdef MODEL_PARSER const auto problo = geom.ProbLoArray(); #endif - const auto geomdata = goem.data(); + const auto geomdata = geom.data(); #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) From 6f591529fdb94b6c9d45c29f8ac208abc2b0cb2c Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 16:54:58 -0500 Subject: [PATCH 06/10] more fix compilation error --- Source/reactions/Castro_react.cpp | 3 +-- Source/scf/scf_relax.cpp | 9 +-------- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index 7d005e726c..8d39095ff3 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -208,9 +208,8 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra const auto dx = geom.CellSizeArray(); #ifdef MODEL_PARSER - const auto problo = geom.ProbLoArray(); -#endif const auto geomdata = geom.data(); +#endif #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp index 7b1899ee13..674a509f89 100644 --- a/Source/scf/scf_relax.cpp +++ b/Source/scf/scf_relax.cpp @@ -243,11 +243,10 @@ Castro::do_hscf_solve() [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { const auto *dx = geomdata.CellSize(); - const auto *problo = geomdata.ProbLo(); // The below assumes we are rotating on the z-axis. - GpuArrayr; + GpuArrayr = {0.0}; position(i, j, k, geomdata, r); // Do a trilinear interpolation to find the contribution from @@ -369,7 +368,6 @@ Castro::do_hscf_solve() [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple { const auto *dx = geomdata.CellSize(); - const auto *problo = geomdata.ProbLo(); // The below assumes we are rotating on the z-axis. @@ -441,9 +439,6 @@ Castro::do_hscf_solve() // enthalpy + gravitational potential + rotational potential = const // We already have the constant, so our goal is to construct the enthalpy field. - const auto *dx = geomdata.CellSize(); - const auto *problo = geomdata.ProbLo(); - GpuArray r = {0.0}; position(i, j, k, geomdata, r); @@ -618,8 +613,6 @@ Castro::do_hscf_solve() { Real dM = 0.0, dK = 0.0, dU = 0.0, dE = 0.0; - const auto* problo = geomdata.ProbLo(); - GpuArray r = {0.0}; position(i, j, k, geomdata, r); From 6d10cc457292f39001dfc1d796e820f714cd28ab Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 17:02:04 -0500 Subject: [PATCH 07/10] more fixes --- Source/rotation/Rotation.cpp | 4 +--- Source/sources/Castro_sponge.cpp | 4 +--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/Source/rotation/Rotation.cpp b/Source/rotation/Rotation.cpp index 082d7beff8..610b7a0e46 100644 --- a/Source/rotation/Rotation.cpp +++ b/Source/rotation/Rotation.cpp @@ -21,9 +21,7 @@ Castro::fill_rotational_psi(const Box& bx, auto omega = get_omega(); Real denom = omega[0] * omega[0] + omega[1] * omega[1] + omega[2] * omega[2]; - auto problo = geom.ProbLoArray(); - - auto dx = geom.CellSizeArray(); + auto geomdata = geom.data(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept diff --git a/Source/sources/Castro_sponge.cpp b/Source/sources/Castro_sponge.cpp index ebf0fef625..67be2b39a5 100644 --- a/Source/sources/Castro_sponge.cpp +++ b/Source/sources/Castro_sponge.cpp @@ -79,9 +79,7 @@ Castro::apply_sponge(const Box& bx, alpha = 0.0_rt; } - auto dx = geom.CellSizeArray(); - auto problo = geom.ProbLoArray(); - const auto geomdata = geom.data(); + auto geomdata = geom.data(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept From 78c1d0224c8d1a3a845505a3848811c286176adc Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 17:09:52 -0500 Subject: [PATCH 08/10] fix --- Source/reactions/Castro_react.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index 8d39095ff3..c763801402 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -557,7 +557,7 @@ Castro::react_state(Real time, Real dt) int lsdc_iteration = sdc_iteration; const auto dx = geom.CellSizeArray(); - const auto problo = geom.ProbLoArray(); + const geomdata = geom.data(); #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) From a6974022cadfa2a4cf936aac7c5480e77eaa4367 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 17:14:37 -0500 Subject: [PATCH 09/10] fix --- Source/reactions/Castro_react.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index c763801402..69625998ec 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -557,7 +557,7 @@ Castro::react_state(Real time, Real dt) int lsdc_iteration = sdc_iteration; const auto dx = geom.CellSizeArray(); - const geomdata = geom.data(); + const auto geomdata = geom.data(); #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) From 00c76bd0e72960829cedba04c2854b02afb5f551 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 7 Nov 2024 17:21:41 -0500 Subject: [PATCH 10/10] fix compile warning --- Source/gravity/Gravity.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 34b9a74354..a9cc9dbab1 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -1333,7 +1333,6 @@ Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab& const auto dx = geom.CellSizeArray(); const Real dr = dx[0] / static_cast(gravity::drdxfac); - const auto problo = geom.ProbLoArray(); const auto geomdata = geom.data(); #ifdef _OPENMP