Completely removed the gravity flux formula

configure.py

@@ -123,12 +123,6 @@
                     default=False,
                     help='enable super-time-stepping')
 
-# -gflux argument
-parser.add_argument('-gflux',
-                    action='store_true',
-                    default=False,
-                    help='enable gravity momentum flux')
-
 # -s argument
 parser.add_argument('-s',
                     action='store_true',
@@ -424,12 +418,6 @@ def c_to_cpp(arg):
 else:
     definitions['STS_ENABLED'] = '0'
 
-# -gflux argument
-if args['gflux']:
-    definitions['GRAVITY_FLUX_ENABLED'] = '1'
-else:
-    definitions['GRAVITY_FLUX_ENABLED'] = '0'
-
 # -s, -g, and -t arguments
 definitions['RELATIVISTIC_DYNAMICS'] = '1' if args['s'] or args['g'] else '0'
 definitions['GENERAL_RELATIVITY'] = '1' if args['g'] else '0'
@@ -802,7 +790,6 @@ def c_to_cpp(arg):
 print('  General relativity:         ' + ('ON' if args['g'] else 'OFF'))
 print('  Frame transformations:      ' + ('ON' if args['t'] else 'OFF'))
 print('  Self-Gravity:               ' + self_grav_string)
-print('  Gravity Fluxes              ' + ('ON' if args['gflux'] else 'OFF'))
 print('  Super-Time-Stepping:        ' + ('ON' if args['sts'] else 'OFF'))
 print('  Debug flags:                ' + ('ON' if args['debug'] else 'OFF'))
 print('  Code coverage flags:        ' + ('ON' if args['coverage'] else 'OFF'))

src/defs.hpp.in

@@ -45,9 +45,6 @@
 // include self gravity? default=0 (false)
 #define SELF_GRAVITY_ENABLED @SELF_GRAVITY_ENABLED@
 
-// use gravity fluxes? default=0 (false)
-#define GRAVITY_FLUX_ENABLED @GRAVITY_FLUX_ENABLED@
-
 // include radiative transfer? default=0 (false)
 #define RADIATION_ENABLED 0
 

src/gravity/gravity.cpp

@@ -32,7 +32,6 @@ Gravity::Gravity(MeshBlock *pmb, ParameterInput *pin) :
     pmy_block(pmb), phi(pmb->ncells3, pmb->ncells2, pmb->ncells1),
     empty_flux{AthenaArray<Real>(), AthenaArray<Real>(), AthenaArray<Real>()},
     four_pi_G(pmb->pmy_mesh->four_pi_G_),
-    grav_mean_rho(pmb->pmy_mesh->grav_mean_rho_),
     gbvar(pmb, &phi, nullptr, empty_flux) {
   if (four_pi_G == 0.0) {
     std::stringstream msg;
@@ -43,15 +42,6 @@ Gravity::Gravity(MeshBlock *pmb, ParameterInput *pin) :
     return;
   }
 
-  if (grav_mean_rho == -1.0 && GRAVITY_FLUX_ENABLED) {
-    std::stringstream msg;
-    msg << "### FATAL ERROR in Gravity::Gravity" << std::endl
-        << "Background Mean Density must be set in the Mesh::InitUserMeshData "
-        << "using the SetMeanDensity function when using gravity flux." << std::endl;
-    ATHENA_ERROR(msg);
-    return;
-  }
-
   // using Gravity as an example of: containing full object members instead of pointer
   // memebers, construting BoundaryVariaable composite obj (no default ctor) in Gravity
   // ctor initializer list, avoiding dynamically-managed memory and the need for a

src/gravity/gravity.hpp

@@ -37,7 +37,6 @@ class Gravity {
   AthenaArray<Real> phi;   // gravitational potential
   AthenaArray<Real> empty_flux[3];
   Real gconst, four_pi_G;
-  Real grav_mean_rho;
   bool srcterm;
 
   // TODO(felker): consider creating a CellCentered.. derived class, and changing to

src/gravity/mg_gravity.cpp

@@ -127,7 +127,6 @@ void MGGravityDriver::Solve(int stage) {
   for (Multigrid* pmg : vmg_) {
     Gravity *pgrav = pmg->pmy_block_->pgrav;
     pmg->RetrieveResult(pgrav->phi, 0, NGHOST);
-    pgrav->grav_mean_rho = mean_rho;
   }
   return;
 }

src/hydro/calculate_fluxes.cpp

@@ -372,9 +372,6 @@ void Hydro::CalculateFluxes(AthenaArray<Real> &w, FaceField &b,
     } // end if (order == 4)
   }
 
-  if (SELF_GRAVITY_ENABLED && GRAVITY_FLUX_ENABLED)
-    AddGravityFlux();
-
   if (!STS_ENABLED)
     AddDiffusionFluxes();