Merge pull request #159 from AetherModel/run_in_1d

Run in 1d

Author: aaronjridley

include/inputs.h

@@ -156,9 +156,6 @@ class Inputs {
      \param iOutput int specifying which output file type to report on
    **/
   std::string get_type_output(int iOutput);
-  
-
-
 
   /**********************************************************************
      \brief returns settings["Euv"]["dt"]
@@ -488,6 +485,7 @@ class Inputs {
    **/
   std::string get_advection_neutrals_vertical();
   bool get_advection_neutrals_bulkwinds();
+  bool get_advection_neutrals_implicitfriction();
 
 
   /**********************************************************************

include/neutrals.h

@@ -60,8 +60,11 @@ class Neutrals {
     std::vector<arma_cube> velocity_vcgc;
     std::vector<arma_cube> newVelocity_vcgc;
 
-    /// Acceleration of each species (m/s^2)
+    /// Acceleration of each species (m/s^2) due to friction term
     std::vector<arma_cube> acc_neutral_friction;
+
+    /// Coefficient for the friction term (sum of friction coefs with others)
+    arma_cube neutral_friction_coef;
 
     /// Acceleration of each species based on Eddy contribution.
     /// Only in vertical direction.
@@ -515,7 +518,6 @@ class Neutrals {
      \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);
 
@@ -524,14 +526,29 @@ class Neutrals {
      \param grid The grid to define the neutrals on
      \param time contains information about the current time
    **/
-
   bool advect_vertical(Grid grid, Times time);
 
+  /**********************************************************************
+     \brief Calculate the neutral friction in one cell using an implicit solver
+     \param iLong the longitude index
+     \param iLat the lat index
+     \param iAlt the altitude index
+     \param dt time step
+     \param vels updated velocity, which acts as a source term for the implicit solve
+   **/
   arma_vec calc_friction_one_cell(int64_t iLong, int64_t iLat, int64_t iAlt,
-				   arma_vec &vels);
+				   precision_t dt, arma_vec &vels);
 
-  void calc_neutral_friction();  
-  
+  /**********************************************************************
+     \brief Calculate the neutral friction in all cells (calls one_cell above)
+     \param dt time step
+   **/
+  void calc_neutral_friction_implicit(precision_t dt);  
+
+  /**********************************************************************
+     \brief Calculate the neutral friction coefficients for semi-implicit solver
+   **/
+  void calc_neutral_friction_coefs();    
 };
 
 #endif  // INCLUDE_NEUTRALS_H_

share/run/UA/inputs/defaults.json

@@ -24,10 +24,11 @@
 	"Neutrals" : {
 	    "Vertical" : "rusanov",
 	    "Horizontal" : "default",
-		"useBulkWinds" : true},
+	    "useBulkWinds" : true,
+	    "useImplicitFriction" : true},
 	"Ions" : {
-	"Along" : "rusanov",
-	"Across" : "default"} },
+	    "Along" : "rusanov",
+	    "Across" : "default"} },
 
     "Student" : {
 	"name" : "",
@@ -120,7 +121,7 @@
     "AuroraFile" : "UA/inputs/aurora_earth.csv",
     "IndicesLookupFile" : "UA/inputs/indices_lookup.json",
 
-    "OnmiwebFile" : [""],
+    "OmniwebFile" : ["UA/inputs/omni_20110319.txt"],
 
     "Logfile" : {
 	"name" : ["UA/output/log_geo.txt", "UA/output/log_mag.txt"],

share/run/UA/inputs/earth.in

@@ -28,10 +28,10 @@ BC is a density that is used in the lowest boundary cell if you
 
 #NEUTRALS
 name, mass, vibration, thermal_cond, thermal_exp, advect, BC
-O, 16, 5, 5.6e-4, 0.69, 1, 5.0e17
+O, 16, 5, 5.6e-4, 0.69, 1, 7.0e17
 N, 14, 5, 5.6e-4, 0.69, 1, 2.5e11
-O2, 32, 7, 3.6e-4, 0.69, 1, 4.0e18
-N2, 28, 7, 3.6e-4, 0.69, 1, 1.7e19
+O2, 32, 7, 3.6e-4, 0.69, 1, 2.0e18
+N2, 28, 7, 3.6e-4, 0.69, 1, 1.0e19
 NO, 30, 7, 3.6e-4, 0.69, 1, 5.4e14
 He, 4, 5, 5.6e-4, 0.69, 1, 1.5e14
 N_2D, 14, 5, 5.6e-4, 0.69, 0, 2.5e9

src/add_sources.cpp

@@ -19,62 +19,63 @@ void Neutrals::add_sources(Times time, Planets planet, Grid grid) {
   precision_t dt = time.get_dt();
 
   heating_sources_total = heating_euv_scgc
-                        + heating_chemical_scgc
-                        + heating_ion_friction_scgc
-                        //+ heating_ion_heat_transfer_scgc
-                        - O_cool_scgc
-                        - NO_cool_scgc;
-  
+                          + heating_chemical_scgc
+                          + heating_ion_friction_scgc
+                          //+ heating_ion_heat_transfer_scgc
+                          - O_cool_scgc
+                          - NO_cool_scgc;
+
   // Solve the laplace equations using the source terms,
   // updating the neutral temperature:
   update_temperature(grid, time);
 
   int64_t iDir, iSpec, iSpecies;
   double tSim = time.get_simulation_time();
-  
+
   // Horizontal winds use bulk winds:
-  if (input.get_use_coriolis()) 
+  if (input.get_use_coriolis())
     acc_coriolis = coriolis(velocity_vcgc, planet.get_omega(), grid.geoLat_scgc);
 
-/*
-  // Vertical winds use species winds:
-  for (iSpec = 0; iSpec < nSpeciesAdvect; iSpec++) {
-    // Pick out the advected neutral species:
-    species_chars & advected_neutral = species[species_to_advect[iSpec]];
-
-    iDir = 2;
-    // update velocities based on acceleration:
-    // reduce neutral friction until solver is added
-    advected_neutral.velocity_vcgc[iDir] =
-      advected_neutral.velocity_vcgc[iDir] +
-      dt * (ramp * grid.cent_acc_vcgc[iDir] + 
-	    ramp * acc_coriolis[iDir] + 
-	    advected_neutral.acc_neutral_friction[iDir] / 4.0 +
-	    advected_neutral.acc_ion_drag[iDir] +
-	    advected_neutral.acc_eddy);
-  }
+  /*
+    // Vertical winds use species winds:
+    for (iSpec = 0; iSpec < nSpeciesAdvect; iSpec++) {
+      // Pick out the advected neutral species:
+      species_chars & advected_neutral = species[species_to_advect[iSpec]];
 
-  calc_mass_density();
-  // Calculate bulk vertical winds:
-  velocity_vcgc[2].zeros();
-  for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
-    if (species[iSpecies].DoAdvect) {
-      velocity_vcgc[2] = velocity_vcgc[2] + 
-        species[iSpecies].mass * species[iSpecies].density_scgc % 
-        species[iSpecies].velocity_vcgc[2] / rho_scgc;
+      iDir = 2;
+      // update velocities based on acceleration:
+      // reduce neutral friction until solver is added
+      advected_neutral.velocity_vcgc[iDir] =
+        advected_neutral.velocity_vcgc[iDir] +
+        dt * (ramp * grid.cent_acc_vcgc[iDir] +
+        ramp * acc_coriolis[iDir] +
+        advected_neutral.acc_neutral_friction[iDir] / 4.0 +
+        advected_neutral.acc_ion_drag[iDir] +
+        advected_neutral.acc_eddy);
     }
-  
-  */
+
+    calc_mass_density();
+    // Calculate bulk vertical winds:
+    velocity_vcgc[2].zeros();
+    for (iSpecies = 0; iSpecies < nSpecies; iSpecies++)
+      if (species[iSpecies].DoAdvect) {
+        velocity_vcgc[2] = velocity_vcgc[2] +
+          species[iSpecies].mass * species[iSpecies].density_scgc %
+          species[iSpecies].velocity_vcgc[2] / rho_scgc;
+      }
+
+    */
 
   // Add Velocity sources to bulk winds:
   for (iDir = 0; iDir < 2; iDir++) {
     velocity_vcgc[iDir] =
       velocity_vcgc[iDir] + dt * (
-				  grid.cent_acc_vcgc[iDir] +
-				  acc_coriolis[iDir] + 
-				  acc_ion_collisions[iDir]);
+        grid.cent_acc_vcgc[iDir] +
+        acc_coriolis[iDir] +
+        acc_ion_collisions[iDir]);
     acc_sources_total[iDir].zeros();
   }
+
   // Apply Viscosity:
   update_horizontal_velocity(grid, time);
 
@@ -83,17 +84,17 @@ void Neutrals::add_sources(Times time, Planets planet, Grid grid) {
     for (iDir = 0; iDir < 2; iDir++)
       species[iSpecies].velocity_vcgc[iDir] = velocity_vcgc[iDir];
 
-  /*  
+  /*
   // If we only consider the bulk winds in the horizontal direction:
   if (input.get_advection_neutrals_bulkwinds()) {
     // Calculate Coriolis:
-    if (input.get_use_coriolis()) 
+    if (input.get_use_coriolis())
       acc_coriolis = coriolis(velocity_vcgc, planet.get_omega(), grid.geoLat_scgc);
     // Add Velocity sources to bulk winds:
     for (int iDir = 0; iDir < 3; iDir++) {
       velocity_vcgc[iDir] = velocity_vcgc[iDir] + dt * (
         grid.cent_acc_vcgc[iDir] +
-        acc_coriolis[iDir] + 
+        acc_coriolis[iDir] +
         acc_ion_collisions[iDir]);
       acc_sources_total[iDir].zeros();
     }
@@ -104,18 +105,18 @@ void Neutrals::add_sources(Times time, Planets planet, Grid grid) {
       // Pick out the advected neutral species:
       species_chars & advected_neutral = species[species_to_advect[iSpec]];
       // Calculate Coriolis:
-      if (input.get_use_coriolis()) 
-        acc_coriolis = coriolis(advected_neutral.velocity_vcgc, 
-                                planet.get_omega(), 
+      if (input.get_use_coriolis())
+        acc_coriolis = coriolis(advected_neutral.velocity_vcgc,
+                                planet.get_omega(),
                                 grid.geoLat_scgc);
 
       for (int iDir = 0; iDir < 2; iDir++) {
         // update velocities based on acceleration:
         // reduce neutral friction until solver is added
         advected_neutral.velocity_vcgc[iDir] =
           advected_neutral.velocity_vcgc[iDir] +
-          dt * (grid.cent_acc_vcgc[iDir] + 
-                acc_coriolis[iDir] + 
+          dt * (grid.cent_acc_vcgc[iDir] +
+                acc_coriolis[iDir] +
                 advected_neutral.acc_neutral_friction[iDir] / 4.0 +
                 advected_neutral.acc_ion_drag[iDir]);
         // eddy acceleration is only in the vertical direction:

src/advance.cpp

@@ -56,7 +56,7 @@ bool advance(Planets &planet,
   neutrals.calc_bulk_velocity();
   neutrals.calc_kappa_eddy();
   neutrals.calc_viscosity();
-  //neutrals.calc_cMax();
+  neutrals.calc_cMax();
 
   ions.fill_electrons();
   ions.calc_sound_speed();
@@ -65,7 +65,7 @@ bool advance(Planets &planet,
   precision_t dtNeutral = calc_dt(gGrid, neutrals.cMax_vcgc);
   precision_t dtIon = calc_dt(gGrid, ions.cMax_vcgc);
   time.calc_dt(dtNeutral, dtIon);
-  
+
   neutralsMag.calc_mass_density();
   neutralsMag.calc_mean_major_mass();
   neutralsMag.calc_specific_heat();
@@ -85,21 +85,23 @@ bool advance(Planets &planet,
   neutralsMag.calc_scale_height(mGrid);
 
   if (didWork)
-    didWork = neutrals.set_bcs(gGrid, time, indices); 
+    didWork = neutrals.set_bcs(gGrid, time, indices);
+
   if (didWork)
     didWork = ions.set_bcs(gGrid, time, indices);
+
   if (didWork)
     didWork = neutralsMag.set_bcs(mGrid, time, indices);
 
   if (gGrid.get_nAlts(false) > 1)
     neutrals.advect_vertical(gGrid, time);
 
   neutrals.exchange_old(gGrid);
-  //advect(gGrid, time, neutrals);
+  advect(gGrid, time, neutrals);
 
   if (didWork & input.get_check_for_nans())
     didWork = neutrals.check_for_nonfinites("After Horizontal Advection");
-  
+
   // ------------------------------------
   // Calculate source terms next:
 
@@ -111,6 +113,7 @@ bool advance(Planets &planet,
                        neutrals,
                        ions,
                        indices);
+
   if (didWork)
     didWork = calc_euv(planet,
                        mGrid,
@@ -126,6 +129,7 @@ bool advance(Planets &planet,
                                      time,
                                      indices,
                                      ions);
+
   if (didWork)
     didWork = electrodynamics.update(planet,
                                      mGrid,
@@ -154,8 +158,10 @@ bool advance(Planets &planet,
     calc_ion_collisions(neutrals, ions);
 
     neutrals.add_sources(time, planet, gGrid);
+
     if (didWork & input.get_check_for_nans())
       didWork = neutrals.check_for_nonfinites("After Add Sources");
+
     neutralsMag.add_sources(time, planet, mGrid);
 
     ions.calc_ion_temperature(neutrals, gGrid, time);
@@ -180,11 +186,13 @@ bool advance(Planets &planet,
 
   if (didWork)
     didWork = output(neutrals, ions, gGrid, time, planet);
+
   if (didWork)
     didWork = output(neutralsMag, ionsMag, mGrid, time, planet);
 
   if (didWork)
     didWork = logfile.write_logfile(indices, neutrals, ions, gGrid, time);
+
   if (didWork)
     didWork = logfileMag.write_logfile(indices, neutralsMag, ionsMag, mGrid, time);
 

src/aurora.cpp

@@ -144,7 +144,7 @@ void calc_aurora(Grid grid,
   };
 
   static std::vector<std::vector<precision_t>> CiArray;
-  static bool IsFirstTime = 1;
+  static bool IsFirstTime = true;
 
   // ENERGY BINS AND DE (E in eV)
   static precision_t min = 100;
@@ -158,12 +158,11 @@ void calc_aurora(Grid grid,
 
   if (!neutrals.auroraInitialized) {
     // Initialize the aurora using the auroral csv file
+    report.print(1, "Reading aurora file");
     read_aurora(neutrals, ions);
   }
 
   if (IsFirstTime) {
-    // Initialize the aurora using the auroral csv file
-    //read_aurora(neutrals, ions);
 
     precision_t lnE;
 
@@ -192,12 +191,10 @@ void calc_aurora(Grid grid,
 
       CiArray.push_back(Ci);
     }
-
-    IsFirstTime = 0;
+    IsFirstTime = false;
   }
 
-  if (report.test_verbose(4))
-    std::cout << "aurora - done with init!\n";
+  report.print(4, "aurora - done with init!");
 
   arma_vec rhoH1d;
   arma_cube scale_height;
@@ -225,8 +222,7 @@ void calc_aurora(Grid grid,
   arma_vec diff_energy_flux;
   bool DoDebug = false;
 
-  if (report.test_verbose(4))
-    std::cout << "aurora - starting main loop!\n";
+  report.print(4, "aurora - starting main loop!");
 
   // loop through each altitude and calculate ionization
   for (iLon = 0; iLon < nLons ; iLon++) {