Remove grav_source_type and rot_source_type parameters

CHANGES.md

@@ -239,6 +239,8 @@
 
 # 21.01
 
+   * grav_source_type and rot_source_type have been removed. (#580)
+
    * The minimum C++ standard supported by Castro is now C++17. Most modern compilers
      support C++17; the notable exception is RHEL 7 and its derivatives like CentOS 7,
      where the default compiler is gcc 4.8. In that case a newer compiler must be loaded,

Docs/source/FlowChart.rst

@@ -421,9 +421,7 @@ In the code, the objective is to evolve the state from the old time,
                \Delta\phi^n = 4\pi G\rho^n,
 
       The construction of the form of the gravity source for the
-      momentum and energy equation is dependent on the parameter
-      ``castro.grav_source_type``. Full details of the gravity
-      solver are given in Chapter :ref:`ch:gravity`.
+      momentum and energy equation are given in Chapter :ref:`ch:gravity`.
 
 
    H. [``ROTATION``] rotation
@@ -432,9 +430,7 @@ In the code, the objective is to evolve the state from the old time,
       and the rotational acceleration (for use in the momentum
       equation). This includes the Coriolis and centrifugal terms in a
       constant-angular-velocity co-rotating frame. The form of the
-      rotational source that is constructed then depends on the
-      parameter ``castro.rot_source_type``. More details are
-      given in Chapter :ref:`ch:rotation`.
+      rotational source that is constructed is given in Chapter :ref:`ch:rotation`.
 
    The source terms here are evaluated using the post-burn state,
    :math:`\Ub^\star` (``Sborder``), and later corrected by using the

Docs/source/gravity.rst

@@ -568,47 +568,18 @@ where :math:`M_{enclosed}` has the same meaning as with the
 Hydrodynamics Source Terms
 ==========================
 
-There are several options to incorporate the effects of gravity into
-the hydrodynamics system. The main parameter here is
-``castro.grav_source_type``.
-
-- ``castro.grav_source_type`` = 1 : we use a standard
-  predictor-corrector formalism for updating the momentum and
-  energy. Specifically, our first update is equal to :math:`\Delta t
-  \times \mathbf{S}^n` , where :math:`\mathbf{S}^n` is the value of
-  the source terms at the old-time (which is usually called time-level
-  :math:`n`). At the end of the timestep, we do a corrector step where
-  we subtract off :math:`\Delta t / 2 \times \mathbf{S}^n` and add on
-  :math:`\Delta t / 2 \times \mathbf{S}^{n+1}`, so that at the end of
-  the timestep the source term is properly time centered.
-
-- ``castro.grav_source_type`` = 2 : we do something very similar
-  to 1. The major difference is that when evaluating the energy source
-  term at the new time (which is equal to :math:`\mathbf{u} \cdot
-  \mathbf{S}^{n+1}_{\rho \mathbf{u}}`, where the latter is the
-  momentum source term evaluated at the new time), we first update the
-  momentum, rather than using the value of :math:`\mathbf{u}` before
-  entering the gravity source terms. This permits a tighter coupling
-  between the momentum and energy update and we have seen that it
-  usually results in a more accurate evolution.
-
-- ``castro.grav_source_type`` = 3 : we do the same momentum update as
-  the previous two, but for the energy update, we put all of the work
-  into updating the kinetic energy alone. In particular, we explicitly
-  ensure that :math:`(\rho e)` remains the same, and update
-  :math:`(\rho K)` with the work due to gravity, adding the new kinetic
-  energy to the old internal energy to determine the final total gas
-  energy. The physical motivation is that work should be done on the
-  velocity, and should not directly update the temperature—only
-  indirectly through things like shocks.
-
-- ``castro.grav_source_type`` = 4 : the energy update is done in a
-  “conservative” fashion. The previous methods all evaluate the value
-  of the source term at the cell center, but this method evaluates the
-  change in energy at cell edges, using the hydrodynamical mass
-  fluxes, permitting total energy to be conserved (excluding possible
-  losses at open domain boundaries). See
-  :cite:`katzthesis` for some more details.
+We use a standard predictor-corrector formalism for updating the momentum.
+Our first update is equal to :math:`\Delta t \times \mathbf{S}^n`, where
+:math:`\mathbf{S}^n` is the value of the source terms at the old-time
+(which is usually called time-level :math:`n`). At the end of the timestep,
+we do a corrector step where we subtract off :math:`\Delta t / 2 \times \mathbf{S}^n`
+and add on :math:`\Delta t / 2 \times \mathbf{S}^{n+1}`, so that at the end of
+the timestep the source term is properly time centered.
+
+The energy update is done in an explicitly “conservative” fashion.
+We evaluate the change in energy at cell edges, using the hydrodynamical mass
+fluxes, permitting total energy to be conserved (excluding possible losses at
+open domain boundaries). See :cite:`katzthesis` for some more details.
 
 .. [2]
    Note: The GR

Docs/source/rotation.rst

@@ -68,9 +68,6 @@ The main parameters that affect rotation are:
    include the forcing from the time derivative of the rotation
    frequency (default: 1)
 
--  ``castro.rot_source_type`` : method of updating the
-   energy during a rotation update (default: 4)
-
 -  ``castro.implicit_rotation_update`` : for the Coriolis
    term, which mixes momenta in the source term, whether we should
    solve for the update implicitly (default: 1)
@@ -282,59 +279,15 @@ reference:
 Adding the forcing to the hydrodynamics
 =======================================
 
-There are several ways to incorporate the effect of the rotation
-forcing on the hydrodynamical evolution. We control this through the
-use of the runtime parameter castro.rot_source_type. This
-is an integer with values currently ranging from 1 through 4, and
-these values are all analogous to the way that gravity is used to
-update the momentum and energy. For the most part, the differences are
-in how the energy update is done:
-
-* ``castro.rot_source_type = 1`` : we use a standard
-  predictor-corrector formalism for updating the momentum and
-  energy. Specifically, our first update is equal to :math:`\Delta t \times \mathbf{S}^n` ,
-  where :math:`\mathbf{S}^n` is the value of
-  the source terms at the old-time (which is usually called time-level
-  :math:`n`). At the end of the timestep, we do a corrector step where
-  we subtract off :math:`\Delta t / 2 \times \mathbf{S}^n` and add on
-  :math:`\Delta t / 2 \times \mathbf{S}^{n+1}`, so that at the end of
-  the timestep the source term is properly time centered.
-
-* ``castro.rot_source_type = 2`` : we do something very similar
-  to 1. The major difference is that when evaluating the energy source
-  term at the new time (which is equal to
-  :math:`\mathbf{u} \cdot \mathbf{S}^{n+1}_{\rho \mathbf{u}}`, where the latter is the
-  momentum source term evaluated at the new time), we first update the
-  momentum, rather than using the value of :math:`\mathbf{u}` before
-  entering the rotation source terms. This permits a tighter coupling
-  between the momentum and energy update and we have seen that it
-  usually results in a more accurate evolution.
-
-* ``castro.rot_source_type = 3`` : we do the same momentum update as
-  the previous two, but for the energy update, we put all of the work
-  into updating the kinetic energy alone. In particular, we explicitly
-  ensure that :math:`(rho e)` maintains the same, and update
-  :math:`(rho K)` with the work due to rotation, adding the new
-  kinetic energy to the old internal energy to determine the final
-  total gas energy. The physical motivation is that work should be
-  done on the velocity, and should not directly update the temperature
-  – only indirectly through things like shocks.
-
-* ``castro.rot_source_type = 4`` : the energy update is done in a
-   “conservative” fashion. The previous methods all evaluate the value
-   of the source term at the cell center, but this method evaluates
-   the change in energy at cell edges, using the hydrodynamical mass
-   fluxes, permitting total energy to be conserved (excluding possible
-   losses at open domain boundaries). Additionally, the velocity
-   update is slightly different—for the corrector step, we note that
-   there is an implicit coupling between the velocity components, and
-   we directly solve this coupled equation, which results in a
-   slightly better coupling and a more accurate evolution.
-
-The other major option is ``castro.implicit_rotation_update``.
-This does the update of the Coriolis term in the momentum equation
-implicitly (e.g., the velocity in the Coriolis force for the zone
-depends on the updated momentum). The energy update is unchanged.
-
-A detailed discussion of these options and some verification
+The momentum update is done using a standard cell-centered formulation.
+The energy update is done in a “conservative” fashion. We evaluate
+the change in energy at cell edges, using the hydrodynamical mass
+fluxes, permitting total energy to be conserved (excluding possible
+losses at open domain boundaries). Additionally, for the corrector step
+in the velocity update, we note that there is an implicit coupling between
+the velocity components, and we can directly solve this coupled equation,
+which results in a slightly better coupling and a more accurate evolution.
+This is controlled using ``castro.implicit_rotation_update``.
+
+A detailed discussion of the rotational forcing and some verification
 tests is presented in :cite:`katz:2016`.

Source/driver/_cpp_parameters

@@ -484,10 +484,6 @@ do_grav                      int          -1
 # to we recompute the center used for the multipole gravity solve each step?
 moving_center                int           0
 
-# determines how the gravitational source term is added to the momentum and
-# energy state variables.
-grav_source_type             int           4
-
 # permits rotation calculation to be turned on and off
 do_rotation                  int          -1
 
@@ -500,10 +496,6 @@ rotation_include_centrifugal int           1                  n        ROTATION
 # permits the Coriolis terms in the rotation to be turned on and off
 rotation_include_coriolis    int           1                  n        ROTATION
 
-# determines how the rotation source terms are added to the momentum and
-# energy equations
-rot_source_type              int           4                  n        ROTATION
-
 # we can do a implicit solution of the rotation update to allow
 # for better coupling of the Coriolis terms
 implicit_rotation_update     int           1                  n        ROTATION

Source/gravity/Castro_gravity.cpp

@@ -258,8 +258,6 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in,
     GeometryData geomdata = geom.data();
 #endif
 
-    AMREX_ALWAYS_ASSERT(castro::grav_source_type >= 1 && castro::grav_source_type <= 4);
-
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
@@ -288,13 +286,6 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in,
                 src[n] = 0.0_rt;
             }
 
-            // Gravitational source options for how to add the work to (rho E):
-            // grav_source_type =
-            // 1: Original version ("does work")
-            // 2: Modification of type 1 that updates the momentum before constructing the energy corrector
-            // 3: Puts all gravitational work into KE, not (rho e)
-            // 4: Conservative energy formulation
-
             Real rho    = uold(i,j,k,URHO);
             Real rhoInv = 1.0_rt / rho;
 
@@ -330,33 +321,15 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in,
             }
 #endif
 
-            Real SrE;
-
-            if (castro::grav_source_type == 1 || castro::grav_source_type == 2) {
-
-                // Src = rho u dot g, evaluated with all quantities at t^n
-
-                SrE = (uold(i,j,k,UMX) * Sr[0] + uold(i,j,k,UMY) * Sr[1] + uold(i,j,k,UMZ) * Sr[2]) * rhoInv;
-
-            } else if (castro::grav_source_type == 3) {
-
-                Real new_ke = 0.5_rt * (snew[UMX] * snew[UMX] + snew[UMY] * snew[UMY] + snew[UMZ] * snew[UMZ]) * rhoInv;
-                SrE = new_ke - old_ke;
-
-            } else if (castro::grav_source_type == 4) {
+            // Our conservative energy formulation does not strictly require
+            // any energy source-term here, because it depends only on the
+            // fluid motions from the hydrodynamical fluxes which we will only
+            // have when we get to the 'corrector' step. Nevertheless we add a
+            // predictor energy source term in the way that the other methods
+            // do, for consistency. We will fully subtract this predictor value
+            // during the corrector step, so that the final result is correct.
 
-                // The conservative energy formulation does not strictly require
-                // any energy source-term here, because it depends only on the
-                // fluid motions from the hydrodynamical fluxes which we will only
-                // have when we get to the 'corrector' step. Nevertheless we add a
-                // predictor energy source term in the way that the other methods
-                // do, for consistency. We will fully subtract this predictor value
-                // during the corrector step, so that the final result is correct.
-                // Here we use the same approach as grav_source_type == 2.
-
-                SrE = (uold(i,j,k,UMX) * Sr[0] + uold(i,j,k,UMY) * Sr[1] + uold(i,j,k,UMZ) * Sr[2]) * rhoInv;
-
-            }
+            Real SrE = (uold(i,j,k,UMX) * Sr[0] + uold(i,j,k,UMY) * Sr[1] + uold(i,j,k,UMZ) * Sr[2]) * rhoInv;
 
             src[UEDEN] = SrE;
 
@@ -418,8 +391,6 @@ void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old,
     GeometryData geomdata = geom.data();
 #endif
 
-    AMREX_ALWAYS_ASSERT(castro::grav_source_type >= 1 && castro::grav_source_type <= 4);
-
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
@@ -445,13 +416,6 @@ void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old,
 
                 Real hdtInv = 0.5_rt / dt;
 
-                // Gravitational source options for how to add the work to (rho E):
-                // grav_source_type =
-                // 1: Original version ("does work")
-                // 2: Modification of type 1 that updates the U before constructing SrEcorr
-                // 3: Puts all gravitational work into KE, not (rho e)
-                // 4: Conservative gravity approach (discussed in first white dwarf merger paper).
-
                 Real rhoo    = uold(i,j,k,URHO);
                 Real rhooinv = 1.0_rt / uold(i,j,k,URHO);
 
@@ -528,73 +492,41 @@ void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old,
 
                 // Correct energy
 
-                Real SrEcorr;
-
-                if (castro::grav_source_type == 1) {
-
-                    // If grav_source_type == 1, then we calculated SrEcorr before updating the velocities.
+                // First, subtract the predictor step we applied earlier.
 
-                    SrEcorr = 0.5_rt * (SrE_new - SrE_old);
+                Real SrEcorr = - SrE_old;
 
-                } else if (castro::grav_source_type == 2) {
+                // For an explanation of this approach, see wdmerger paper I.
+                // The main idea is that we are evaluating the change of the
+                // potential energy at zone edges and applying that in an equal
+                // and opposite sense to the gas energy. The physics is described
+                // in Section 2.4; we are using a version of the formula similar to
+                // Equation 94 in Springel (2010) based on the gradient rather than
+                // the potential because the gradient-version works for all forms
+                // of gravity we use, some of which do not explicitly calculate phi.
 
-                    // For this source type, we first update the momenta
-                    // before we calculate the energy source term.
+                // Construct the time-averaged edge-centered gravity.
 
-                    for (int n = 0; n < 3; ++n) {
-                        vnew[n] = snew[UMX+n] * rhoninv;
-                    }
-                    SrE_new = vnew[0] * Sr_new[0] + vnew[1] * Sr_new[1] + vnew[2] * Sr_new[2];
-
-                    SrEcorr = 0.5_rt * (SrE_new - SrE_old);
-
-                } else if (castro::grav_source_type == 3) {
-
-                    // Instead of calculating the energy source term explicitly,
-                    // we simply update the kinetic energy.
-
-                    Real new_ke = 0.5_rt * (snew[UMX] * snew[UMX] + snew[UMY] * snew[UMY] + snew[UMZ] * snew[UMZ]) * rhoninv;
-                    SrEcorr = new_ke - old_ke;
-
-                } else if (castro::grav_source_type == 4) {
-
-                    // First, subtract the predictor step we applied earlier.
-
-                    SrEcorr = - SrE_old;
-
-                    // For an explanation of this approach, see wdmerger paper I.
-                    // The main idea is that we are evaluating the change of the
-                    // potential energy at zone edges and applying that in an equal
-                    // and opposite sense to the gas energy. The physics is described
-                    // in Section 2.4; we are using a version of the formula similar to
-                    // Equation 94 in Springel (2010) based on the gradient rather than
-                    // the potential because the gradient-version works for all forms
-                    // of gravity we use, some of which do not explicitly calculate phi.
-
-                    // Construct the time-averaged edge-centered gravity.
-
-                    GpuArray<Real, 3> g;
-                    for (int n = 0; n < 3; ++n) {
-                        g[n] = 0.5_rt * (gnew(i,j,k,n) + gold(i,j,k,n));
-                    }
-
-                    Real gxl = 0.5_rt * (g[0] + 0.5_rt * (gnew(i-1*dg0,j,k,0) + gold(i-1*dg0,j,k,0)));
-                    Real gxr = 0.5_rt * (g[0] + 0.5_rt * (gnew(i+1*dg0,j,k,0) + gold(i+1*dg0,j,k,0)));
+                GpuArray<Real, 3> g;
+                for (int n = 0; n < 3; ++n) {
+                    g[n] = 0.5_rt * (gnew(i,j,k,n) + gold(i,j,k,n));
+                }
 
-                    Real gyl = 0.5_rt * (g[1] + 0.5_rt * (gnew(i,j-1*dg1,k,1) + gold(i,j-1*dg1,k,1)));
-                    Real gyr = 0.5_rt * (g[1] + 0.5_rt * (gnew(i,j+1*dg1,k,1) + gold(i,j+1*dg1,k,1)));
+                Real gxl = 0.5_rt * (g[0] + 0.5_rt * (gnew(i-1*dg0,j,k,0) + gold(i-1*dg0,j,k,0)));
+                Real gxr = 0.5_rt * (g[0] + 0.5_rt * (gnew(i+1*dg0,j,k,0) + gold(i+1*dg0,j,k,0)));
 
-                    Real gzl = 0.5_rt * (g[2] + 0.5_rt * (gnew(i,j,k-1*dg2,2) + gold(i,j,k-1*dg2,2)));
-                    Real gzr = 0.5_rt * (g[2] + 0.5_rt * (gnew(i,j,k+1*dg2,2) + gold(i,j,k+1*dg2,2)));
+                Real gyl = 0.5_rt * (g[1] + 0.5_rt * (gnew(i,j-1*dg1,k,1) + gold(i,j-1*dg1,k,1)));
+                Real gyr = 0.5_rt * (g[1] + 0.5_rt * (gnew(i,j+1*dg1,k,1) + gold(i,j+1*dg1,k,1)));
 
-                    SrEcorr += hdtInv * (flux0(i      ,j,k) * gxl * dx[0] +
-                                         flux0(i+1*dg0,j,k) * gxr * dx[0] +
-                                         flux1(i,j      ,k) * gyl * dx[1] +
-                                         flux1(i,j+1*dg1,k) * gyr * dx[1] +
-                                         flux2(i,j,k      ) * gzl * dx[2] +
-                                         flux2(i,j,k+1*dg2) * gzr * dx[2]) / vol(i,j,k);
+                Real gzl = 0.5_rt * (g[2] + 0.5_rt * (gnew(i,j,k-1*dg2,2) + gold(i,j,k-1*dg2,2)));
+                Real gzr = 0.5_rt * (g[2] + 0.5_rt * (gnew(i,j,k+1*dg2,2) + gold(i,j,k+1*dg2,2)));
 
-                }
+                SrEcorr += hdtInv * (flux0(i      ,j,k) * gxl * dx[0] +
+                                     flux0(i+1*dg0,j,k) * gxr * dx[0] +
+                                     flux1(i,j      ,k) * gyl * dx[1] +
+                                     flux1(i,j+1*dg1,k) * gyr * dx[1] +
+                                     flux2(i,j,k      ) * gzl * dx[2] +
+                                     flux2(i,j,k+1*dg2) * gzr * dx[2]) / vol(i,j,k);
 
                 src[UEDEN] = SrEcorr;