Merge pull request #132 from ibutsky/cosmic_rays

Anisotropic cosmic rays

doc/manual/source/parameters/index.rst

@@ -3393,18 +3393,39 @@ For details on the cosmic ray solver in Enzo see :ref:`cosmic_rays`.
     Switches on diffusion of the cosmic ray energy density. Default: 0
 
     0. Off
-    1. On with constant coefficient (``CRkappa``)
+    1. On; Isotropic diffusion with constant coefficient (``CRkappa``). Compatible with HydroMethod = 2 or 4.
+    2. On; Anisotropic diffusion around magnetic field lines with constant coefficient (``CRkappa``). Compatible with HydroMethod = 4. 
 
 
 ``CRkappa`` (external)
-    Cosmic ray diffusion coefficient in CGS units (cm^2/s), Default: 0.0. For MW-like galaxies: 1E28.
+    Cosmic ray diffusion coefficient in CGS units (cm^2/s), Default: 0.0. For MW-like galaxies: 1E28-3E29.
+
+``CRStreaming`` (external)
+    Switches on streaming of the cosmic ray energy density. Default: 0
+
+   0. Off
+   1. On; Anisotropic streaming around magnetic field lines. Compatible with HydroMethod = 4. 
+
+``CRHeating`` (external)
+    This parameter adds heating of gas from streaming of cosmic rays. Physically, this should
+    only be turned on for runs with cosmic ray streaming. However, this will work for any combination of CR
+    parameters. Only compatible with HydroMethod = 4. Default: 0. 
+
+``CRStreamingStabilityFactor`` (external)
+    This is used to preserve stability in cosmic ray streaming according to the regularization scheme
+    described in Sharma et al. 2009. This number should be greater than 1 and calibrated for each simulation. 
+    Default: 1e2
+
+``CRStreamingVelocityFactor`` (external)
+    This sets the CR streaming velocity in units of the alfven velocity. Default: 1.0; 
 
 ``CRCourantSafetyNumber`` (external)
     Multiplies CR diffusion timestep, for stability should be <= 0.5. Default: 0.5
 
 ``CRFeedback`` (external)
     Specify fraction of star formation feedback energy should be diverted into the cosmic
-    ray energy density. implemented ONLY for star_maker3 (feedback method 2). Default: 0.0
+    ray energy density. Implemented ONLY for star_maker3 (feedback method 2) and star_maker2 
+    (feedback method 1). Default: 0.0
 
 ``CRdensFloor`` (external)
     Floor in gas density, can be imposed, for speed purposes (default 0.0 = off). Any value

doc/manual/source/physics/cosmic_rays.rst

@@ -6,9 +6,12 @@ Cosmic Ray Two-Fluid Model
 
 This section documents the two-fluid cosmic ray model implemented in
 Enzo, which was first used (and is described in detail) in `Salem &
-Bryan 2014 <http://adsabs.harvard.edu/abs/2014MNRAS.437.3312S>`_ .
+Bryan 2014 <http://adsabs.harvard.edu/abs/2014MNRAS.437.3312S>`_ . 
+The implementation of anisotropic cosmic ray physics is described and tested
+in `Butsky and Quinn 2018 <https://ui.adsabs.harvard.edu/abs/2018ApJ...868..108B/abstract>`_.
 For relevant parameters, please also see
-:ref:`cosmic_ray_two_fluid_model_parameters`.  The bulk of the code itself can be found in *Grid_ZeusSolver.C*
+:ref:`cosmic_ray_two_fluid_model_parameters`.  The bulk of the advection code itself can be found in *Grid_ZeusSolver.C* and *hydro_rk/Riemann_LLF_MHD.C*. Cosmic ray transport is solved in *Grid_ComputeCRDiffusion.C*, *Grid_ComputeAnisotropicCRDiffusion.C*, and *Grid_ComputeCRStreaming.C*. 
+
 
 This module models the dynamical role of cosmic rays via a set of two-fluid hydro equations
 (see `Jun et. al. 1994
@@ -20,19 +23,21 @@ star formation. See :ref:`cosmic_ray_two_fluid_model_parameters` for information
 these options. But most important are:
 
 
-  - ``CRModel`` - Switches on the CR physics (0 = off, 1 = on)
+  - ``CRModel`` - Switches on the CR physics (0 = off, 1 = on). Compatible with HydroMethod = 2 or 4.
 
   - ``CRgamma`` - For polytropic equation of state. 4/3 = relativistic, adiabatic gas (default)
 
-  - ``CRDiffusion`` - turns on diffusion of CREnergyDensity field
+  - ``CRDiffusion`` - Turns on diffusion of CREnergyDensity field (0 = off, 1 = isotropic, 2 = anisotropic). CRDiffusion = 2 is only Compatible with HydroMethod = 4. 
 
   - ``CRkappa`` - Diffusion coefficient (currently constant, isotropic)
 
-  - ``CRFeedback`` - Controls production of rays in star forming regions
+  - ``CRStreaming`` - Turns on CR streaming along mangetic field lines (0 = off, 1 = on). Compatible with HydroMethod = 4. 
+
+  - ``CRFeedback`` - Controls production of rays in star forming regions (0 = off, 1 = on) 
 
 
-For this model to run properly you *must be running the Zeus Hydro 
-Solver!* (``HydroMethod = 2``). The model has not yet been implemented for
+For this model to run properly you *must be running either the Zeus Hydro 
+Solver!* (``HydroMethod = 2``) or the Dedner MHD Solver (``HydroMethod = 4``). The model has not yet been implemented for
 any of the other fluid solvers in Enzo.
 
 If you plan on including cosmic rays, definitely first verify the solver is working by running

src/enzo/CRShockTubesInitialize.C

@@ -35,6 +35,10 @@ int CRShockTubesInitialize(FILE *fptr, FILE *Outfptr,
   char *Vel1Name = "x-velocity";
   char *Vel2Name = "y-velocity";
   char *Vel3Name = "z-velocity";
+  char *BxName = "Bx";
+  char *ByName = "By";
+  char *BzName = "Bz";
+  char *PhiName = "Phi";
   char *CRName = "CREnergyDensity";
   char *ColourName = "colour";
 
@@ -52,7 +56,10 @@ int CRShockTubesInitialize(FILE *fptr, FILE *Outfptr,
     LeftVelocityY = 0.0, RightVelocityY = 0.0, CenterVelocityY = 0.0,
     LeftVelocityZ = 0.0, RightVelocityZ = 0.0, CenterVelocityZ = 0.0,
     LeftPressure = 1.0, RightPressure = 1.0, CenterPressure = 1.0,
-    LeftCRDensity = 1.0, RightCRDensity = 1.0, CenterCRDensity = 1.0;
+    LeftCRDensity = 1.0, RightCRDensity = 1.0, CenterCRDensity = 1.0, 
+    LeftBx = 0.0, RightBx = 0.0, CenterBx = 0.0, 
+    LeftBy = 0.0, RightBy = 0.0, CenterBy = 0.0, 
+    LeftBz = 0.0, RightBz = 0.0, CenterBz = 0.0;
 
   /* read input from file */
 
@@ -104,6 +111,24 @@ int CRShockTubesInitialize(FILE *fptr, FILE *Outfptr,
       &CenterDensity);
     ret += sscanf(line, "HydroShockTubesCenterCREnDensity = %"FSYM,
       &CenterCRDensity);
+    ret += sscanf(line, "HydroShockTubesLeftBx = %"FSYM,
+                  &LeftBx);
+    ret += sscanf(line, "HydroShockTubesRightBx = %"FSYM,
+                  &RightBx);
+    ret += sscanf(line, "HydroShockTubesCenterBx = %"FSYM,
+                  &CenterBx);
+    ret += sscanf(line, "HydroShockTubesLeftBy = %"FSYM,
+                  &LeftBy);
+    ret += sscanf(line, "HydroShockTubesRightBy = %"FSYM,
+                  &RightBy);
+    ret += sscanf(line, "HydroShockTubesCenterBy = %"FSYM,
+                  &CenterBy);
+    ret += sscanf(line, "HydroShockTubesLeftBz = %"FSYM,
+                  &LeftBz);
+    ret += sscanf(line, "HydroShockTubesRightBz = %"FSYM,
+                  &RightBz);
+    ret += sscanf(line, "HydroShockTubesCenterBz = %"FSYM,
+                  &CenterBz);
 
     /* if the line is suspicious, issue a warning */
 
@@ -131,7 +156,10 @@ int CRShockTubesInitialize(FILE *fptr, FILE *Outfptr,
 				  LeftVelocityY,  RightVelocityY, CenterVelocityY,
 				  LeftVelocityZ,  RightVelocityZ, CenterVelocityZ,
 				  LeftPressure,   RightPressure,  CenterPressure,
-				  LeftCRDensity,  RightCRDensity, CenterCRDensity);
+			       LeftCRDensity,  RightCRDensity, CenterCRDensity, 
+			       LeftBx, RightBx, CenterBx, 
+			       LeftBy, RightBy, CenterBy, 
+			       LeftBz, RightBz, CenterBz);
 
   /* Convert minimum initial overdensity for refinement to mass
      (unless MinimumMass itself was actually set). */
@@ -175,7 +203,10 @@ int CRShockTubesInitialize(FILE *fptr, FILE *Outfptr,
 	   LeftVelocityY,  RightVelocityY, CenterVelocityY,
 	   LeftVelocityZ,  RightVelocityZ, CenterVelocityZ,
 	   LeftPressure,   RightPressure,  CenterPressure,
-	   LeftCRDensity,  RightCRDensity, CenterCRDensity);
+	   LeftCRDensity,  RightCRDensity, CenterCRDensity, 
+	   LeftBx, RightBx, CenterBx, 
+	   LeftBy, RightBy, CenterBy, 
+	   LeftBz, RightBz, CenterBz);
 	Temp = Temp->NextGridThisLevel;
       }
     } // end: loop over levels
@@ -209,11 +240,16 @@ int CRShockTubesInitialize(FILE *fptr, FILE *Outfptr,
   DataLabel[count++] = Vel2Name;
   DataLabel[count++] = Vel3Name;
   DataLabel[count++] = TEName;
-  DataLabel[count++] = CRName;
   if (DualEnergyFormalism) {
     DataLabel[count++] = GEName;
   }
-
+  if (HydroMethod == MHD_RK) {
+    DataLabel[count++] =  BxName;
+    DataLabel[count++] =  ByName;
+    DataLabel[count++] =  BzName;
+    DataLabel[count++] =  PhiName;
+  }
+  DataLabel[count++] = CRName;
   for (i = 0; i < count; i++)
     DataUnits[i] = NULL;
 
@@ -263,7 +299,24 @@ int CRShockTubesInitialize(FILE *fptr, FILE *Outfptr,
 	    CenterPressure);
     fprintf(Outfptr, "HydroShockTubesCenterCREnDensity     = %"FSYM"\n",
       CenterCRDensity);
-
+    fprintf(Outfptr, "HydroShockTubesLeftBx     = %"FSYM"\n",
+	    LeftBx);
+    fprintf(Outfptr, "HydroShockTubesRightBx     = %"FSYM"\n",
+            RightBx);
+    fprintf(Outfptr, "HydroShockTubesCenterBx     = %"FSYM"\n",
+            CenterBx);
+    fprintf(Outfptr, "HydroShockTubesLeftBy     = %"FSYM"\n",
+            LeftBy);
+    fprintf(Outfptr, "HydroShockTubesRightBy     = %"FSYM"\n",
+            RightBy);
+    fprintf(Outfptr, "HydroShockTubesCenterBy     = %"FSYM"\n",
+            CenterBy);
+    fprintf(Outfptr, "HydroShockTubesLeftBz     = %"FSYM"\n",
+            LeftBz);
+    fprintf(Outfptr, "HydroShockTubesRightBz     = %"FSYM"\n",
+            RightBz);
+    fprintf(Outfptr, "HydroShockTubesCenterBz     = %"FSYM"\n",
+            CenterBz);
   }
 
   return SUCCESS;