Skip to content

Commit

Permalink
update local area and volume calculations for 2D spherical geometry (#…
Browse files Browse the repository at this point in the history
…2953)

This adds local area and volume calculations for 2D spherical geometry. We assume dx[0] = dr and dx[1] = dtheta.
  • Loading branch information
zhichen3 authored Sep 10, 2024
1 parent 0070766 commit a12d278
Showing 1 changed file with 39 additions and 7 deletions.
46 changes: 39 additions & 7 deletions Source/driver/Castro_util.H
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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

Expand All @@ -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<Real>(i) * dx[0];
Real r_r = geomdata.ProbLo()[0] + static_cast<Real>(i+1) * dx[0];
Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(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
Expand Down Expand Up @@ -290,8 +305,7 @@ Real area(const int& i, const int& j, const int& k,
a = dx[0];
}

}
else {
} else if (coord == 1) {

// Cylindrical

Expand All @@ -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<Real>(i) * dx[0];
Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(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<Real>(i) + 0.5_rt) * dx[0];
Real theta = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];

a = 2.0_rt * M_PI * std::sin(theta) * r * dx[0];
}

}

#else
Expand Down

0 comments on commit a12d278

Please sign in to comment.