address rotation source in spherical 2d coordinate (#2967)

Docs/source/rotation.rst

@@ -40,6 +40,10 @@ The "vertical direction" is defined as follows:
 
   * ``coord_sys = 1``, (r,z): along the z-axis
 
+  * ``coord_sys = 2``, (r, :math:`\theta`): along the z-axis after converting
+    to the spherical coordinate, i.e.
+    :math:`\hat{z} = \cos{\theta} \hat{r} = \sin{\theta} \hat{\theta}`.
+
 * 3D
 
   * ``coord_sys = 0``, (x,y,z): along the z-axis
@@ -76,7 +80,9 @@ The main parameters that affect rotation are:
    solve for the update implicitly (default: 1)
 
 -  ``castro.rot_axis`` : rotation axis (default: 3
-   (Cartesian); 2 (cylindrical))
+   (Cartesian); 2 (cylindrical)). This parameter doesn't affect
+   spherical coordinate since rotation axis is automatically set
+   to z-axis.
 
 For completeness, we show below a derivation of the source terms that
 appear in the momentum and total energy evolution equations upon

Exec/hydro_tests/rotating_torus/problem_initialize.H

@@ -19,13 +19,12 @@ void problem_initialize ()
     auto omega = get_omega();
 #else
     // Provide a dummy value so that we can compile without rotation.
-    Real omega[3] = {0.0_rt, 0.0_rt, 2.0_rt * M_PI};
+    Real omega = 2.0_rt * M_PI;
 #endif
 
     // Figure out R_0, the maximum pressure radius.
 
-    problem::density_maximum_radius = std::pow(C::Gconst * point_mass /
-                                               (omega[0] * omega[0] + omega[1] * omega[1] + omega[2] * omega[2]), 1.0_rt/3.0_rt);
+    problem::density_maximum_radius = std::pow(C::Gconst * point_mass / (omega * omega), 1.0_rt/3.0_rt);
 
     // Maximum and minimum vertical extent of the torus is the same as the radial extent
 

Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H

@@ -17,7 +17,7 @@ void problem_initialize_state_data (int i, int j, int k,
     const Real* probhi = geomdata.ProbHi();
 
 #ifdef ROTATION
-    auto omega = get_omega();
+    auto omega = get_omega_vec(geomdata, j);
 #else
     // Provide a dummy value so that we can compile without rotation.
     Real omega[3] = {0.0_rt, 0.0_rt, 2.0_rt * M_PI};

Exec/science/wdmerger/Prob.cpp

@@ -798,6 +798,7 @@ Castro::update_relaxation(Real time, Real dt) {
     get_lagrange_points(mass_P, mass_S, com_P, com_S, L1, L2, L3);
 
     const auto dx = geom.CellSizeArray();
+    const auto coord = geom.Coord();
     GeometryData geomdata = geom.data();
 
     MultiFab& phi_new = get_new_data(PhiGrav_Type);
@@ -833,7 +834,8 @@ Castro::update_relaxation(Real time, Real dt) {
 
                 // Compute the effective potential.
 
-                Real phiEff = phi(i,j,k) + rotational_potential(r);
+                auto omega_vec = get_omega_vec(geomdata, j);
+                Real phiEff = phi(i,j,k) + rotational_potential(r, omega_vec, coord);
 
                 for (int n = 0; n < 3; ++n) {
                     r[n] -= L1[n];
@@ -901,7 +903,8 @@ Castro::update_relaxation(Real time, Real dt) {
                 GpuArray<Real, 3> r;
                 position(i, j, k, geomdata, r);
 
-                Real phiEff = phi(i,j,k) + rotational_potential(r);
+                auto omega_vec = get_omega_vec(geomdata, j);
+                Real phiEff = phi(i,j,k) + rotational_potential(r, omega_vec, coord);
 
                 Real done = 0.0_rt;
 

Exec/science/wdmerger/Problem_Derive.cpp

@@ -419,6 +419,8 @@ void ca_derphieff(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -438,8 +440,9 @@ void ca_derphieff(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncomp*/
 #else
         loc[2] = 0.0_rt;
 #endif
+        auto omega = get_omega_vec(geomdata, j);
 
-        der(i,j,k,0) = dat(i,j,k,0) + rotational_potential(loc);
+        der(i,j,k,0) = dat(i,j,k,0) + rotational_potential(loc, omega, coord);
     });
 }
 
@@ -458,6 +461,8 @@ void ca_derphieffpm_p(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -492,7 +497,9 @@ void ca_derphieffpm_p(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
             loc[iloc] -= problem::center[iloc];
         }
 
-        der(i,j,k,0) = -C::Gconst * problem::mass_P / r + rotational_potential(loc);
+        auto omega = get_omega_vec(geomdata, j);
+
+        der(i,j,k,0) = -C::Gconst * problem::mass_P / r + rotational_potential(loc, omega, coord);
     });
 }
 
@@ -508,6 +515,8 @@ void ca_derphieffpm_s(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -542,7 +551,9 @@ void ca_derphieffpm_s(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*nco
             loc[iloc] -= problem::center[iloc];
         }
 
-        der(i,j,k,0) = -C::Gconst * problem::mass_S / r + rotational_potential(loc);
+        auto omega = get_omega_vec(geomdata, j);
+
+        der(i,j,k,0) = -C::Gconst * problem::mass_S / r + rotational_potential(loc, omega, coord);
     });
 }
 
@@ -572,6 +583,8 @@ void ca_derrhophiRot(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncom
 
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
+    const auto coord = geom.Coord();
+    const auto geomdata = geom.data();
 
     amrex::ParallelFor(bx,
     [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
@@ -592,7 +605,9 @@ void ca_derrhophiRot(const Box& bx, FArrayBox& derfab, int /*dcomp*/, int /*ncom
         loc[2] = 0.0_rt;
 #endif
 
-        der(i,j,k,0) = dat(i,j,k,0) * rotational_potential(loc);
+        auto omega = get_omega_vec(geomdata, j);
+
+        der(i,j,k,0) = dat(i,j,k,0) * rotational_potential(loc, omega, coord);
     });
 }
 

Exec/science/wdmerger/problem_initialize_state_data.H

@@ -16,7 +16,7 @@ void problem_initialize_state_data (int i, int j, int k,
     // or secondary (in which case interpolate from the respective model)
     // or if we are in an ambient zone.
 
-    auto omega = get_omega();
+    auto omega = get_omega_vec(geomdata, j);
 
     GpuArray<Real, 3> loc;
     position(i, j, k, geomdata, loc);

Exec/science/wdmerger/wdmerger_util.H

@@ -89,7 +89,9 @@ Real stellar_mask (int i, int j, int k,
         loc[iloc] -= problem::center[iloc];
     }
 
-    Real phi_rot = rotational_potential(loc);
+    auto omega = get_omega_vec(geomdata, j);
+
+    Real phi_rot = rotational_potential(loc, omega, geomdata.Coord());
 
     Real phi_p = -C::Gconst * problem::mass_P / r_P + phi_rot;
     Real phi_s = -C::Gconst * problem::mass_S / r_S + phi_rot;

Source/driver/_cpp_parameters

@@ -490,7 +490,12 @@ rot_source_type              int           4                  ROTATION
 # for better coupling of the Coriolis terms
 implicit_rotation_update     bool           1                  ROTATION
 
-# the coordinate axis (:math:`x=1`, :math:`y=2`, :math:`z=3`) for the rotation vector
+# the coordinate axis for the rotation vector
+# For Cartesian: (:math:`x=1`, :math:`y=2`, :math:`z=3`)
+# For non-Cartesian coordinates, this parameter doesn't do anything because:
+# For RZ (Cylindrical 2D), it is automatically set to z-axis (rot_axis = 2)
+# For Spherical2D, it is also assumed to be in the z-axis
+# i.e. cos(theta) r_hat - sin(theta) theta_hat in Spherical coordinate.
 rot_axis                     int           3                  ROTATION
 
 # include a central point mass

Source/driver/sum_utils.cpp

@@ -400,7 +400,8 @@ Castro::gwstrain (Real time,
 
             GpuArray<Real, 3> pos{r};
 #ifdef ROTATION
-            pos = inertial_rotation(r, time);
+            auto omega = get_omega_vec(geomdata, j);
+            pos = inertial_rotation(r, omega, time);
 #endif
 
             // For constructing the velocity in the inertial frame, we need to
@@ -431,7 +432,7 @@ Castro::gwstrain (Real time,
 
             GpuArray<Real, 3> inertial_g{g};
 #ifdef ROTATION
-            inertial_g = inertial_rotation(g, time);
+            inertial_g = inertial_rotation(g, omega, time);
 #endif
 
             // Absorb the factor of 2 outside the integral into the zone mass, for efficiency.