From 39a11ffed98739bde382de9c2c8bb0ada820a09e Mon Sep 17 00:00:00 2001 From: Zhi Date: Sun, 8 Sep 2024 17:52:01 -0400 Subject: [PATCH 1/5] update local area and volume calculations --- Source/driver/Castro_util.H | 47 +++++++++++++++++++++++++++++++------ 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index 4202132493..43871fc78e 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,8 +318,27 @@ 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 = M_PI * std::sin(theta) * dx[0] * (2.0_rt * r + dx[0]); + } + } + #else // Cartesian From 23e64a7d4f2e9f2d38488e6c02fc0768f1f39567 Mon Sep 17 00:00:00 2001 From: Zhi Date: Sun, 8 Sep 2024 18:02:02 -0400 Subject: [PATCH 2/5] fix theta area --- Source/driver/Castro_util.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index 43871fc78e..3c4c40aed7 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -333,7 +333,7 @@ Real area(const int& i, const int& j, const int& k, Real r = geomdata.ProbLo()[0] + (static_cast(i) + 0.5_rt) * dx[0]; Real theta = geomdata.ProbLo()[1] + static_cast(j) * dx[1]; - a = M_PI * std::sin(theta) * dx[0] * (2.0_rt * r + dx[0]); + a = 2.0_rt * M_PI * std::sin(theta) * dx[0] * r; } } From 59ca30f68e0de3407c2345cbcef8dfd173961984 Mon Sep 17 00:00:00 2001 From: Zhi Date: Sun, 8 Sep 2024 18:06:43 -0400 Subject: [PATCH 3/5] remove empty line, cleanup --- Source/driver/Castro_util.H | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index 3c4c40aed7..e7b8751d12 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -333,12 +333,11 @@ Real area(const int& i, const int& j, const int& k, 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) * dx[0] * r; + a = 2.0_rt * M_PI * std::sin(theta) * r * dx[0]; } } - #else // Cartesian From 9423c9a90ff586b66ecce10f704f593093dbd204 Mon Sep 17 00:00:00 2001 From: Zhi Date: Sun, 8 Sep 2024 18:39:50 -0400 Subject: [PATCH 4/5] update divu calculation for spherical 2d geometry --- Source/hydro/advection_util.cpp | 95 ++++++++++++++++++++++----------- 1 file changed, 64 insertions(+), 31 deletions(-) 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; From fea3b4dfdd0a3d424c85b8ccbe8856092bb1a78d Mon Sep 17 00:00:00 2001 From: Zhi Date: Sun, 8 Sep 2024 18:51:30 -0400 Subject: [PATCH 5/5] revert Castro_util.H to development --- Source/driver/Castro_util.H | 46 ++++++------------------------------- 1 file changed, 7 insertions(+), 39 deletions(-) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index e7b8751d12..4202132493 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -157,12 +157,10 @@ 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. - // this function only provides Cartesian in 1D/2D/3D, - // Cylindrical (R-Z) in 2D, - // and Spherical in 1D and 2D (R-THETA). - // + // 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. amrex::ignore_unused(i); amrex::ignore_unused(j); @@ -212,7 +210,8 @@ Real volume(const int& i, const int& j, const int& k, vol = dx[0] * dx[1]; - } else if (coord == 1) { + } + else { // Cylindrical @@ -221,20 +220,6 @@ 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 @@ -305,7 +290,8 @@ Real area(const int& i, const int& j, const int& k, a = dx[0]; } - } else if (coord == 1) { + } + else { // Cylindrical @@ -318,24 +304,6 @@ 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