Merge pull request #156 from AetherModel/advection_v2

Advection v2

Author: aaronjridley

ext/IE/EIE_Initialize.f90

@@ -168,18 +168,23 @@ subroutine EIE_Initialize(iOutputError)
       IsFound_AuroralModel = .true.
    endif
 
-  if (index(EIE_NameOfAuroralModel,'hpi') > 0) then
+   if (index(EIE_NameOfAuroralModel,'hpi') > 0) then
       call read_conductance_model(iError)
       IsFound_AuroralModel = .true.
    endif
-  if (index(EIE_NameOfAuroralModel,'pem') > 0) then
+   if (index(EIE_NameOfAuroralModel,'pem') > 0) then
       call read_conductance_model(iError)
       IsFound_AuroralModel = .true.
    endif
-  if (index(EIE_NameOfAuroralModel,'fta') > 0) then
+   if (index(EIE_NameOfAuroralModel,'fta') > 0) then
+      IsFound_AuroralModel = .true.
+   endif
+   if (index(EIE_NameOfAuroralModel,'zero') > 0) then
+      IsFound_AuroralModel = .true.
+   endif
+   if (index(EIE_NameOfAuroralModel,'none') > 0) then
       IsFound_AuroralModel = .true.
    endif
-
   if (iDebugLevel > 4) write(*,*) "=====> Back from read conductance"
 
   if (index(EIE_NameOfEFieldModel,'amie') > 0) then

ext/IE/EIE_IoLibrary.f90

@@ -604,6 +604,14 @@ subroutine IO_GetElectronDiffuseAurora(EFluxOut, AveEOut, iError)
 
       endif
 
+      if ((index(EIE_NameOfAuroralModel,'zero') > 0) .or. &
+          (index(EIE_NameOfAuroralModel,'none') > 0)) then
+         iError = 0
+         AveEOut = 3.0
+         EFluxOut = 1e-9
+      endif
+
+
   endif
 
 end subroutine IO_GetElectronDiffuseAurora

include/calc_momentum_friction.h

@@ -4,11 +4,8 @@
 #ifndef INCLUDE_CALC_MOMENTUM_FRICTION_H_
 #define INCLUDE_CALC_MOMENTUM_FRICTION_H_
 
-arma_vec neutral_friction_one_cell(int64_t iLong, int64_t iLat, int64_t iAlt,
-				   arma_vec &vels,
-				   Neutrals &neutrals);
+#include "../include/aether.h"
 
-void calc_neutral_friction(Neutrals &neutrals);
 
 /**********************************************************************
   \brief Calculate acceleration due to ion drag
@@ -17,7 +14,6 @@ void calc_neutral_friction(Neutrals &neutrals);
   \param dt The change in time
   \param report allow reporting to occur
 **/
-void calc_ion_collisions(Neutrals &neutrals,
-			 Ions &ions);
+void calc_ion_collisions(Neutrals &neutrals, Ions &ions);
 
 #endif  // INCLUDE_CALC_MOMENTUM_FRICTION_H_

include/grid.h

@@ -127,6 +127,14 @@ class Grid
   arma_cube MeshCoefp1;
   arma_cube MeshCoefp2;
 
+  // This is for a one-sided 3rd order gradient for the bottom boundary:
+
+  arma_cube MeshCoef1s3rdp1;
+  arma_cube MeshCoef1s3rdp2;
+  arma_cube MeshCoef1s3rdp3;
+  arma_cube MeshCoef1s3rdp4;
+  arma_cube MeshCoef1s3rdp5;
+
   arma_cube dlon_center_scgc;
   arma_cube dlon_center_dist_scgc;
 
@@ -137,6 +145,21 @@ class Grid
   // Vector of dx dy of different altitudes
   arma_vec drefx, drefy;
 
+  /// These are switching to the LR and DU directions for generalized coords
+  /// They are also in radians
+
+  arma_cube x_Center, y_Center;
+  arma_cube x_Left, y_Down;
+
+  /// these are center-to-center distances in the LR (X) and DU (Y) directions:
+  arma_cube dx_Center, dy_Center;
+  /// need dx on the lower / upper edges, don't need them on the left/right
+  arma_cube dx_Down;
+  /// need dy on the left / right edges:
+  arma_cube dy_Left;
+  /// cell area (in radians^2)
+  arma_cube cell_area;
+
   std::vector<arma_cube> bfield_vcgc;
   arma_cube bfield_mag_scgc;
   std::vector<arma_cube> bfield_unit_vcgc;

include/inputs.h

@@ -438,13 +438,23 @@ class Inputs {
    **/
   bool get_O_cooling();
 
+  /**********************************************************************
+     \brief returns settings["
+     \param 
+   **/
+  bool get_use_centripetal();
+  
+  /**********************************************************************
+     \brief returns settings["
+     \param 
+   **/
+  bool get_use_coriolis();
 
   /**********************************************************************
      \brief returns settings["
      \param 
    **/
   bool get_cent_acc();
-
 
   /**********************************************************************
      \brief returns settings["
@@ -477,6 +487,7 @@ class Inputs {
      \param 
    **/
   std::string get_advection_neutrals_vertical();
+  bool get_advection_neutrals_bulkwinds();
 
 
   /**********************************************************************

include/ions.h

@@ -20,7 +20,7 @@ class Ions {
     std::string cName;
     precision_t mass;
     int charge;
-
+    int vibe;
     int DoAdvect;
 
     std::vector<precision_t> nu_ion_neutral_coef;
@@ -39,8 +39,10 @@ class Ions {
     // Sources and Losses:
 
     arma_cube density_scgc;
+    arma_cube newDensity_scgc;
     std::vector<arma_cube> par_velocity_vcgc;
     std::vector<arma_cube> perp_velocity_vcgc;
+    std::vector<arma_cube> velocity_vcgc;
 
     arma_cube temperature_scgc;
     arma_cube conduction_scgc;
@@ -57,6 +59,11 @@ class Ions {
   arma_cube temperature_scgc;
   arma_cube conduction_scgc;
   arma_cube electron_temperature_scgc;
+  arma_cube rho_scgc;
+  arma_cube mean_major_mass_scgc;
+  std::vector<arma_cube> cMax_vcgc;
+  arma_cube sound_scgc;
+  arma_cube gamma_scgc;
 
   // This is the vector that will contain all of the different species:
   std::vector<species_chars> species;
@@ -98,6 +105,12 @@ class Ions {
   void init_ion_temperature(Neutrals neutrals, Grid grid);
   void set_floor();
   void fill_electrons();
+  void calc_sound_speed();
+  void calc_cMax();
+  bool set_bcs(Grid grid, Times time, Indices indices);
+  bool set_upper_bcs(Grid grid);
+  bool set_lower_bcs(Grid grid, Times time, Indices indices);
+
   int get_species_id(std::string name);
   void calc_efield(Grid grid);
   void calc_exb_drift(Grid grid);
@@ -111,5 +124,14 @@ class Ions {
   bool check_for_nonfinites();
   void nan_test(std::string variable);
   bool restart_file(std::string dir, bool DoRead);
+
+  /**********************************************************************
+     \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);
+
 };
 #endif  // INCLUDE_IONS_H_

include/neutrals.h

@@ -72,6 +72,8 @@ class Neutrals {
 
     /// concentration (density of species / total density)
     arma_cube concentration_scgc;
+    // mass concentration (mass * density of species / rho)
+    arma_cube mass_concentration_scgc;
 
     /// Diffusion through other neutral species:
     std::vector<float> diff0;
@@ -169,6 +171,9 @@ class Neutrals {
   /// Eddy Diffusion
   arma_cube kappa_eddy_scgc;
 
+  /// Viscosity
+  arma_cube viscosity_scgc;
+
   /// O cooling 
   arma_cube O_cool_scgc;
 
@@ -183,6 +188,15 @@ class Neutrals {
 
   // Source terms:
 
+  // Bulk acceleration due to collisions with ions:
+  std::vector<arma_cube> acc_ion_collisions;
+
+  // Bulk acceleration due to coriolis
+  std::vector<arma_cube> acc_coriolis;
+
+  // Total bulk acceleration
+  std::vector<arma_cube> acc_sources_total;
+
   /// Bulk neutral thermal conduction temperature change rate (K/s)
   arma_cube conduction_scgc;
 
@@ -192,6 +206,15 @@ class Neutrals {
   /// Bulk neutral chemical heating temperatuare change (K/s)
   arma_cube heating_chemical_scgc;
 
+  // Bulk neutral collisional heating with ions (K/s)
+  arma_cube heating_ion_friction_scgc;
+
+  // Bulk neutral collisional heating with ions (K/s)
+  arma_cube heating_ion_heat_transfer_scgc;
+
+  // Total heating sources
+  arma_cube heating_sources_total;
+
   /// Nuetral gas direct absorption heating efficiency (~5%)
   precision_t heating_efficiency;
 
@@ -294,6 +317,11 @@ class Neutrals {
    **/
   void calc_scale_height(Grid grid);
 
+  /**********************************************************************
+     \brief Calculate the viscosity coefficient
+   **/
+  void calc_viscosity();
+  
   /**********************************************************************
      \brief Calculate the eddy diffusion coefficient in valid pressure
    **/
@@ -303,7 +331,13 @@ class Neutrals {
      \brief Calculate the concentration for each species (species ndensity / total ndensity)
    **/
   void calc_concentration();
-    
+
+  /**********************************************************************
+     \brief Calculate the density of each species from the mass concentration 
+            for each species and rho (ndensity = con * rho / mass)
+   **/
+  void calc_density_from_mass_concentration();
+  
   /**********************************************************************
      \brief Calculate the bulk mean major mass
    **/
@@ -335,25 +369,24 @@ class Neutrals {
   void calc_cMax();
 
   /**********************************************************************
-     \brief Calculate dt (cell size / cMax) in each direction, and take min
-     \param dt returns the neutral time-step
+     \brief Calculate the chapman integrals for the individual species
      \param grid The grid to define the neutrals on
    **/
-  precision_t calc_dt(Grid grid);
-  precision_t calc_dt_cubesphere(Grid grid);
+  void calc_chapman(Grid grid);
 
   /**********************************************************************
-     \brief Calculate the chapman integrals for the individual species
+     \brief Calculate the neutral bulk vertical thermal conduction
      \param grid The grid to define the neutrals on
+     \param time The times within the model (dt is needed)
    **/
-  void calc_chapman(Grid grid);
+  void update_temperature(Grid grid, Times time);
 
   /**********************************************************************
-     \brief Calculate the neutral bulk vertical thermal conduction
+     \brief Calculate the neutral bulk horizontal viscosity
      \param grid The grid to define the neutrals on
      \param time The times within the model (dt is needed)
    **/
-  void calc_conduction(Grid grid, Times time);
+  void update_horizontal_velocity(Grid grid, Times time);
 
   /**********************************************************************
      \brief Calculate the O radiative cooling
@@ -368,8 +401,10 @@ class Neutrals {
   /**********************************************************************
      \brief Add all of the neutral source terms to each of the equations
      \param time The times within the model (dt is needed)
+     \param planet Need things like rotation rate
+     \param grid Need things like radius
    **/
-  void add_sources(Times time);
+  void add_sources(Times time, Planets planet, Grid grid);
 
   /**********************************************************************
      \brief Set boundary conditions for the neutrals
@@ -415,7 +450,7 @@ class Neutrals {
   /*****************************************************************************
       \brief  Checks for nans and +/- infinities in density, temp, and velocity
     **/
-  bool check_for_nonfinites();
+  bool check_for_nonfinites(std::string location);
 
   /**********************************************************************
       \brief Checks for nans in the specified variable
@@ -491,6 +526,11 @@ class Neutrals {
    **/
 
   bool advect_vertical(Grid grid, Times time);
+
+  arma_vec calc_friction_one_cell(int64_t iLong, int64_t iLat, int64_t iAlt,
+				   arma_vec &vels);
+
+  void calc_neutral_friction();  
 
 };