Add passive scalars

inputs/hydro/athinput.slotted_cylinder2d

@@ -0,0 +1,70 @@
+<comment>
+problem   = advection of a 2D slotted cylinder passive scalar profile
+reference = Zalesak, S.T., JCP 31, 335-362 (1979) (Figure 11, section VII)
+reference = Colella, P., et al, JCP 230, 2952–2976 (2011) (Figure 13, section 4.4.5)
+configure = --prob=slotted_cylinder --eos=isothermal --nscalars=1
+
+<job>
+problem_id = SlottedCylinder   # problem ID: basename of output filenames
+
+<output1>
+file_type  = hst       # History data dump
+dt         = 0.01      # time increment between outputs
+data_format = %.17g    # round trip conversion: binary64 <-> decimal
+
+<output2>
+file_type   = vtk      # Binary data dump
+variable    = prim     # variables to be output
+dt          = 0.01     # time increment between outputs
+
+<time>
+cfl_number = 0.4        # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1         # cycle limit
+tlim       = 1.0        # time limit
+integrator  = vl2       # time integration algorithm
+xorder      = 2         # order of spatial reconstruction
+ncycle_out  = 1         # interval for stdout summary info
+correct_ic  = false     # correct midpoint assumption in initial condition
+correct_err = false     # correct midpoint assumption in analytic solution
+
+<mesh>
+nx1        = 100       # Number of zones in X1-direction
+x1min      = -0.5      # minimum value of X1
+x1max      = 0.5       # maximum value of X1
+ix1_bc     = periodic  # inner-X1 boundary flag
+ox1_bc     = periodic  # outer-X1 boundary flag
+
+nx2        = 100       # Number of zones in X2-direction
+x2min      = -0.5      # minimum value of X2
+x2max      = 0.5       # maximum value of X2
+ix2_bc     = periodic  # inner-X2 boundary flag
+ox2_bc     = periodic  # outer-X2 boundary flag
+
+nx3        = 1         # Number of zones in X3-direction
+x3min      = -0.5      # minimum value of X3
+x3max      = 0.5       # maximum value of X3
+ix3_bc     = periodic  # inner-X2 boundary flag
+ox3_bc     = periodic  # outer-X2 boundary flag
+
+num_threads = 1        # maximum number of OMP threads
+refinement  = none
+
+<meshblock>
+nx1        = 100        # Number of zones in X1-direction
+nx2        = 100        # Number of zones in X2-direction
+nx3        = 1        # Number of zones in X3-direction
+
+<hydro>
+iso_sound_speed = 1.0     # isothermal sound speed
+active    = background    # compute Hydro fluxes, etc. but do NOT evolve the fluid
+
+<problem>
+compute_error = true   # set value to 'true' to compute L1 error compared to initial data
+radius     = 0.15      # Radius of the cylinder
+center_x1  = 0.00      # x1 position of cylinder center
+center_x2  = 0.25      # x2 position of cylinder center
+omega      = 1.0       # angular velocity of the advecting cylinder
+omega_x1   = 0.00      # x1 position of axis of rotation
+omega_x2   = 0.00      # x2 position of axis of rotation
+s_width    = 0.05      # Width of the slot
+s_height   = 0.25      # Height of the slot

inputs/hydro/athinput.slotted_cylinder2d_refined

@@ -0,0 +1,99 @@
+<comment>
+problem   = advection of a 2D slotted cylinder passive scalar profile (with SMR or AMR)
+reference = Zalesak, S.T., JCP 31, 335-362 (1979) (Figure 11, section VII)
+reference = Colella, P., et al, JCP 230, 2952–2976 (2011) (Figure 13, section 4.4.5)
+configure = --prob=slotted_cylinder --eos=isothermal --nscalars=1 -hdf5
+
+<job>
+problem_id = SlottedCylinder   # problem ID: basename of output filenames
+
+<output1>
+file_type  = hst       # History data dump
+dt         = 0.01      # time increment between outputs
+data_format = %.17g    # round trip conversion: binary64 <-> decimal
+
+<output2>
+file_type   = hdf5     # HDF5 data dump
+variable    = prim     # variables to be output
+dt          = 0.01     # time increment between outputs
+
+<time>
+cfl_number = 0.4        # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1         # cycle limit
+tlim       = 1.0        # time limit
+integrator  = vl2       # time integration algorithm
+xorder      = 2         # order of spatial reconstruction
+ncycle_out  = 1         # interval for stdout summary info
+correct_ic  = false     # correct midpoint assumption in initial condition
+correct_err = false     # correct midpoint assumption in analytic solution
+
+<mesh>
+nx1        = 100       # Number of zones in X1-direction
+x1min      = -0.5      # minimum value of X1
+x1max      = 0.5       # maximum value of X1
+ix1_bc     = periodic  # inner-X1 boundary flag
+ox1_bc     = periodic  # outer-X1 boundary flag
+
+nx2        = 100       # Number of zones in X2-direction
+x2min      = -0.5      # minimum value of X2
+x2max      = 0.5       # maximum value of X2
+ix2_bc     = periodic  # inner-X2 boundary flag
+ox2_bc     = periodic  # outer-X2 boundary flag
+
+nx3        = 1         # Number of zones in X3-direction
+x3min      = -0.5      # minimum value of X3
+x3max      = 0.5       # maximum value of X3
+ix3_bc     = periodic  # inner-X2 boundary flag
+ox3_bc     = periodic  # outer-X2 boundary flag
+
+num_threads = 1        # maximum number of OMP threads
+refinement  = none
+numlevel    = 4
+deref_count = 5
+
+<meshblock>
+nx1        = 4        # Number of zones in X1-direction
+nx2        = 4        # Number of zones in X2-direction
+nx3        = 1        # Number of zones in X3-direction
+
+# The following refinement regions are used in regression tests.  By default,
+# they are ignored.  To use, set MeshBlock size to 4x4x1. (or 20x20x1, but then the default
+# static refinement regions below will be larger than indicated in the comments)
+
+# enclose the cylinder in a 0.3x0.3 square region (assuming default <problem> parameter values)
+<refinement1>
+x1min = -0.15
+x1max = 0.15
+x2min = 0.1
+x2max = 0.4
+x3min = -0.5
+x3max = 0.5
+level = 1
+
+# cover the slot with a 0.08x0.3 rectangular region (assuming default <problem> parameter values)
+<refinement2>
+# for 4x4x1 MeshBlocks with 5 cell slot on root grid, this will refine all blocks wholly
+# within slot, +1x on both sides
+x1min = -0.040
+x1max = 0.040
+x2min = 0.1
+x2max = 0.4
+x3min = -0.5
+x3max = 0.5
+level = 2
+
+<hydro>
+iso_sound_speed = 1.0     # isothermal sound speed
+active    = background    # compute Hydro fluxes, etc. but do NOT evolve the fluid
+
+<problem>
+compute_error = true   # set value to 'true' to compute L1 error compared to initial data
+radius     = 0.15      # Radius of the cylinder
+center_x1  = 0.00      # x1 position of cylinder center
+center_x2  = 0.25      # x2 position of cylinder center
+omega      = 1.0       # angular velocity of the advecting cylinder
+omega_x1   = 0.00      # x1 position of axis of rotation
+omega_x2   = 0.00      # x2 position of axis of rotation
+s_width    = 0.05      # Width of the slot
+s_height   = 0.25      # Height of the slot
+thr        = 0.2       # scalar gradient threshold for AMR

src/athena.hpp

@@ -78,7 +78,7 @@ struct LogicalLocation { // aggregate and POD type
 //! \struct RegionSize
 //  \brief physical size and number of cells in a Mesh or a MeshBlock
 
-struct RegionSize {
+struct RegionSize {  // aggregate and POD type; do NOT reorder member declarations:
   Real x1min, x2min, x3min;
   Real x1max, x2max, x3max;
   Real x1rat, x2rat, x3rat; // ratio of x(i)/x(i-1)
@@ -158,6 +158,8 @@ enum CoordinateDirection {X1DIR=0, X2DIR=1, X3DIR=2};
 enum class BoundaryQuantity {cc, fc, cc_flcor, fc_flcor, mggrav, mggrav_f};
 enum class HydroBoundaryQuantity {cons, prim};
 enum class BoundaryCommSubset {mesh_init, gr_amr, all};
+// TODO(felker): consider generalizing/renaming to QuantityFormulation
+enum class FluidFormulation {evolve, background, disabled}; // rename background -> fixed?
 
 //----------------------------------------------------------------------------------------
 // function pointer prototypes for user-defined modules set at runtime

src/bvals/bvals_base.cpp

@@ -81,7 +81,7 @@ BoundaryBase::BoundaryBase(Mesh *pm, LogicalLocation iloc, RegionSize isize,
   loc = iloc;
   block_size_ = isize;
   pmy_mesh_ = pm;
-  if (called_ == false) {
+  if (!called_) {
     int dim = 1;
     if (block_size_.nx2 > 1) dim = 2;
     if (block_size_.nx3 > 1) dim = 3;

src/bvals/bvals_var.cpp

@@ -1,3 +1,4 @@
+
 //========================================================================================
 // Athena++ astrophysical MHD code
 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
@@ -223,7 +224,7 @@ bool BoundaryVariable::ReceiveBoundaryBuffers() {
         int test;
         MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &test, MPI_STATUS_IGNORE);
         MPI_Test(&(bd_var_.req_recv[nb.bufid]), &test, MPI_STATUS_IGNORE);
-        if (static_cast<bool>(test) == false) {
+        if (!static_cast<bool>(test)) {
           bflag = false;
           continue;
         }

src/bvals/cc/bvals_cc.cpp

@@ -534,7 +534,7 @@ void CellCenteredBoundaryVariable::SetupPersistentMPI() {
   MeshBlock* pmb = pmy_block_;
   int &mylevel = pmb->loc.level;
 
-  int f2 = pmy_mesh_->f2_, f3 = pmy_mesh_->f3_;
+  int f2 = pmy_mesh_->f2, f3 = pmy_mesh_->f3;
   int cng, cng1, cng2, cng3;
   cng  = cng1 = pmb->cnghost;
   cng2 = cng*f2;

src/bvals/cc/bvals_shear_cc.cpp

@@ -332,7 +332,7 @@ bool CellCenteredBoundaryVariable::ReceiveShearingBoxBoundaryBuffers() {
             MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &test,
                        MPI_STATUS_IGNORE);
             MPI_Test(&shear_bd_var_[upper].req_recv[n], &test, MPI_STATUS_IGNORE);
-            if (static_cast<bool>(test) == false) {
+            if (!static_cast<bool>(test)) {
               flag[upper] = false;
               continue;
             }

src/bvals/cc/flux_correction_cc.cpp

@@ -173,7 +173,7 @@ bool CellCenteredBoundaryVariable::ReceiveFluxCorrection() {
           MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &test,
                      MPI_STATUS_IGNORE);
           MPI_Test(&(bd_var_flcor_.req_recv[nb.bufid]), &test, MPI_STATUS_IGNORE);
-          if (static_cast<bool>(test) == false) {
+          if (!static_cast<bool>(test)) {
             bflag = false;
             continue;
           }

src/bvals/cc/mg/bvals_mg.cpp

@@ -44,15 +44,15 @@ class MultigridDriver;
 
 MGBoundaryValues::MGBoundaryValues(Multigrid *pmg, BoundaryFlag *input_bcs,
                                    MGBoundaryFunc *MGBoundary)
-    : BoundaryBase(pmg->pmy_driver_->pmy_mesh_, pmg->loc_, pmg->size_, input_bcs) {
-  pmy_mg_ = pmg;
+    : BoundaryBase(pmg->pmy_driver_->pmy_mesh_, pmg->loc_, pmg->size_, input_bcs),
+      pmy_mg_(pmg) {
 #ifdef MPI_PARALLEL
   mgcomm_ = pmg->pmy_driver_->MPI_COMM_MULTIGRID;
   // currently assuming that Multigrid gravity is the only "int physid" associated with
   // the owning MultigridDriver object:
   mg_grav_phys_id_ = pmg->pmy_driver_->mg_phys_id_;
 #endif
-  if (pmy_mg_->root_flag_ == true) {
+  if (pmy_mg_->root_flag_) {
     for (int i=0; i<6; i++)
       MGBoundaryFunction_[i] = MGBoundary[i];
   } else {
@@ -441,7 +441,7 @@ bool MGBoundaryValues::ReceiveMultigridBoundaryBuffers(AthenaArray<Real> &dst,
         int test;
         MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, mgcomm_, &test, MPI_STATUS_IGNORE);
         MPI_Test(&(pbd->req_recv[nb.bufid]), &test, MPI_STATUS_IGNORE);
-        if (static_cast<bool>(test) == false) {
+        if (!static_cast<bool>(test)) {
           bflag = false;
           continue;
         }

src/bvals/fc/bvals_fc.cpp

@@ -249,7 +249,7 @@ int FaceCenteredBoundaryVariable::ComputeVariableBufferSize(const NeighborIndexe
   int nx1 = pmb->block_size.nx1;
   int nx2 = pmb->block_size.nx2;
   int nx3 = pmb->block_size.nx3;
-  int f2 = pmy_mesh_->f2_, f3 = pmy_mesh_->f3_;
+  int f2 = pmy_mesh_->f2, f3 = pmy_mesh_->f3;
   int cng1, cng2, cng3;
   cng1 = cng;
   cng2 = cng*f2;
@@ -1207,7 +1207,7 @@ void FaceCenteredBoundaryVariable::SetupPersistentMPI() {
   int nx3 = pmb->block_size.nx3;
   int &mylevel = pmb->loc.level;
 
-  int f2 = pmy_mesh_->f2_, f3 = pmy_mesh_->f3_;
+  int f2 = pmy_mesh_->f2, f3 = pmy_mesh_->f3;
   int cng, cng1, cng2, cng3;
   cng  = cng1 = pmb->cnghost;
   cng2 = cng*f2;
@@ -1331,7 +1331,7 @@ void FaceCenteredBoundaryVariable::SetupPersistentMPI() {
       if (nb.snb.level == mylevel) { // the same level
         if ((nb.ni.type == NeighborConnect::face)
             || ((nb.ni.type == NeighborConnect::edge)
-                && (edge_flag_[nb.eid] == true))) {
+                && (edge_flag_[nb.eid]))) {
           tag = pbval_->CreateBvalsMPITag(nb.snb.lid, nb.targetid, fc_flx_phys_id_);
           if (bd_var_flcor_.req_send[nb.bufid] != MPI_REQUEST_NULL)
             MPI_Request_free(&bd_var_flcor_.req_send[nb.bufid]);
@@ -1416,7 +1416,7 @@ void FaceCenteredBoundaryVariable::StartReceiving(BoundaryCommSubset phase) {
         if ((nb.snb.level > mylevel) ||
             ((nb.snb.level == mylevel) && ((nb.ni.type == NeighborConnect::face)
                                            || ((nb.ni.type == NeighborConnect::edge)
-                                               && (edge_flag_[nb.eid] == true)))))
+                                               && (edge_flag_[nb.eid])))))
           MPI_Start(&(bd_var_flcor_.req_recv[nb.bufid]));
       }
     }
@@ -1466,7 +1466,7 @@ void FaceCenteredBoundaryVariable::ClearBoundary(BoundaryCommSubset phase) {
           else if ((nb.snb.level == mylevel)
                    && ((nb.ni.type == NeighborConnect::face)
                        || ((nb.ni.type == NeighborConnect::edge)
-                           && (edge_flag_[nb.eid] == true))))
+                           && (edge_flag_[nb.eid]))))
             MPI_Wait(&(bd_var_flcor_.req_send[nb.bufid]), MPI_STATUS_IGNORE);
         }
       }

src/bvals/fc/bvals_shear_emf.cpp

@@ -272,7 +272,7 @@ bool FaceCenteredBoundaryVariable::ReceiveEMFShearingBoxBoundaryCorrection() {
             MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &test,
                        MPI_STATUS_IGNORE);
             MPI_Test(&shear_bd_emf_[upper].req_recv[n], &test, MPI_STATUS_IGNORE);
-            if (static_cast<bool>(test) == false) {
+            if (!static_cast<bool>(test)) {
               flag[upper] = false;
               continue;
             }

src/bvals/fc/bvals_shear_fc.cpp

@@ -330,7 +330,7 @@ bool FaceCenteredBoundaryVariable::ReceiveShearingBoxBoundaryBuffers() {
             MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &test,
                        MPI_STATUS_IGNORE);
             MPI_Test(&shear_bd_var_[upper].req_recv[n], &test, MPI_STATUS_IGNORE);
-            if (static_cast<bool>(test) == false) {
+            if (!static_cast<bool>(test)) {
               flag[upper] = false;
               continue;
             }