Merge pull request #127 from AetherModel/vertical_2nd_try

Incorporate a vertical advection scheme into Aether

include/grid.h

@@ -107,6 +107,12 @@ class Grid {
   arma_cube dalt_ratio_scgc;
   arma_cube dalt_ratio_sq_scgc;
 
+  arma_cube MeshCoefm2;
+  arma_cube MeshCoefm1;
+  arma_cube MeshCoefp0;
+  arma_cube MeshCoefp1;
+  arma_cube MeshCoefp2;
+
   arma_cube dlon_center_scgc;
   arma_cube dlon_center_dist_scgc;
 

include/inputs.h

@@ -72,7 +72,9 @@ class Inputs {
 
   json get_initial_condition_types();
   json get_boundary_condition_types();
-
+
+  std::string get_advection_neutrals_vertical();
+
   // ------------------------------
   // Grid inputs:
 

include/ions.h

@@ -90,6 +90,7 @@ class Ions {
   species_chars create_species(Grid grid);
   int read_planet_file(Planets planet);
   void init_ion_temperature(Neutrals neutrals, Grid grid);
+  void set_floor();
   void fill_electrons();
   int get_species_id(std::string name);
   void calc_efield(Grid grid);

include/neutrals.h

@@ -50,13 +50,15 @@ class Neutrals {
 
     /// Number density of species (/m3)
     arma_cube density_scgc;
+    arma_cube newDensity_scgc;
 
     /// Velocity of each species (m/s). For all below:
     /// Index 0 = longitudinal component of velocity
     /// Index 1 = latitudinal
     /// Index 2 = altitudinal
     std::vector<arma_cube> velocity_vcgc;
-
+    std::vector<arma_cube> newVelocity_vcgc;
+
     /// Acceleration of each species (m/s^2)
     std::vector<arma_cube> acc_neutral_friction;
 
@@ -135,8 +137,12 @@ class Neutrals {
   /// bulk velocity (m/s)
   std::vector<arma_cube> velocity_vcgc;
 
+  /// sound speed + abs(bulk velocity (m/s))
+  std::vector<arma_cube> cMax_vcgc;
+
   /// bunk temperature (K)
   arma_cube temperature_scgc;
+  arma_cube newTemperature_scgc;
 
   /// bulk mass density (kg/m3)
   arma_cube rho_scgc;
@@ -170,6 +176,9 @@ class Neutrals {
 
   /// Vector of all species-specific items:
   std::vector<species_chars> species;
+
+  /// when computing dt, derive a dt for neutrals:
+  precision_t dt;
 
   /// Maximum Chapman integral (will give nearly infinite tau in EUV)
   precision_t max_chapman = 1.0e26;
@@ -253,13 +262,18 @@ class Neutrals {
      temperature.  It is temporary until we get a vertical solver.
 
      \param iSpecies The species to fill (optional)
+     \param iStart The starting altitude to work with
+     \param iEnd The ending altitude to work with (NOT INCLUDED!!)
      \param grid The grid to define the neutrals on
    **/
-  void fill_with_hydrostatic(Grid grid);
+  void fill_with_hydrostatic(int64_t iStart,
+			     int64_t iEnd,
+			     Grid grid);
 
   void fill_with_hydrostatic(int64_t iSpecies,
+			     int64_t iStart,
+			     int64_t iEnd,
 			     Grid grid);
-
 
   /**********************************************************************
      \brief Calculate the bulk mass density from individual species densities
@@ -302,6 +316,18 @@ class Neutrals {
    **/
   void calc_specific_heat();
 
+  /**********************************************************************
+     \brief Calculate speed of sound + abs(velocity) in all 3 directions
+   **/
+  void calc_cMax();
+
+  /**********************************************************************
+     \brief Calculate dt (cell size / cMax) in each direction, and take min
+     \param dt returns the neutral time-step
+     \param grid The grid to define the neutrals on
+   **/
+  precision_t calc_dt(Grid grid);
+
   /**********************************************************************
      \brief Calculate the chapman integrals for the individual species
      \param grid The grid to define the neutrals on
@@ -433,6 +459,24 @@ class Neutrals {
 		       bool DoReverseX,
 		       bool DoReverseY,
 		       bool XbecomesY);
+
+  /**********************************************************************
+     \brief Vertical advection solver - Rusanov 
+     \param grid The grid to define the neutrals on
+     \param time contains information about the current time
+   **/
+
+  void solver_vertical_rusanov(Grid grid,
+			       Times time);
+
+  /**********************************************************************
+     \brief Call the correct vertical advection scheme
+     \param grid The grid to define the neutrals on
+     \param time contains information about the current time
+   **/
+
+  bool advect_vertical(Grid grid, Times time);
+
 };
 
 #endif  // INCLUDE_NEUTRALS_H_

include/solvers.h

@@ -56,6 +56,7 @@ arma_cube calc_gradient_lon(arma_cube value, Grid grid);
 arma_cube calc_gradient_lat(arma_cube value, Grid grid);
 arma_cube calc_gradient_alt(arma_cube value, Grid grid);
 std::vector<arma_cube> calc_gradient_vector(arma_cube value_scgc, Grid grid);
+arma_cube calc_gradient_alt_4th(arma_cube value, Grid grid);
 
 // interpolation in 1D
 precision_t linear_interpolation(const precision_t y0,
@@ -68,4 +69,8 @@ precision_t interpolate_unit_cube(const arma_cube &data,
                                   const precision_t yRatio,
                                   const precision_t zRatio);
 
+precision_t limiter_mc(precision_t dUp,
+		       precision_t dDown,
+		       precision_t beta);
+
 #endif  // INCLUDE_SOLVERS_H_

include/times.h

@@ -135,12 +135,13 @@ class Times {
 
   /**************************************************************
      \brief Calculates the delta time in the code.
-
+     \param dtNeutral dt from the neutrals
+     \param dtIon dt from the ions
      This calculates the delta-time in Aether for the current state.
    At this moment, this is a stand-in code. In reality, we need to pass
    in the grid, neutrals, ions, and inputs for calculation.
    **/
-  void calc_dt();
+  void calc_dt(precision_t dtNeutral, precision_t dtIon);
 
   /**************************************************************
      \brief Get the current time as an array

include/tools.h

@@ -27,7 +27,19 @@ void display_vector(arma_vec vec);
 
 bool sync_across_all_procs(bool value);
 
-// ----------------------------------------------------------------------
+// ----------------------------------------------------------------------------
+// Find max across all processors and return value to everyone
+// ----------------------------------------------------------------------------
+
+precision_t sync_max_across_all_procs(precision_t value);
+
+// ----------------------------------------------------------------------------
+// Find min across all processors and return value to everyone
+// ----------------------------------------------------------------------------
+
+precision_t sync_min_across_all_procs(precision_t value);
+
+// ----------------------------------------------------------------------------
 // Calculate the average value across all processors
 // ----------------------------------------------------------------------
 

share/run/UA/inputs/defaults.json

@@ -21,6 +21,14 @@
     "BoundaryConditions" : {
 	"type" : "Planet"},
 
+    "Advection" : {
+	"Neutrals" : {
+	    "Vertical" : "rusanov",
+	    "Horizontal" : "default"},
+	"Ions" : {
+	    "Along" : "rusanov",
+	    "Across" : "default"} },
+
     "Student" : {
 	"name" : "",
 	"is" : false },
@@ -36,11 +44,6 @@
 	"HeatingEfficiency" : 0.05,
 	"dt" : 60.0},
 
-    "Sources": {
-	"Neutrals" : {
-	    "NO_cool" : false,
-	    "O_cool" : false } },
-
     "DoCalcBulkIonTemp" : false,
 
     "Eddy" : {

share/run/UA/inputs/ion_neutral_collision_frequencies.csv

@@ -19,7 +19,7 @@ CO2+,0.47,0.51,2,1.76,3.4,3.22,3.18,R,,,,,
 Ion,Neutral,Temp,Coef1,Tn Frac,Ti Frac,Coef2,,,,,,,
 H+,H,50,2.65E-16,0.5,0.5,0.083,,,,,,,
 He+,He,50,8.73E-17,0.5,0.5,0.093,,,,,,,
-N+,Neutral,275,3.83E-17,0.5,0.5,0.063,,,,,,,
+N+,N,275,3.83E-17,0.5,0.5,0.063,,,,,,,
 O+,O,235,3.67E-17,0.5,0.5,0.064,,,,,,,
 O+_2D,O,235,3.67E-17,0.5,0.5,0.064,,,,,,,
 O+_2P,O,235,3.67E-17,0.5,0.5,0.064,,,,,,,

src/add_sources.cpp

@@ -19,27 +19,29 @@ void Neutrals::add_sources(Times time) {
   precision_t dt = time.get_dt();
 
   temperature_scgc = temperature_scgc +
-    dt * (heating_euv_scgc
-          + heating_chemical_scgc
-          + conduction_scgc
-          - O_cool_scgc
-          - NO_cool_scgc);
+                     dt * (heating_euv_scgc
+                           + heating_chemical_scgc
+                           + conduction_scgc
+                           - O_cool_scgc
+                           - NO_cool_scgc);
 
   for (int64_t iSpec = 0; iSpec < nSpecies; iSpec++) {
     for (int iDir = 0; iDir < 3; iDir++) {
       // update velocities based on acceleration:
       // reduce neutral friction until solver is added
       species[iSpec].velocity_vcgc[iDir] =
-	species[iSpec].velocity_vcgc[iDir] +
-	dt * (species[iSpec].acc_neutral_friction[iDir]/4.0 +
-	      species[iSpec].acc_ion_drag[iDir]);
+        species[iSpec].velocity_vcgc[iDir] +
+        dt * (species[iSpec].acc_neutral_friction[iDir] / 4.0 +
+              species[iSpec].acc_ion_drag[iDir]);
+
       // eddy acceleration is only in the vertical direction:
       if (iDir == 2)
-	species[iSpec].velocity_vcgc[iDir] =
-	  species[iSpec].velocity_vcgc[iDir] =
-	  dt * species[iSpec].acc_eddy;
+        species[iSpec].velocity_vcgc[iDir] =
+          species[iSpec].velocity_vcgc[iDir] +
+          dt * species[iSpec].acc_eddy;
     }
   }
+
   calc_bulk_velocity();
 
   report.exit(function);

src/advance.cpp

@@ -43,7 +43,24 @@ int advance(Planets &planet,
   neutrals.calc_pressure();
   neutrals.calc_bulk_velocity();
   neutrals.calc_kappa_eddy();
-  time.calc_dt();
+
+  neutrals.calc_cMax();
+  precision_t dtNeutral = neutrals.calc_dt(gGrid);
+  precision_t dtIon = 100.0;
+  time.calc_dt(dtNeutral, dtIon);
+
+  // ------------------------------------
+  // Do advection first :
+
+  // Upper BCs requires the scale height to be calculated, so do that
+  // first
+
+  neutrals.calc_scale_height(gGrid);
+  neutrals.set_bcs(gGrid, time, indices);
+  neutrals.advect_vertical(gGrid, time);
+
+  // ------------------------------------
+  // Calculate source terms next:
 
   iErr = calc_euv(planet,
                   gGrid,
@@ -65,11 +82,12 @@ int advance(Planets &planet,
   // Calculate some neutral source terms:
   neutrals.calc_conduction(gGrid, time);
   chemistry.calc_chemistry(neutrals, ions, time, gGrid);
+
   if (input.get_O_cooling())
     neutrals.calc_O_cool();
+
   if (input.get_NO_cooling())
     neutrals.calc_NO_cool();
-  neutrals.add_sources(time);
 
   neutrals.calc_conduction(gGrid, time);
   chemistry.calc_chemistry(neutrals, ions, time, gGrid);
@@ -79,12 +97,10 @@ int advance(Planets &planet,
   calc_neutral_friction(neutrals);
 
   neutrals.add_sources(time);
+
   ions.calc_ion_temperature(neutrals, gGrid, time);
   ions.calc_electron_temperature(neutrals, gGrid);
 
-  neutrals.set_bcs(gGrid, time, indices);
-  neutrals.calc_scale_height(gGrid);
-  neutrals.fill_with_hydrostatic(gGrid);
   neutrals.exchange(gGrid);
 
   time.increment_time();

src/calc_chemical_sources.cpp

@@ -54,16 +54,19 @@ void Chemistry::calc_chemical_sources(Neutrals &neutrals,
     //    reaction rate * loss den 1 * loss den 2 (* loss den 3 if needed)
 
     change3d.fill(rate);
+
     // Check for type of temperature dependence and calculate
     if (reactions[iReaction].type > 0) {
       // Determined which temperature to use in equation:
       // use Tn by default
       arma_cube temp = neutrals.temperature_scgc;
       std::string denom = reactions[iReaction].denominator;
+
       if (denom == "Te")
         temp = ions.electron_temperature_scgc;
       else if (denom == "Ti")
         temp = ions.temperature_scgc;
+
       // Calculate reaction rate:
       if (reactions[iReaction].numerator &&
           reactions[iReaction].type == 1) {
@@ -108,7 +111,8 @@ void Chemistry::calc_chemical_sources(Neutrals &neutrals,
 
       }
     }
-    // if temperature dependence is piecewise, only operate on cells
+
+    // If temperature dependence is piecewise, only operate on cells
     // within temperature range:
     if (reactions[iReaction].min || reactions[iReaction].max) {
       // Figure out which temperature is the limiter.  Default to ions:
@@ -126,7 +130,7 @@ void Chemistry::calc_chemical_sources(Neutrals &neutrals,
       if (reactions[iReaction].max > 0)
         change3d = change3d % (change3d <= reactions[iReaction].max);
     }
-    
+
     // Now that the reaction rate is calculated, multiply by the
     // densities on the left side of the equation (loss terms):
     for (iLoss = 0; iLoss < reactions[iReaction].nLosses; iLoss++) {

src/calc_chemistry.cpp

@@ -76,7 +76,7 @@ void Chemistry::calc_chemistry(Neutrals &neutrals,
   // ---------------------------------------------------------
   // Recalculate electrons
   // ---------------------------------------------------------
-
+  ions.set_floor();
   ions.fill_electrons();
 
   report.exit(function);