Skip to content

Commit

Permalink
subtract off center in position and create separate distance function
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhichen3
zhichen3 committed Nov 7, 2024
1 parent ee78fcd commit dbd972b
Show file tree
Hide file tree
Showing 17 changed files with 86 additions and 270 deletions.
19 changes: 6 additions & 13 deletions Source/driver/Castro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real, 3> loc;

loc[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0];
#if AMREX_SPACEDIM >= 2
loc[1] = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1];
#endif
#if AMREX_SPACEDIM == 3
loc[2] = problo[2] + (static_cast<Real>(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;
Expand Down Expand Up @@ -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;
}

}
Expand Down
25 changes: 23 additions & 2 deletions Source/driver/Castro_util.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real, 3>& loc,
Expand Down Expand Up @@ -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) {
Expand All @@ -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<Real, 3>& 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
{

Expand Down
116 changes: 24 additions & 92 deletions Source/driver/Derive.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -584,17 +584,9 @@ extern "C"
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{

Real x = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
#if AMREX_SPACEDIM >= 2
Real y = problo[1] + (static_cast<Real>(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<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
#else
Real z = 0.0_rt;
#endif
GpuArray<Real, 3> loc;

position(i, j, k, geomdata, loc);

if (domain_is_plane_parallel) {
#if AMREX_SPACEDIM == 2
Expand All @@ -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 );
}

});
Expand All @@ -639,17 +631,9 @@ extern "C"
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{

Real x = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
#if AMREX_SPACEDIM >= 2
Real y = problo[1] + (static_cast<Real>(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<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
#else
Real z = 0.0_rt;
#endif
GpuArray<Real, 3> loc;

position(i, j, k, geomdata, loc);

if (domain_is_plane_parallel) {
#if AMREX_SPACEDIM == 2
Expand All @@ -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
Expand All @@ -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));
}
Expand Down Expand Up @@ -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<Real, 3> 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)
Expand Down Expand Up @@ -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<Real, 3> 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);
Expand Down Expand Up @@ -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<Real, 3> 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);
Expand Down
4 changes: 0 additions & 4 deletions Source/driver/sum_utils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -384,10 +384,6 @@ Castro::gwstrain (Real time,
GpuArray<Real, 3> 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);
Expand Down
8 changes: 1 addition & 7 deletions Source/gravity/Castro_gravity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -375,10 +375,7 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in,

#ifdef HYBRID_MOMENTUM
GpuArray<Real, 3> 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<Real, 3> hybrid_src;

Expand Down Expand Up @@ -574,9 +571,6 @@ void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old,
#ifdef HYBRID_MOMENTUM
GpuArray<Real, 3> loc;
position(i, j, k, geomdata, loc);
for (int n = 0; n < 3; ++n) {
loc[n] -= problem::center[n];
}

GpuArray<Real, 3> hybrid_src;

Expand Down
28 changes: 6 additions & 22 deletions Source/gravity/Gravity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real, 3> loc;

loc[0] = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];

#if AMREX_SPACEDIM >= 2
loc[1] = problo[1] + (static_cast<Real>(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<Real>(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<int>(r / dr);

Expand Down Expand Up @@ -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<Real>(i) + 0.5_rt) * dx[0] - problem::center[0];
Real lo_i = problo[0] + static_cast<Real>(i) * dx[0] - problem::center[0];

Real yc = problo[1] + (static_cast<Real>(j) + 0.5_rt) * dx[1] - problem::center[1];
Real lo_j = problo[1] + static_cast<Real>(j) * dx[1] - problem::center[1];

Real zc = problo[2] + (static_cast<Real>(k) + 0.5_rt) * dx[2] - problem::center[2];
Real lo_k = problo[2] + static_cast<Real>(k) * dx[2] - problem::center[2];

Real r = std::sqrt(xc * xc + yc * yc + zc * zc);
GpuArray<Real ,3> loc;
position(i, j, k, geomdata, loc);
Real r = distance(geomdata, loc);

int index = static_cast<int>(r * drinv);

// We may be coming in here with a masked out zone (in a zone on a coarse
Expand Down
3 changes: 0 additions & 3 deletions Source/hydro/Castro_ctu_hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real>::min());
Real RInv = 1.0_rt / R;
Expand Down
Loading

0 comments on commit dbd972b

Please sign in to comment.