From 34d2b9a4b815f73da08619562f2ddc79e15ae2f1 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Fri, 20 Dec 2024 15:30:55 -0500 Subject: [PATCH] address rotation source in spherical 2d coordinate (#2967) --- Docs/source/rotation.rst | 8 +- .../rotating_torus/problem_initialize.H | 5 +- .../problem_initialize_state_data.H | 2 +- Exec/science/wdmerger/Prob.cpp | 7 +- Exec/science/wdmerger/Problem_Derive.cpp | 23 +++- .../wdmerger/problem_initialize_state_data.H | 2 +- Exec/science/wdmerger/wdmerger_util.H | 4 +- Source/driver/_cpp_parameters | 7 +- Source/driver/sum_utils.cpp | 5 +- Source/rotation/Rotation.H | 130 ++++++++++++++++-- Source/rotation/Rotation.cpp | 14 +- Source/rotation/rotation_sources.cpp | 69 ++-------- Source/scf/scf_relax.cpp | 12 +- 13 files changed, 194 insertions(+), 94 deletions(-) diff --git a/Docs/source/rotation.rst b/Docs/source/rotation.rst index 2bb60ea4a4..485f6983d9 100644 --- a/Docs/source/rotation.rst +++ b/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 diff --git a/Exec/hydro_tests/rotating_torus/problem_initialize.H b/Exec/hydro_tests/rotating_torus/problem_initialize.H index 479408fb46..0e4c3c4650 100644 --- a/Exec/hydro_tests/rotating_torus/problem_initialize.H +++ b/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 diff --git a/Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H b/Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H index 5ec65d6b8a..53a4630134 100644 --- a/Exec/hydro_tests/rotating_torus/problem_initialize_state_data.H +++ b/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}; diff --git a/Exec/science/wdmerger/Prob.cpp b/Exec/science/wdmerger/Prob.cpp index ba28a87a18..1e5b246f0d 100644 --- a/Exec/science/wdmerger/Prob.cpp +++ b/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 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; diff --git a/Exec/science/wdmerger/Problem_Derive.cpp b/Exec/science/wdmerger/Problem_Derive.cpp index 019b56a728..776b7078c5 100644 --- a/Exec/science/wdmerger/Problem_Derive.cpp +++ b/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); }); } diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 3b19019039..f9f52414f8 100644 --- a/Exec/science/wdmerger/problem_initialize_state_data.H +++ b/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 loc; position(i, j, k, geomdata, loc); diff --git a/Exec/science/wdmerger/wdmerger_util.H b/Exec/science/wdmerger/wdmerger_util.H index d2c28b3bfa..9a8c4fe327 100644 --- a/Exec/science/wdmerger/wdmerger_util.H +++ b/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; diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 4983d32b86..3ba6449ca9 100644 --- a/Source/driver/_cpp_parameters +++ b/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 diff --git a/Source/driver/sum_utils.cpp b/Source/driver/sum_utils.cpp index 6423b50b70..f19a02f2d6 100644 --- a/Source/driver/sum_utils.cpp +++ b/Source/driver/sum_utils.cpp @@ -400,7 +400,8 @@ Castro::gwstrain (Real time, GpuArray 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 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. diff --git a/Source/rotation/Rotation.H b/Source/rotation/Rotation.H index 878ae4a6b0..06f2a04490 100644 --- a/Source/rotation/Rotation.H +++ b/Source/rotation/Rotation.H @@ -6,36 +6,68 @@ #include #include + /// -/// Return the omega vector corresponding to the current rotational period +/// Return the magnitude of omega along the rotating axis. /// AMREX_GPU_HOST_DEVICE AMREX_INLINE -GpuArray get_omega() +Real get_omega() { - GpuArray omega = {0.0_rt}; + Real omega = 0.0_rt; // If rotational_period is less than zero, that means rotation is disabled, and so we should effectively // shut off the source term by setting omega = 0. if (castro::rotational_period > 0.0_rt) { - omega[castro::rot_axis - 1] = 2.0_rt * M_PI / castro::rotational_period; + omega = 2.0_rt * M_PI / castro::rotational_period; } return omega; } +/// +/// Return the omega vector corresponding to the current rotational period +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +GpuArray get_omega_vec(const GeometryData& geomdata, + const int j) +{ + GpuArray omega_vec = {0.0_rt}; + + const auto coord = geomdata.Coord(); + const auto problo = geomdata.ProbLo(1); + const auto dtheta = geomdata.CellSize(1); + + auto omega = get_omega(); + + // If coord == 2, i.e. Spherical2D, then it is assumed that rotational axis is in z-axis + // then we convert back to r and theta direction. + if (coord == 2) { + Real theta = problo + (static_cast(j) + 0.5_rt) * dtheta; + omega_vec[0] = std::cos(theta) * omega; + omega_vec[1] = -std::sin(theta) * omega; + } else { + omega_vec[castro::rot_axis - 1] = omega; + } + + return omega_vec; +} + /// /// Compute the rotational acceleration for a single zone /// (Coriolis and centrifugal) /// /// @param r distance from origin of rotation vector /// @param v velocity +/// @param omega angular velocity vector array +/// @param coord coordinate system /// @param coriolis do we include the Coriolis force /// @param Sr rotational acceleration /// AMREX_GPU_HOST_DEVICE AMREX_INLINE void -rotational_acceleration(GpuArray& r, GpuArray& v, +rotational_acceleration(const GpuArray& r, const GpuArray& v, + const GpuArray& omega, const int coord, const bool coriolis, Real* Sr) { // Given a position and velocity, calculate @@ -48,8 +80,6 @@ rotational_acceleration(GpuArray& r, GpuArray& v, Sr[1] = 0.0; Sr[2] = 0.0; - auto omega = get_omega(); - // Allow the various terms to be turned off. This is often used // for diagnostic purposes, but there are genuine science cases // for disabling certain terms in some cases (in particular, when @@ -63,12 +93,21 @@ rotational_acceleration(GpuArray& r, GpuArray& v, bool c2 = (castro::rotation_include_coriolis == 1 && coriolis) ? true : false; + auto r_vec = r; + GpuArray omega_cross_v; cross_product(omega, v, omega_cross_v); if (c1) { + // For Spherical coordinate, the physical position vector is just r rhat, + // but the input r will be in the format of (r,theta,0). So change it to (r,0,0) + + if (coord == 2) { + r_vec[1] = 0.0_rt; + r_vec[2] = 0.0_rt; + } GpuArray omega_cross_r; - cross_product(omega, r, omega_cross_r); + cross_product(omega, r_vec, omega_cross_r); GpuArray omega_cross_omega_cross_r; cross_product(omega, omega_cross_r, omega_cross_omega_cross_r); @@ -88,14 +127,19 @@ rotational_acceleration(GpuArray& r, GpuArray& v, AMREX_GPU_HOST_DEVICE AMREX_INLINE Real -rotational_potential(GpuArray& r) { +rotational_potential(const GpuArray& r, const GpuArray& omega, + const int coord) { // Construct rotational potential, phi_R = -1/2 | omega x r |**2 // Real phi = 0.0_rt; + auto r_vec = r; - auto omega = get_omega(); + if (coord == 2) { + r_vec[1] = 0.0_rt; + r_vec[2] = 0.0_rt; + } if (rotation_include_centrifugal == 1) { @@ -133,7 +177,12 @@ inertial_to_rotational_velocity (const int i, const int j, const int k, loc[dir] -= problem::center[dir]; } - auto omega = get_omega(); + if (geomdata.Coord() == 2) { + loc[1] = 0.0_rt; + loc[2] = 0.0_rt; + } + + auto omega = get_omega_vec(geomdata, j); // do the cross product Omega x loc v[0] += -(omega[1]*loc[2] - omega[2]*loc[1]); @@ -163,7 +212,12 @@ rotational_to_inertial_velocity (const int i, const int j, const int k, loc[dir] -= problem::center[dir]; } - auto omega = get_omega(); + if (geomdata.Coord() == 2) { + loc[1] = 0.0_rt; + loc[2] = 0.0_rt; + } + + auto omega = get_omega_vec(geomdata, j); // do the cross product Omega x loc v[0] += (omega[1]*loc[2] - omega[2]*loc[1]); @@ -177,15 +231,14 @@ rotational_to_inertial_velocity (const int i, const int j, const int k, // since the beginning of the simulation. AMREX_GPU_HOST_DEVICE AMREX_INLINE -GpuArray inertial_rotation(const GpuArray& vec, Real time) +GpuArray inertial_rotation(const GpuArray& vec, + const GpuArray& omega, Real time) { GpuArray vec_i{}; Array1D theta{}; Array2D rot_matrix{}; - auto omega = get_omega(); - if (castro::do_rotation == 1) { for (int n = 0; n < 3; ++n) { theta(n) = omega[n] * time; @@ -219,4 +272,51 @@ GpuArray inertial_rotation(const GpuArray& vec, Real time) return vec_i; } + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void fill_dt_omega_matrix(const Real dt, const GpuArray& omega, + Array2D& dt_omega_matrix) +{ + // Fill in dt_omega_matrix. This is used in function corrrsrc + + Real dt_omega[3]; + + // Don't do anything here if we've got the Coriolis force disabled. + + if (castro::rotation_include_coriolis == 1) { + + for (int idir = 0; idir < 3; idir++) { + dt_omega[idir] = dt * omega[idir]; + } + + } else { + + for (auto& e : dt_omega) { + e = 0.0_rt; + } + + } + + dt_omega_matrix(0, 0) = 1.0_rt + dt_omega[0] * dt_omega[0]; + dt_omega_matrix(0, 1) = dt_omega[0] * dt_omega[1] + dt_omega[2]; + dt_omega_matrix(0, 2) = dt_omega[0] * dt_omega[2] - dt_omega[1]; + + dt_omega_matrix(1, 0) = dt_omega[1] * dt_omega[0] - dt_omega[2]; + dt_omega_matrix(1, 1) = 1.0_rt + dt_omega[1] * dt_omega[1]; + dt_omega_matrix(1, 2) = dt_omega[1] * dt_omega[2] + dt_omega[0]; + + dt_omega_matrix(2, 0) = dt_omega[2] * dt_omega[0] + dt_omega[1]; + dt_omega_matrix(2, 1) = dt_omega[2] * dt_omega[1] - dt_omega[0]; + dt_omega_matrix(2, 2) = 1.0_rt + dt_omega[2] * dt_omega[2]; + + for (int l = 0; l < 3; l++) { + for (int m = 0; m < 3; m++) { + dt_omega_matrix(l, m) /= (1.0_rt + dt_omega[0] * dt_omega[0] + + dt_omega[1] * dt_omega[1] + + dt_omega[2] * dt_omega[2]); + } + } + +} + #endif diff --git a/Source/rotation/Rotation.cpp b/Source/rotation/Rotation.cpp index 9a5f858fec..5de2aecd4d 100644 --- a/Source/rotation/Rotation.cpp +++ b/Source/rotation/Rotation.cpp @@ -18,12 +18,13 @@ Castro::fill_rotational_psi(const Box& bx, // routine uniquely determines the rotation law. For the other // rotation laws, we would simply divide by v_0^2 or j_0^2 instead. - auto omega = get_omega(); - Real denom = omega[0] * omega[0] + omega[1] * omega[1] + omega[2] * omega[2]; + const auto omega = get_omega(); + Real denom = omega * omega; - auto problo = geom.ProbLoArray(); - - auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const auto dx = geom.CellSizeArray(); + const auto coord = geom.Coord(); + const auto geomdata = geom.data(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -44,7 +45,8 @@ Castro::fill_rotational_psi(const Box& bx, #endif if (denom != 0.0) { - psi(i,j,k) = rotational_potential(r) / denom; + auto omega_vec = get_omega_vec(geomdata, j); + psi(i,j,k) = rotational_potential(r, omega_vec, coord) / denom; } else { psi(i,j,k) = 0.0; diff --git a/Source/rotation/rotation_sources.cpp b/Source/rotation/rotation_sources.cpp index 4c034810d8..d116da5321 100644 --- a/Source/rotation/rotation_sources.cpp +++ b/Source/rotation/rotation_sources.cpp @@ -12,6 +12,7 @@ Castro::rsrc(const Box& bx, const Real dt) { GeometryData geomdata = geom.data(); + const auto coord = geomdata.Coord(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -31,6 +32,8 @@ Castro::rsrc(const Box& bx, loc[dir] -= problem::center[dir]; } + auto omega = get_omega_vec(geomdata, j); + Real rho = uold(i,j,k,URHO); Real rhoInv = 1.0_rt / rho; @@ -47,7 +50,7 @@ Castro::rsrc(const Box& bx, v[2] = uold(i,j,k,UMZ) * rhoInv; bool coriolis = true; - rotational_acceleration(loc, v, coriolis, Sr); + rotational_acceleration(loc, v, omega, coord, coriolis, Sr); for (auto& e : Sr) { e *= rho; @@ -167,6 +170,7 @@ Castro::corrrsrc(const Box& bx, // is the new time at time-level n+1. GeometryData geomdata = geom.data(); + const auto coord = geomdata.Coord(); Real hdtInv = 0.5_rt / dt; @@ -178,53 +182,6 @@ Castro::corrrsrc(const Box& bx, dx[i] = 0.0_rt; } - auto omega = get_omega(); - - Real dt_omega[3]; - - Array2D dt_omega_matrix = {}; - - if (implicit_rotation_update == 1) { - - // Don't do anything here if we've got the Coriolis force disabled. - - if (rotation_include_coriolis == 1) { - - for (int idir = 0; idir < 3; idir++) { - dt_omega[idir] = dt * omega[idir]; - } - - } else { - - for (auto& e : dt_omega) { - e = 0.0_rt; - } - - } - - - dt_omega_matrix(0, 0) = 1.0_rt + dt_omega[0] * dt_omega[0]; - dt_omega_matrix(0, 1) = dt_omega[0] * dt_omega[1] + dt_omega[2]; - dt_omega_matrix(0, 2) = dt_omega[0] * dt_omega[2] - dt_omega[1]; - - dt_omega_matrix(1, 0) = dt_omega[1] * dt_omega[0] - dt_omega[2]; - dt_omega_matrix(1, 1) = 1.0_rt + dt_omega[1] * dt_omega[1]; - dt_omega_matrix(1, 2) = dt_omega[1] * dt_omega[2] + dt_omega[0]; - - dt_omega_matrix(2, 0) = dt_omega[2] * dt_omega[0] + dt_omega[1]; - dt_omega_matrix(2, 1) = dt_omega[2] * dt_omega[1] - dt_omega[0]; - dt_omega_matrix(2, 2) = 1.0_rt + dt_omega[2] * dt_omega[2]; - - for (int l = 0; l < 3; l++) { - for (int m = 0; m < 3; m++) { - dt_omega_matrix(l, m) /= (1.0_rt + dt_omega[0] * dt_omega[0] + - dt_omega[1] * dt_omega[1] + - dt_omega[2] * dt_omega[2]); - } - } - - } - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -245,6 +202,8 @@ Castro::corrrsrc(const Box& bx, loc[dir] -= problem::center[dir]; } + auto omega = get_omega_vec(geomdata, j); + Real rhoo = uold(i,j,k,URHO); Real rhooinv = 1.0_rt / uold(i,j,k,URHO); @@ -267,7 +226,7 @@ Castro::corrrsrc(const Box& bx, vold[2] = uold(i,j,k,UMZ) * rhooinv; bool coriolis = true; - rotational_acceleration(loc, vold, coriolis, Sr_old); + rotational_acceleration(loc, vold, omega, coord, coriolis, Sr_old); for (auto& e : Sr_old) { e *= rhoo; @@ -284,7 +243,7 @@ Castro::corrrsrc(const Box& bx, vnew[1] = unew(i,j,k,UMY) * rhoninv; vnew[2] = unew(i,j,k,UMZ) * rhoninv; - rotational_acceleration(loc, vnew, coriolis, Sr_new); + rotational_acceleration(loc, vnew, omega, coord, coriolis, Sr_new); for (auto& e : Sr_new) { e *= rhon; @@ -309,7 +268,7 @@ Castro::corrrsrc(const Box& bx, Real acc[3]; coriolis = false; - rotational_acceleration(loc, vnew, coriolis, acc); + rotational_acceleration(loc, vnew, omega, coord, coriolis, acc); Real new_mom_tmp[3]; for (int n = 0; n < 3; n++) { @@ -328,6 +287,9 @@ Castro::corrrsrc(const Box& bx, // of the dt_omega_matrix. It also has the correct form if we have disabled // the Coriolis force entirely; at that point it reduces to the identity matrix. + Array2D dt_omega_matrix = {}; + fill_dt_omega_matrix(dt, omega, dt_omega_matrix); + Real new_mom[3] = {}; // new_mom = matmul(dt_omega_matrix, new_mom) @@ -398,7 +360,7 @@ Castro::corrrsrc(const Box& bx, Real acc[3]; coriolis = true; - rotational_acceleration(loc, vnew, coriolis, acc); + rotational_acceleration(loc, vnew, omega, coord, coriolis, acc); Sr_new[0] = rhon * acc[0]; Sr_new[1] = rhon * acc[1]; @@ -456,7 +418,7 @@ Castro::corrrsrc(const Box& bx, Real temp_Sr[3]; coriolis = false; - rotational_acceleration(loc, temp_vel, coriolis, temp_Sr); + rotational_acceleration(loc, temp_vel, omega, coord, coriolis, temp_Sr); edge_Sr[dir][edge] = temp_Sr[dir]; @@ -487,4 +449,3 @@ Castro::corrrsrc(const Box& bx, }); } - diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp index 22935bb4e9..380b2723f1 100644 --- a/Source/scf/scf_relax.cpp +++ b/Source/scf/scf_relax.cpp @@ -377,6 +377,7 @@ Castro::do_hscf_solve() { const auto *dx = geomdata.CellSize(); const auto *problo = geomdata.ProbLo(); + const auto coord = geomdata.Coord(); // The below assumes we are rotating on the z-axis. @@ -413,7 +414,8 @@ Castro::do_hscf_solve() scale = (1.0_rt - rr[0]) * (1.0_rt - rr[1]) * (1.0_rt - rr[2]); } - Real bernoulli_zone = scale * (phi_arr(i,j,k) + rotational_potential(r)); + auto omega = get_omega_vec(geomdata, j); + Real bernoulli_zone = scale * (phi_arr(i,j,k) + rotational_potential(r, omega, coord)); return {bernoulli_zone}; }); @@ -457,6 +459,7 @@ Castro::do_hscf_solve() const auto *dx = geomdata.CellSize(); const auto *problo = geomdata.ProbLo(); + const auto coord = geomdata.Coord(); GpuArray r = {0.0}; @@ -468,7 +471,8 @@ Castro::do_hscf_solve() r[2] = problo[2] + (static_cast(k) + 0.5_rt) * dx[2] - problem::center[2]; #endif - enthalpy_arr(i,j,k) = bernoulli - phi_arr(i,j,k) - rotational_potential(r); + auto omega = get_omega_vec(geomdata, j); + enthalpy_arr(i,j,k) = bernoulli - phi_arr(i,j,k) - rotational_potential(r, omega, coord); }); } @@ -640,6 +644,7 @@ Castro::do_hscf_solve() Real dM = 0.0, dK = 0.0, dU = 0.0, dE = 0.0; const auto* problo = geomdata.ProbLo(); + const auto coord = geomdata.Coord(); GpuArray r = {0.0}; @@ -655,7 +660,8 @@ Castro::do_hscf_solve() { dM = state_arr(i,j,k,URHO) * dV; - dK = rotational_potential(r) * dM; + auto omega = get_omega_vec(geomdata, j); + dK = rotational_potential(r, omega, coord) * dM; dU = 0.5_rt * phi_arr(i,j,k) * dM;