Fully-implicit cosmic-ray diffusion module using Multigrid

Makefile.in

@@ -62,6 +62,7 @@ SRC_FILES := $(wildcard src/*.cpp) \
              $(wildcard src/nr_radiation/implicit/*.cpp) \
              $(wildcard src/cr/*.cpp) \
              $(wildcard src/cr/integrators/*.cpp) \
+             $(wildcard src/crdiffusion/*.cpp) \
 	     src/hydro/rsolvers/$(RSOLVER_DIR)$(RSOLVER_FILE) \
 	     $(wildcard src/inputs/*.cpp) \
 	     $(wildcard src/mesh/*.cpp) \

configure.py

@@ -46,6 +46,7 @@
 #   -nr_radiation        turn on non-relativistic radiation transport
 #   -implicit_radiation  implicit radiation transport module
 #   -cr                  enable cosmic ray transport
+#   -crdiff              enable cosmic ray diffusion with Multigrid
 # ----------------------------------------------------------------------------------------
 
 # Modules
@@ -267,6 +268,12 @@
                     default=False,
                     help='enable cosmic ray transport')
 
+# -cosmic ray diffusion argument
+parser.add_argument('-crdiff',
+                    action='store_true',
+                    default=False,
+                    help='enable implicit cosmic ray diffusion')
+
 # The main choices for --cxx flag, using "ctype[-suffix]" formatting, where "ctype" is the
 # major family/suite/group of compilers and "suffix" may represent variants of the
 # compiler version and/or predefined sets of compiler options. The C++ compiler front ends
@@ -558,6 +565,12 @@ def c_to_cpp(arg):
 else:
     definitions['CR_ENABLED'] = '0'
 
+# -crdiff argument
+if args['crdiff']:
+    definitions['CRDIFFUSION_ENABLED'] = '1'
+else:
+    definitions['CRDIFFUSION_ENABLED'] = '0'
+
 
 # --cxx=[name] argument
 if args['cxx'] == 'g++':
@@ -1025,6 +1038,7 @@ def output_config(opt_descr, opt_choice, filehandle=None):
 output_config('Radiative Transfer', ('ON' if args['nr_radiation'] else 'OFF'), flog)
 output_config('Implicit Radiation', ('ON' if args['implicit_radiation'] else 'OFF'), flog)
 output_config('Cosmic Ray Transport', ('ON' if args['cr'] else 'OFF'), flog)
+output_config('Cosmic Ray Diffusion', ('ON' if args['crdiff'] else 'OFF'), flog)
 output_config('Frame transformations', ('ON' if args['t'] else 'OFF'), flog)
 output_config('Self-Gravity', self_grav_string, flog)
 output_config('Super-Time-Stepping', ('ON' if args['sts'] else 'OFF'), flog)

inputs/cosmic_ray/athinput.cr_diffusion_mg

@@ -0,0 +1,75 @@
+<comment>
+problem   = CR diffusion test with Multigrid
+reference = 
+configure = --prob=cr_diffusion_mg -crdiff
+
+<job>
+problem_id = MGCRDiffusion # problem ID: basename of output filenames
+
+<nooutput1>
+file_type  = hst        # History data dump
+dt         = 0.01       # time increment between outputs
+
+<output1>
+file_type  = vtk        # Binary data dump
+variable   = prim       # variables to be output
+dcycle        = 1       # time increment between outputs
+
+<time>
+cfl_number = 0.3        # 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
+
+<mesh>
+nx1        = 64         # Number of zones in X1-direction
+x1min      = -0.5       # minimum value of X1
+x1max      = 0.5        # maximum value of X1
+ix1_bc     = outflow    # inner-X1 boundary flag
+ox1_bc     = outflow    # outer-X1 boundary flag
+
+nx2        = 64         # Number of zones in X2-direction
+x2min      = -0.5       # minimum value of X2
+x2max      = 0.5        # maximum value of X2
+ix2_bc     = outflow    # inner-X2 boundary flag
+ox2_bc     = outflow    # outer-X2 boundary flag
+
+nx3        = 64         # Number of zones in X3-direction
+x3min      = -0.5       # minimum value of X3
+x3max      = 0.5        # maximum value of X3
+ix3_bc     = outflow    # inner-X3 boundary flag
+ox3_bc     = outflow    # outer-X3 boundary flag
+
+<meshblock>
+nx1        = 64
+nx2        = 64
+nx3        = 64
+
+<hydro>
+gamma           = 1.666666666667 # gamma = C_p/C_v
+iso_sound_speed = 0.4082482905   # equivalent to sqrt(gamma*p/d) for p=0.1, d=1
+
+<crdiffusion>
+mgmode          = FMG
+fas             = true
+npresmooth      = 2
+npostsmooth     = 2
+omega           = 1.3
+niteration      = 10
+#threshold       = 0.00001
+output_defect   = true
+ix1_bc          = user 
+ox1_bc          = user
+ix2_bc          = user
+ox2_bc          = user
+ix3_bc          = zerograd
+ox3_bc          = zerograd
+
+Dpara           = 100.0
+Dperp           = 1.0
+Lambda          = 100.0
+zeta_factor     = 1.0
+
+<problem>

src/athena.hpp

@@ -50,7 +50,7 @@ class Coordinates;
 class ParameterInput;
 class HydroDiffusion;
 class FieldDiffusion;
-struct MGCoordinates;
+class MGCoordinates;
 class OrbitalAdvection;
 class NRRadiation;
 class IMRadiation;
@@ -166,8 +166,8 @@ enum CoordinateDirection {X1DIR=0, X2DIR=1, X3DIR=2};
 // KGF: Except for the 2x MG* enums, these may be unnessary w/ the new class inheritance
 // Now, only passed to BoundaryVariable::InitBoundaryData(); could replace w/ bool switch
 // TODO(tomo-ono): consider necessity of orbita_cc and orbital_fc
-enum class BoundaryQuantity {cc, fc, cc_flcor, fc_flcor, mggrav,
-                             mggrav_f, orbital_cc, orbital_fc};
+enum class BoundaryQuantity {cc, fc, cc_flcor, fc_flcor, mg, mg_faceonly, mg_coeff,
+                             orbital_cc, orbital_fc};
 enum class HydroBoundaryQuantity {cons, prim};
 enum class BoundaryCommSubset {mesh_init, gr_amr, all, orbital, radiation, radhydro};
 // TODO(felker): consider generalizing/renaming to QuantityFormulation
@@ -213,7 +213,7 @@ using FieldDiffusionCoeffFunc = void (*)(
     const AthenaArray<Real> &w,
     const AthenaArray<Real> &bmag,
     int is, int ie, int js, int je, int ks, int ke);
-using MGSourceMaskFunc = void (*)(AthenaArray<Real> &src,
+using MGMaskFunc = void (*)(AthenaArray<Real> &dat,
     int is, int ie, int js, int je, int ks, int ke, const MGCoordinates &coord);
 using OrbitalVelocityFunc = Real (*)(
     OrbitalAdvection *porb, Real x1, Real x2, Real x3);

src/bvals/bvals.hpp

@@ -124,8 +124,8 @@ class BoundaryValues : public BoundaryBase, //public BoundaryPhysics,
   std::vector<BoundaryVariable *> bvars_main_int;
   //! subset of bvars that are exchanged in the SuperTimeStepTaskList
   std::vector<BoundaryVariable *> bvars_sts;
-  //! Pointer to the Gravity Boundary Variable
-  CellCenteredBoundaryVariable *pgbvar;
+  //! Pointer to the Gravity and CRDiffusion Boundary Variable
+  CellCenteredBoundaryVariable *pgbvar, *pcrbvar;
 
   // inherited functions (interface shared with BoundaryVariable objects):
   // ------
@@ -150,8 +150,8 @@ class BoundaryValues : public BoundaryBase, //public BoundaryPhysics,
   void ProlongateBoundaries(const Real time, const Real dt,
                             std::vector<BoundaryVariable *> bvars_subset);
 
-  // temporary workaround for self-gravity
-  void ProlongateGravityBoundaries(const Real time, const Real dt);
+  // temporary workaround for Multigrid
+  void ProlongateBoundariesPostMG(CellCenteredBoundaryVariable* pbvar);
 
   // compute the shear at each integrator stage
   //! \todo (felker):
@@ -212,10 +212,11 @@ class BoundaryValues : public BoundaryBase, //public BoundaryPhysics,
   void ProlongateGhostCells(const NeighborBlock& nb,
                             int si, int ei, int sj, int ej, int sk, int ek);
 
-  // temporary workaround for self-gravity
-  void RestrictGravityGhostCellsOnSameLevel(const NeighborBlock& nb,
-                                            int nk, int nj, int ni);
-  void ProlongateGravityGhostCells(int si, int ei, int sj, int ej, int sk, int ek);
+  // temporary workaround for Multigrid
+  void RestrictGhostCellsOnSameLevelPostMG(CellCenteredBoundaryVariable* pbvar,
+                                    const NeighborBlock& nb, int nk, int nj, int ni);
+  void ProlongateGhostCellsPostMG(CellCenteredBoundaryVariable* pbvar,
+                                  int si, int ei, int sj, int ej, int sk, int ek);
 
   void DispatchBoundaryFunctions(
       MeshBlock *pmb, Coordinates *pco, Real time, Real dt,

src/bvals/bvals_refine_gravity.cpp → src/bvals/bvals_refine_postmg.cpp

@@ -24,13 +24,14 @@
 #include "bvals.hpp"
 
 
-class CellCenteredBoundaryVariable;
+class CellCenteredCellCenteredBoundaryVariable;
 
 //----------------------------------------------------------------------------------------
-//! \fn void BoundaryValues::ProlongateGravityBoundaries(const Real time, const Real dt)
-//! \brief Prolongate boundaries
+//! \fn void BoundaryValues::ProlongateBoundariesPostMG(
+//!                                              CellCenteredBoundaryVariable* pbvar)
+//! \brief Prolongate boundaries after Multigrid assuming physical boundaries are filled
 
-void BoundaryValues::ProlongateGravityBoundaries(const Real time, const Real dt) {
+void BoundaryValues::ProlongateBoundariesPostMG(CellCenteredBoundaryVariable* pbvar) {
   MeshBlock *pmb = pmy_block_;
   int &mylevel = loc.level;
 
@@ -69,7 +70,7 @@ void BoundaryValues::ProlongateGravityBoundaries(const Real time, const Real dt)
 
           // this neighbor block is on the same level
           // and needs to be restricted for prolongation
-          RestrictGravityGhostCellsOnSameLevel(nb, nk, nj, ni);
+          RestrictGhostCellsOnSameLevelPostMG(pbvar, nb, nk, nj, ni);
         }
       }
     }
@@ -108,18 +109,20 @@ void BoundaryValues::ProlongateGravityBoundaries(const Real time, const Real dt)
 
 
     // Step 3. Finally, the ghost-ghost zones are ready for prolongation:
-    ProlongateGravityGhostCells(si, ei, sj, ej, sk, ek);
+    ProlongateGhostCellsPostMG(pbvar, si, ei, sj, ej, sk, ek);
   } // end loop over nneighbor
   return;
 }
 
 
 //----------------------------------------------------------------------------------------
-//! \fn void BoundaryValues::RestrictGravityGhostCellsOnSameLevel(const NeighborBlock& nb,
-//!                                                        int nk, int nj, int ni)
-//! \brief Restrict ghost cells on same level
-void BoundaryValues::RestrictGravityGhostCellsOnSameLevel(const NeighborBlock& nb,
-                                                          int nk, int nj, int ni) {
+//! \fn void BoundaryValues::RestrictGhostCellsOnSameLevelPostMG(
+//!                          CellCenteredBoundaryVariable* pbvar, const NeighborBlock& nb,
+//!                          int nk, int nj, int ni)
+//! \brief Restrict ghost cells on same level after Multigrid
+void BoundaryValues::RestrictGhostCellsOnSameLevelPostMG(
+                             CellCenteredBoundaryVariable* pbvar, const NeighborBlock& nb,
+                             int nk, int nj, int ni) {
   MeshBlock *pmb = pmy_block_;
   MeshRefinement *pmr = pmb->pmr;
 
@@ -152,22 +155,23 @@ void BoundaryValues::RestrictGravityGhostCellsOnSameLevel(const NeighborBlock& n
     rks = pmb->cks - 1, rke = pmb->cks - 1;
   }
 
-  pmb->pmr->RestrictCellCenteredValues(*(pgbvar->var_cc), *(pgbvar->coarse_buf), 0, 0,
+  pmb->pmr->RestrictCellCenteredValues(*(pbvar->var_cc), *(pbvar->coarse_buf), 0, 0,
                                        ris, rie, rjs, rje, rks, rke);
 
   return;
 }
 
 
 //----------------------------------------------------------------------------------------
-//! \fn void BoundaryValues::ProlongateGravityGhostCells(int si, int ei, int sj, int ej,
-//!                                                      int sk, int ek)
-//! \brief Prolongate gravity ghost cells
+//! \fn void BoundaryValues::ProlongateGhostCellsPostMG(
+//!                                    CellCenteredBoundaryVariable* pbvar,
+//!                                    int si, int ei, int sj, int ej, int sk, int ek)
+//! \brief Prolongate ghost cells after Multigrid
 
-void BoundaryValues::ProlongateGravityGhostCells(int si, int ei, int sj, int ej,
-                                                 int sk, int ek) {
+void BoundaryValues::ProlongateGhostCellsPostMG(CellCenteredBoundaryVariable* pbvar,
+                               int si, int ei, int sj, int ej, int sk, int ek) {
   MeshRefinement *pmr = pmy_block_->pmr;
-  pmr->ProlongateCellCenteredValues(*(pgbvar->coarse_buf), *(pgbvar->var_cc), 0, 0,
+  pmr->ProlongateCellCenteredValues(*(pbvar->coarse_buf), *(pbvar->var_cc), 0, 0,
                                     si, ei, sj, ej, sk, ek);
   return;
 }