Corrected problem with testing suite in response to JOSS review. Tryi…

…ng to extract dictionary that dosnt exist anymore and isnt required for testing. Also noticed the testing suite was using data from a previous setup and thus slightly different. This has been changed.

.gitignore

@@ -30,6 +30,12 @@
 *.tar
 *.zip
 
+# Text files #
+######################
+*.txt
+!mechanism_files/*.txt
+!Docker_README.txt
+
 # Logs and databases #
 ######################
 *.log
@@ -60,7 +66,7 @@ Thumbs.db
 !images/Example_SOA_simulation.png
 !images/Example_deafult_gas_simulation.png
 *.jpeg
-
+*.gif
 
 
 # Folders

Aerosol/Property_calculation.py

@@ -33,6 +33,7 @@
 from umansysprop import critical_properties
 from umansysprop import liquid_densities
 from umansysprop import partition_models
+from umansysprop import groups
 from umansysprop import activity_coefficient_models_dev as aiomfac
 from umansysprop.forms import CoreAbundanceField
 import pdb
@@ -68,6 +69,8 @@ def Pure_component1(num_species,species_dict,species_dict2array,Pybel_object_dic
 
     y_density_array=[1000.0]*num_species
     y_mw=[200.0]*num_species
+    o_c=[0.0]*num_species
+    h_c=[0.0]*num_species
     sat_vp=[100.0]*num_species
     sat_vp_org=dict() #Recorded seperately and will not include the extension for water. This is used for any checks with equilibrium partitioning predictions
     y_gas=[0.0]*num_species
@@ -114,6 +117,10 @@ def Pure_component1(num_species,species_dict,species_dict2array,Pybel_object_dic
             #y_density_array.append(1400.0)
             y_mw[species_dict2array[compound]]=(Pybel_object_dict[SMILES_dict[compound]].molwt)
             #In the following you will need to select which vapour pressure method you like.
+
+            groups_dict=groups.composition(Pybel_object_dict[SMILES_dict[compound]])
+            o_c[species_dict2array[compound]]=groups_dict['O']/groups_dict['C']
+            h_c[species_dict2array[compound]]=groups_dict['H']/groups_dict['C']
 
             # Calculate boiling point
             b = boiling_point(Pybel_object_dict[SMILES_dict[compound]])
@@ -142,6 +149,8 @@ def Pure_component1(num_species,species_dict,species_dict2array,Pybel_object_dic
     return_dict=dict()
     return_dict['y_density_array']=y_density_array
     return_dict['y_mw']=y_mw
+    return_dict['o_c']=o_c
+    return_dict['h_c']=h_c
     return_dict['sat_vp']=sat_vp
     return_dict['Delta_H']=Delta_H
     return_dict['Latent_heat_gas']=Latent_heat_gas   

ODE_solver.py

@@ -76,7 +76,7 @@ def dydt_func(t,y):
         #pdb.set_trace()
         #Here we use the pre-created Numba based functions to arrive at our value for dydt
         # Calculate time of day
-        time_of_day_seconds=start_time+t
+        time_of_day_seconds=start_time+t 
 
         # make sure the y array is not a list. Assimulo uses lists
         y_asnumpy=numpy.array(y)
@@ -133,6 +133,41 @@ def dydt_func(t,y):
 
         return dydt
 
+    #-------------------------------------------------------------------------------------
+    #-------------------------------------------------------------------------------------
+    # define jacobian function to be called
+    def jac(t,y):    
+
+        """
+        This function defines Jacobian of the ordinary differential equations [ODEs] to be solved
+        input:
+        • t - time variable [internal to solver]
+        • y - array holding concentrations of all compounds in both gas and particulate [molecules/cc]
+        output:
+        dydt_dydt - the N_compounds x N_compounds matrix of Jacobian values
+        """
+
+        # Different solvers might call jacobian at different stages, so we have to redo some calculations here 
+        # Calculate time of day
+        time_of_day_seconds=start_time+t
+
+        # make sure the y array is not a list. Assimulo uses lists
+        y_asnumpy=numpy.array(y)
+
+        #Calculate the concentration of RO2 species, using an index file created during parsing
+        RO2=numpy.sum(y[RO2_indices])
+
+        #Calculate reaction rate for each equation.
+        # Note that H2O will change in parcel mode
+        rates=evaluate_rates(time_of_day_seconds,RO2,H2O,temp,numpy.zeros((equations)),numpy.zeros((63)))
+        #pdb.set_trace()
+        # Now use reaction rates with the loss_gain matrix to calculate the final dydt for each compound
+        # With the assimulo solvers we need to output numpy arrays
+        dydt_dydt=jacobian_function(rates,y_asnumpy,numpy.zeros((len(y_asnumpy),len(y_asnumpy))))
+        #pdb.set_trace()
+
+        return dydt_dydt
+
 
     #-------------------------------------------------------------------------------------
 
@@ -141,6 +176,7 @@ def dydt_func(t,y):
     from Rate_coefficients_numba import evaluate_rates 
     from Reactants_conc_numba import reactants as reactant_product
     from Loss_Gain_numba import dydt as dydt_eval
+    from Jacobian_numba import jacobian as jacobian_function
 
     # 'Unpack' variables from input_dict
     species_dict=input_dict['species_dict']
@@ -210,7 +246,7 @@ def dydt_func(t,y):
 
         # Define ODE parameters. 
         # Initial steps might be slower than mid-simulation. It varies.
-        #exp_mod.jac = dydt_jac
+        exp_mod.jac = jac
         # Define which ODE solver you want to use
         exp_sim = CVode(exp_mod) 
         tol_list=[1.0e-3]*num_species
@@ -222,7 +258,7 @@ def dydt_func(t,y):
         # Use of a jacobian makes a big differece in simulation time. This is relatively 
         # easy to define for a gas phase - not sure for an aerosol phase with composition
         # dependent processes. 
-        exp_sim.usejac = False # To be provided as an option in future update. See Fortran variant for use of Jacobian
+        exp_sim.usejac = True # To be provided as an option in future update. See Fortran variant for use of Jacobian
         #exp_sim.fac1 = 0.05
         #exp_sim.fac2 = 50.0
         exp_sim.report_continuously = True

f2py/ODE_solver.py

@@ -191,11 +191,12 @@ def jacobian(t,y):
         # Define ODE parameters. 
         # Initial steps might be slower than mid-simulation. It varies.
         #exp_mod.jac = dydt_jac
+        exp_mod.jac = jacobian
         # Define which ODE solver you want to use
         exp_sim = CVode(exp_mod) 
         tol_list=[1.0e-3]*num_species
         exp_sim.atol = tol_list #Default 1e-6
-        exp_sim.rtol = 0.03 #Default 1e-6
+        exp_sim.rtol = 1e-6 #Default 1e-6
         exp_sim.inith = 1.0e-6 #Initial step-size
         #exp_sim.discr = 'Adams'
         exp_sim.maxh = 100.0

test/test_modules.py

@@ -107,12 +107,13 @@ def setup(filename):
     # you will need to change the reference here
     outputdict=Parse_eqn_file.extract_smiles_species(outputdict,'../Aerosol/MCM.xml')
 
+    #pdb.set_trace()
     # Collect the dictionaries generated
-    reaction_dict=outputdict['reaction_dict']
+    #reaction_dict=outputdict['reaction_dict']
     rate_dict=outputdict['rate_dict']
     rate_dict_fortran=rate_coeff_conversion.convert_rate_mcm_fortran(rate_dict)
     rate_dict_reactants=outputdict['rate_dict_reactants']
-    rate_def=outputdict['rate_def']
+    #rate_def=outputdict['rate_def']
     loss_dict=outputdict['loss_dict']
     gain_dict=outputdict['gain_dict']
     stoich_dict=outputdict['stoich_dict']
@@ -171,7 +172,7 @@ def setup(filename):
     #Create a number concentration for a lognormal distribution
     N_perbin,x=Size_distributions.lognormal(num_bins,total_conc,meansize,std,lowersize,uppersize)
 
-    openMP=True
+    openMP=False
 
     # Convert the rate coefficient expressions into Fortran commands
     print("Converting rate coefficient operation into Fortran file")
@@ -181,23 +182,27 @@ def setup(filename):
     Parse_eqn_file.write_rate_file_fortran(filename,rate_dict_fortran,openMP)    
     print("Compiling rate coefficient file using f2py")
     #Parse_eqn_file.write_rate_file(filename,rate_dict,mcm_constants_dict)
-    os.system("python f2py_rate_coefficient.py build_ext --inplace")
+    #os.system("python f2py_rate_coefficient.py build_ext --inplace")
+    os.system('f2py -c -m rate_coeff_f2py Rate_coefficients.f90 --f90flags="-O3 -ffast-math -fopenmp" -lgomp')
 
     # Create Fortran file for calculating prodcts all of reactants for all reactions
     print("Creating Fortran file to calculate reactant contribution to equation")
     Parse_eqn_file.write_reactants_indices_fortran(filename,equations,species_dict2array,rate_dict_reactants,loss_dict,openMP)
     print("Compiling reactant product file using f2py")
-    os.system("python f2py_reactant_conc.py build_ext --inplace")
+    #os.system("python f2py_reactant_conc.py build_ext --inplace")
+    os.system('f2py -c -m reactants_conc_f2py Reactants_conc.f90 --f90flags="-O3 -ffast-math -fopenmp" -lgomp')
 
     # Create Fortran file for calculating dy_dt
     print("Creating Fortran file to calculate dy_dt for each reaction")
     Parse_eqn_file.write_loss_gain_fortran(filename,equations,num_species,loss_dict,gain_dict,species_dict2array,openMP)
     print("Compiling dydt file using f2py")
-    os.system("python f2py_loss_gain.py build_ext --inplace")
+    #os.system("python f2py_loss_gain.py build_ext --inplace")
+    os.system('f2py -c -m loss_gain_f2py Loss_Gain.f90 --f90flags="-O3 -ffast-math -fopenmp" -lgomp')
 
     Parse_eqn_file.write_gas_jacobian_fortran(filename,equations,num_species,loss_dict,gain_dict,species_dict2array,rate_dict_reactants,openMP)
     print("Compiling jacobian function using f2py")      
-    os.system("python f2py_jacobian.py build_ext --inplace")
+    #os.system("python f2py_jacobian.py build_ext --inplace")
+    os.system('f2py -c -m jacobian_f2py Jacobian.f90 --f90flags="-O3 -ffast-math -fopenmp" -lgomp')
 
     # Create .npy file with indices for all RO2 species
     print("Creating file that holds RO2 species indices")
@@ -396,6 +401,11 @@ def test_rates_numba(self):
         rates_numba_base=numpy.load('./data/'+filename+'rates_numba_base.npy')
         rates_numba=numpy.load(filename+'rates_numba.npy')
         npt.assert_almost_equal(rates_numba_base, rates_numba)
+
+    def rates_fortran_numba(self):
+        rates_fortran_base=numpy.load('./data/'+filename+'rates_fortran_base.npy')
+        rates_numba_base=numpy.load('./data/'+filename+'rates_numba_base.npy')
+        npt.assert_almost_equal(rates_fortran_base, rates_numba_base)
 
     def test_reactants_numba(self):
         reactants_numba_base=numpy.load('./data/'+filename+'reactants_numba_base.npy')