Merge pull request Planet-INC#2 from SylvainPlessis/planet_grins

Planet grins

src/core/include/planet/atmospheric_mixture.h

@@ -470,8 +470,8 @@ namespace Planet
       antioch_assert_equal_to(upper_fluxes.size(),molar_concentrations.size());
       for(unsigned int s = 0; s < _neutral_composition.n_species(); s++)
       {
-          CoeffType ms = _neutral_composition.M(s) * 1e-3 / Antioch::Constants::Avogadro<CoeffType>(); //to kg.mol-1 then kg
-          upper_fluxes[s] = this->Jeans_flux(ms,molar_concentrations[s],_temperature.neutral_temperature(_zmax),_zmax); // km-2/s
+          CoeffType ms = _neutral_composition.M(s) * Antioch::constant_clone(ms,1e-3) / Antioch::Constants::Avogadro<CoeffType>(); //to kg.mol-1 then kg
+          upper_fluxes[s] = - this->Jeans_flux(ms,molar_concentrations[s],_temperature.neutral_temperature(_zmax),_zmax) * Antioch::constant_clone(ms,1e-3); // cm-3.km/s, escaping flux, term < 0
       }
   }
 
@@ -483,7 +483,7 @@ namespace Planet
       antioch_assert_less(s,molar_concentrations.size());
 
       CoeffType ms = _neutral_composition.M(s) * 1e-3 / Antioch::Constants::Avogadro<CoeffType>(); //to kg.mol-1 then kg
-      CoeffType value = this->Jeans_flux(ms, molar_concentrations[s],_temperature.neutral_temperature(_zmax),_zmax) * 1e12; // to km-3
+      CoeffType value = - this->Jeans_flux(ms, molar_concentrations[s],_temperature.neutral_temperature(_zmax),_zmax) * Antioch::constant_clone(ms,1e-3); // to cm-3.km.s-1, escaping flux, term < 0
       return value;
   }
 

src/grins_interface/include/planet/planet_physics.h

@@ -201,7 +201,6 @@ namespace Planet
     PlanetPhysicsEvaluator<CoeffType,VectorCoeffType,MatrixCoeffType> evaluator(_helper);
 
 //    std::cout << "Element #" << context.get_elem().id() << std::endl;
-
     for (unsigned int qp=0; qp != n_qpoints; qp++)
       {
         const libMesh::Number r = s_qpoint[qp](0);
@@ -211,10 +210,12 @@ namespace Planet
 
         std::vector<libMesh::Number> molar_concentrations(this->_n_species, 0);
         std::vector<libMesh::Number> dmolar_concentrations_dz(this->_n_species, 0);
+
         for(unsigned int s=0; s < this->_n_species; s++ )
           {
             molar_concentrations[s] = context.interior_value(this->_species_vars[s],qp);
             dmolar_concentrations_dz[s] = context.interior_gradient(this->_species_vars[s],qp)(0);
+
           }
 
         evaluator.compute(molar_concentrations, dmolar_concentrations_dz, // {n}_s, {dn_dz}_s
@@ -230,13 +231,13 @@ namespace Planet
             libMesh::Real omega = evaluator.diffusion_term(s);
 
             libMesh::Real omega_dot = evaluator.chemical_term(s);
-//std::cout <<"z = " << z << ", omega = " << omega << ", omega_dot = " << omega_dot << std::endl;
+//std::cout << "omega = " << omega << ", omega_dot = " << omega_dot << std::endl;
 
             for(unsigned int i=0; i != n_s_dofs; i++)
               {
-                Fs(i) += (  omega_dot*s_phi[i][qp] 
- //                         + 2*omega*s_phi[i][qp] 
-                            - omega*s_grad_phi[i][qp](0) )*jac;
+                Fs(i) += (  omega_dot * s_phi[i][qp]  // chemistry
+                          + omega * s_grad_phi[i][qp](0) //diffusion
+                          )*jac;
 
                 if( compute_jacobian )
                   {

src/grins_interface/include/planet/planet_physics_helper.h

@@ -165,7 +165,7 @@ namespace Planet
                           const std::string& file_flyby,
                           const std::string& root_input) const;
 
-    void read_crossSection( const std::string &file,
+    void read_cross_section( const std::string &file,
                             VectorCoeffType &lambda, VectorCoeffType &sigma ) const;
 
     void read_hv_flux(VectorCoeffType& lambda, VectorCoeffType& phy1AU, const std::string &file) const;
@@ -442,8 +442,8 @@ namespace Planet
 
     this->read_hv_flux(_lambda_hv, _phy1AU, input_hv);
 
-    this->read_crossSection(input_N2,  lambda_N2,  sigma_N2);
-    this->read_crossSection(input_CH4, lambda_CH4, sigma_CH4);
+    this->read_cross_section(input_N2,  lambda_N2,  sigma_N2);
+    this->read_cross_section(input_CH4, lambda_CH4, sigma_CH4);
 
     /* here only N2 and CH4 absorb */
     _tau->add_cross_section( lambda_N2, sigma_N2, Antioch::Species::N2,
@@ -760,7 +760,7 @@ namespace Planet
                                                                                         const std::string& file_flyby,
                                                                                         const std::string& root_input) const
   {
-//to do, change to a GetPot
+//TODO, change to a GetPot
     std::ifstream flyby(file_flyby.c_str());
     if( !flyby )
       {
@@ -801,7 +801,7 @@ namespace Planet
   }
 
   template<typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
-  void PlanetPhysicsHelper<CoeffType,VectorCoeffType,MatrixCoeffType>::read_crossSection( const std::string &file,
+  void PlanetPhysicsHelper<CoeffType,VectorCoeffType,MatrixCoeffType>::read_cross_section( const std::string &file,
                             VectorCoeffType &lambda, VectorCoeffType &sigma ) const
   {
     std::string line;

test/input/solver_test.in.in

@@ -3,13 +3,15 @@
 
 medium = 'N2 CH4'
 
-neutral_species = 'H H2 C CH (1)CH2 (3)CH2 CH3 CH4 C2H2 C2H4 C2H6 N(4S) N2 NH CN
+neutral_species = 'N2 CH4 H H2 C CH (1)CH2 (3)CH2 CH3 C2H2 C2H4 C2H6 N(4S) NH CN
                    HCN H2CN C2N NCCN HC3N N(2D) N(2P) C2H2** HC3N* Ar C2 C2H
                    CH2C2H2 CH3C2H C4H2 C2H5 C3H2 C3H3 C3H5 C3H6 C3H7 C3H8 C4H3
                    C4H4 C4H6 C4H10 C6H2 C6H5 cC6H6 nC6H6 C6H6 nC6H7 cC6H7 C2H3
                    C2H3CN CH2NH CHCN CH3CN H2C3N HC5N C2N2 HC2N2 C4N2 [14]N[15]N
                    [15]N HC[15]N H2C[15]N C[15]N [15]NH [15]N(2D)'
 
+sparse_neutral_species = 'N2 CH4 H'
+
 # This example, the ionic species and neutral species are the same.
 # If ionic species were different, would just list them here
 ionic_species = '${Planet/neutral_species}'
@@ -22,6 +24,7 @@ input_N2 = '@abs_top_srcdir@/test/input/N2_hv_cross-sections.dat'
 input_CH4 = '@abs_top_srcdir@/test/input/CH4_hv_cross-sections.dat'
 input_hv = '@abs_top_srcdir@/test/input/hv_SSI.dat'
 input_reactions_elem = '@abs_top_srcdir@/test/input/neutral_reactions.bimol'
+sparse_input_reactions_elem = '@abs_top_srcdir@/test/input/neutral_reactions.solv_test'
 input_reactions_fall = '@abs_top_srcdir@/test/input/neutral_reactions.falloff'
 
 zmin = '600.0'
@@ -74,10 +77,12 @@ use_numerical_jacobians_only = 'true'
 
 # Options for time solvers
 [unsteady-solver]
-transient = 'false' 
+transient = 'true' 
 theta = 1.0
-n_timesteps = 700
-deltat = 50.0
+n_timesteps = 100
+deltat = 1000000
+target_tolerance = -1
+max_growth = 100
 
 []
 

test/physics_helper_unit.C

@@ -115,7 +115,9 @@ void parse_equation(std::vector<std::string> &reactants, std::vector<std::string
    if(molecules.size() != 2)antioch_error();
 
    Antioch::SplitString(molecules[0],"+",reactants,false);
+   if(reactants.empty())reactants.push_back(molecules[0]);
    Antioch::SplitString(molecules[1],"+",products,false);
+   if(products.empty())products.push_back(molecules[1]);
    shave_strings(reactants);
    shave_strings(products);
    stoi_reac.resize(reactants.size(),1);