Merge pull request #145 from AetherModel/IEaddition

Add Fortran IE models such as Weimer and FTA

Author: aaronjridley

include/neutrals.h

@@ -178,9 +178,6 @@ 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;
 

include/tools.h

@@ -15,6 +15,13 @@ struct mat_2x2{
     arma_mat A22;
 };
 
+// ----------------------------------------------------------------------------
+// Fix corners in an arma cube
+//   - basically fill in the corners with values near them
+// ----------------------------------------------------------------------------
+
+void fill_corners(arma_cube &values, int64_t nGCs);
+
 // -----------------------------------------------------------------------------
 // add cMember into a string just before last period
 // -----------------------------------------------------------------------------

src/advance.cpp

@@ -35,6 +35,9 @@ bool advance(Planets &planet,
     report.print(-1, "(1) What function is this " +
                  input.get_student_name() + "?");
 
+  if (didWork & input.get_check_for_nans())
+    didWork = neutrals.check_for_nonfinites();
+
   gGrid.calc_sza(planet, time);
   neutrals.calc_mass_density();
   neutrals.calc_mean_major_mass();
@@ -56,11 +59,17 @@ bool advance(Planets &planet,
   // first
 
   neutrals.calc_scale_height(gGrid);
-  didWork = neutrals.set_bcs(gGrid, time, indices);
+
+  if (didWork)
+    didWork = neutrals.set_bcs(gGrid, time, indices);
 
   if (input.get_nAltsGeo() > 1)
     neutrals.advect_vertical(gGrid, time);
 
+  if (didWork & input.get_check_for_nans())
+    didWork = neutrals.check_for_nonfinites();
+
+
   // ------------------------------------
   // Calculate source terms next:
 
@@ -97,9 +106,6 @@ bool advance(Planets &planet,
     if (input.get_NO_cooling())
       neutrals.calc_NO_cool();
 
-    neutrals.calc_conduction(gGrid, time);
-    chemistry.calc_chemistry(neutrals, ions, time, gGrid);
-
     neutrals.vertical_momentum_eddy(gGrid);
     calc_ion_collisions(neutrals, ions);
     calc_neutral_friction(neutrals);
@@ -124,6 +130,9 @@ bool advance(Planets &planet,
     }
   }
 
+  if (didWork & input.get_check_for_nans())
+    didWork = neutrals.check_for_nonfinites();
+
   if (didWork)
     didWork = output(neutrals, ions, gGrid, time, planet);
 

src/calc_chemical_sources.cpp

@@ -39,7 +39,7 @@ void Chemistry::calc_chemical_sources(Neutrals &neutrals,
 
   for (iReaction = 0; iReaction < nReactions; iReaction++) {
 
-    if (report.test_verbose(3)) {
+    if (report.test_verbose(4)) {
       std::cout << "===> Reaction : " << iReaction
                 << " of " << nReactions << "\n";
       display_reaction(reactions[iReaction]);

src/calc_neutral_derived.cpp

@@ -337,6 +337,8 @@ precision_t Neutrals::calc_dt(Grid grid) {
   static int iFunction = -1;
   report.enter(function, iFunction);
 
+  precision_t dt;
+
   if (input.get_is_cubesphere())
     dt = calc_dt_cubesphere(grid);
   else {
@@ -382,6 +384,7 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) {
   static int iFunction = -1;
   report.enter(function, iFunction);
 
+  precision_t dt;
   int iDir;
 
   arma_vec dta(4);
@@ -390,7 +393,7 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) {
   int64_t nAlts = grid.get_nAlts();
   int64_t nXs = grid.get_nLons();
   int64_t nYs = grid.get_nLats();
-  
+
   // dtx dty for reference coordinate system
   arma_cube dtx(size(cMax_vcgc[0]));
   arma_cube dty(size(cMax_vcgc[0]));
@@ -401,11 +404,18 @@ precision_t Neutrals::calc_dt_cubesphere(Grid grid) {
   // Loop through altitudes
   for (int iAlt = 0; iAlt < nAlts; iAlt++) {
     // Conver cMax to contravariant velocity first
-    arma_mat u1 = cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) + cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt);
-    arma_mat u2 = cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) + cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt);
-
-    dtx.slice(iAlt) = grid.drefx(iAlt)*dummy_1 / u1; 
-    dty.slice(iAlt) = grid.drefy(iAlt)*dummy_1 / u2; 
+    arma_mat u1 = sqrt(
+                    cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) %
+                    cMax_vcgc[0].slice(iAlt) % grid.A11_inv_scgc.slice(iAlt) +
+                    cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt) %
+                    cMax_vcgc[1].slice(iAlt) % grid.A12_inv_scgc.slice(iAlt));
+    arma_mat u2 = sqrt(
+                    cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) %
+                    cMax_vcgc[0].slice(iAlt) % grid.A21_inv_scgc.slice(iAlt) +
+                    cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt) %
+                    cMax_vcgc[1].slice(iAlt) % grid.A22_inv_scgc.slice(iAlt));
+    dtx.slice(iAlt) = grid.drefx(iAlt) * dummy_1 / u1;
+    dty.slice(iAlt) = grid.drefy(iAlt) * dummy_1 / u2;
   }
 
   // simply some things, and just take the bulk value for now:

src/electrodynamics.cpp

@@ -5,7 +5,7 @@
 
 
 // -----------------------------------------------------------------------------
-// Call (fortran) ionospheric electrodynamics if we enable this 
+// Call (fortran) ionospheric electrodynamics if we enable this
 // in the cmake file
 // -----------------------------------------------------------------------------
 
@@ -27,36 +27,40 @@ Electrodynamics::Electrodynamics(Times time) {
   HaveFortranIe = false;
   isSet = false;
 
-  #ifdef FORTRAN
-    // Initialize IE components (reading in the data files):
-    std::string ieDir = input.get_electrodynamics_dir();
-    int* ieDir_iArray = copy_string_to_int(ieDir);
-    std::string efield = input.get_potential_model();
-    std::string aurora = input.get_diffuse_auroral_model();
+#ifdef FORTRAN
+  // Initialize IE components (reading in the data files):
+  std::string ieDir = input.get_electrodynamics_dir();
+  int* ieDir_iArray = copy_string_to_int(ieDir);
+  std::string efield = input.get_potential_model();
+  std::string aurora = input.get_diffuse_auroral_model();
 
-    if (efield.length() == 0 & aurora.length() == 0) {
-      HaveFortranIe = false;
+  if (efield.length() == 0 & aurora.length() == 0)
+    HaveFortranIe = false;
+
+  else {
+    int* efield_iArray = copy_string_to_int(efield);
+    int* aurora_iArray = copy_string_to_int(aurora);
+    ie_init_library(ieDir_iArray, efield_iArray, aurora_iArray, &iError);
+
+    if (iError > 0) {
+      report.print(0, "Error in setting fortran IE!");
+      IsOk = false;
     } else {
-      int* efield_iArray = copy_string_to_int(efield);
-      int* aurora_iArray = copy_string_to_int(aurora);
-      ie_init_library(ieDir_iArray, efield_iArray, aurora_iArray, &iError);
-      if (iError > 0) {
-        report.print(0,"Error in setting fortran IE!");
-        IsOk = false;
-      } else {
-        HaveFortranIe = true;
-        isSet = true;
-      }
+      HaveFortranIe = true;
+      isSet = true;
     }
-    isCompiled = true;
-  #else
-    isCompiled = false;
-  #endif
+  }
+
+  isCompiled = true;
+#else
+  isCompiled = false;
+#endif
 
   // If we don't set IE through Fortran, then try to set it through a file:
 
   if (!HaveFortranIe) {
     std::string electrodynamics_file = input.get_electrodynamics_file();
+
     if (electrodynamics_file.length() > 0 &
         electrodynamics_file != "none") {
       // This function sets HaveElectrodynamicsFile = true.
@@ -83,7 +87,7 @@ Electrodynamics::Electrodynamics(Times time) {
 // -----------------------------------------------------------------------------
 
 void Electrodynamics::set_all_indices_for_ie(Times time,
-                                              Indices &indices) {
+                                             Indices &indices) {
   std::string function = "Electrodynamics::set_all_indices_for_ie";
   static int iFunction = -1;
   report.enter(function, iFunction);
@@ -111,17 +115,18 @@ void Electrodynamics::set_all_indices_for_ie(Times time,
     std::cout << "sw v : " << iVx_ << " " << swv << "\n";
     std::cout << "sw n : " << iN_ << " " << swn << "\n";
   }
-  #ifdef FORTRAN
-    ie_set_time(&time_now);
-    ie_set_imfby(&imfby);
-    ie_set_imfbz(&imfbz);
-    ie_set_swv(&swv);
-    ie_set_swn(&swn);
-    ie_set_ae(&ae);
-    ie_set_au(&au);
-    ie_set_al(&al);
-    ie_set_hp_from_ae(&ae);
-  #endif
+
+#ifdef FORTRAN
+  ie_set_time(&time_now);
+  ie_set_imfby(&imfby);
+  ie_set_imfbz(&imfbz);
+  ie_set_swv(&swv);
+  ie_set_swn(&swn);
+  ie_set_ae(&ae);
+  ie_set_au(&au);
+  ie_set_al(&al);
+  ie_set_hp_from_ae(&ae);
+#endif
 
   report.exit(function);
   return;
@@ -154,60 +159,62 @@ bool Electrodynamics::update(Planets planet,
     ions.eflux.zeros();
     ions.avee.ones();
 
-    #ifdef FORTRAN
-      if (HaveFortranIe) {
-        report.print(3, "Using Fortran Electrodynamics!");
-        set_all_indices_for_ie(time, indices);
-
-        if (!IsAllocated) {
-          int nXs = gGrid.get_nX();
-          ie_set_nxs(&nXs);
-          int nYs = gGrid.get_nY();
-          ie_set_nys(&nYs);
-          int64_t iTotal = nXs * nYs;
-          mlt2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
-          lat2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
-          pot2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
-          eflux2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
-          avee2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
-          IsAllocated = true;
-        }
+#ifdef FORTRAN
+
+    if (HaveFortranIe) {
+      report.print(3, "Using Fortran Electrodynamics!");
+      set_all_indices_for_ie(time, indices);
+
+      if (!IsAllocated) {
+        int nXs = gGrid.get_nX();
+        ie_set_nxs(&nXs);
+        int nYs = gGrid.get_nY();
+        ie_set_nys(&nYs);
+        int64_t iTotal = nXs * nYs;
+        mlt2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
+        lat2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
+        pot2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
+        eflux2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
+        avee2d = static_cast<float*>(malloc(iTotal * sizeof(float)));
+        IsAllocated = true;
+      }
 
-        int64_t nZs = gGrid.get_nZ();
-        int64_t iZ;
+      int64_t nZs = gGrid.get_nZ();
+      int64_t iZ;
 
-        int iError;
+      int iError;
 
-        for (iZ = 0; iZ < nZs; iZ++) {
-          copy_mat_to_array(gGrid.magLocalTime_scgc.slice(iZ), mlt2d, true);
-          copy_mat_to_array(gGrid.magLat_scgc.slice(iZ), lat2d, true);
+      for (iZ = 0; iZ < nZs; iZ++) {
+        copy_mat_to_array(gGrid.magLocalTime_scgc.slice(iZ), mlt2d, true);
+        copy_mat_to_array(gGrid.magLat_scgc.slice(iZ), lat2d, true);
 
-          ie_set_mlts(mlt2d, &iError);
-          ie_set_lats(lat2d, &iError);
-          ie_update_grid(&iError);
+        ie_set_mlts(mlt2d, &iError);
+        ie_set_lats(lat2d, &iError);
+        ie_update_grid(&iError);
 
-          ie_get_potential(pot2d, &iError);
-          copy_array_to_mat(pot2d, ions.potential_scgc.slice(iZ), true);
+        ie_get_potential(pot2d, &iError);
+        copy_array_to_mat(pot2d, ions.potential_scgc.slice(iZ), true);
 
-          if (iZ == nZs-1) {
-            ie_get_electron_diffuse_aurora(eflux2d, avee2d, &iError);
-            copy_array_to_mat(avee2d, ions.avee, true);
-            copy_array_to_mat(eflux2d, ions.eflux, true);
-          }
+        if (iZ == nZs - 1) {
+          ie_get_electron_diffuse_aurora(eflux2d, avee2d, &iError);
+          copy_array_to_mat(avee2d, ions.avee, true);
+          copy_array_to_mat(eflux2d, ions.eflux, true);
         }
       }
-    #endif
+    }
+
+#endif
 
     if (HaveElectrodynamicsFile) {
       report.print(3, "Setting electrodynamics from file!");
       auto electrodynamics_values =
         get_electrodynamics(gGrid.magLat_scgc,
-                             gGrid.magLocalTime_scgc);
+                            gGrid.magLocalTime_scgc);
       ions.potential_scgc = std::get<0>(electrodynamics_values);
       ions.eflux = std::get<1>(electrodynamics_values);
       ions.avee = std::get<2>(electrodynamics_values);
     }
-  } 
+  }
 
   report.exit(function);
   return true;

src/exchange_messages_v2.cpp

@@ -601,16 +601,27 @@ void Grid::exchange(arma_cube &data, const bool pole_inverse) {
 
 bool Neutrals::exchange(Grid &grid) {
   // For each species, exchange if its DoAdvect is true
+  int64_t nGCs = grid.get_nGCs();
+
   for (int i = 0; i < nSpecies; ++i) {
     if (species[i].DoAdvect)
       grid.exchange(species[i].density_scgc, false);
+
+    fill_corners(species[i].density_scgc, nGCs);
   }
 
   // Exchange temperature
   grid.exchange(temperature_scgc, false);
+  // Fix corners:
+  fill_corners(temperature_scgc, nGCs);
+
   // Exchange velocity
   grid.exchange(velocity_vcgc[0], true);
   grid.exchange(velocity_vcgc[1], true);
   grid.exchange(velocity_vcgc[2], false);
+
+  for (int iDir = 0; iDir < 3; iDir++)
+    fill_corners(velocity_vcgc[iDir], nGCs);
+
   return true;
 }

src/init_geo_grid.cpp

@@ -984,6 +984,12 @@ bool Grid::init_geo_grid(Quadtree quadtree,
 
   // Calculate the radius (for spherical or non-spherical)
   fill_grid_radius(planet);
+
+  // Correct the reference grid with correct length scale:
+  // (with R = actual radius)
+  if (input.get_is_cubesphere())
+    correct_xy_grid(planet);
+
   // Calculate grid spacing
   calc_grid_spacing(planet);
   //calculate radial unit vector (for spherical or oblate planet)
@@ -994,11 +1000,6 @@ bool Grid::init_geo_grid(Quadtree quadtree,
   // Calculate magnetic field and magnetic coordinates:
   fill_grid_bfield(planet);
 
-  // Correct the reference grid with correct length scale:
-  // (with R = actual radius)
-  if (input.get_is_cubesphere())
-    correct_xy_grid(planet);
-
   // Throw a little message for students:
   report.student_checker_function_name(input.get_is_student(),
                                        input.get_student_name(),