Added energy and power density to output

bat_can.py

@@ -162,7 +162,6 @@ def model_run(sim):
 
     for sim in parameters['simulations']:
         model = importlib.import_module('.'+sim['type'], package='simulations')
-
         solution = model.plot(an, sep, ca, parameters, sim)
 
 #===========================================================================

bat_can_fit.py

@@ -8,7 +8,7 @@
 
 """
 # Import modules
-from datetime import datetime   
+from datetime import datetime
 from functools import partial
 import importlib # allows us to import from user input string.
 import itertools
@@ -23,7 +23,7 @@
 from submodels.fitting import voltage_capacity as fit
 from bat_can_init import initialize
 
-# This is the main function that runs the model.  We define it this way so it 
+# This is the main function that runs the model.  We define it this way so it
 # is called by "main," below:
 def bat_can(input, cores, print_flag):
     if input is None:
@@ -55,32 +55,32 @@ def bat_can(input, cores, print_flag):
     #===========================================================================
     #   CREATE ELEMENT CLASSES AND INITIAL SOLUTION VECTOR SV_0
     #===========================================================================
-    # For each element (anode 'an', separator 'sep', cathode 'ca') the 'class' 
-    # variable from the inputs tells what kind of anode, separator, or cathode 
-    # it is, and points to a '.py' file in this directory.  We import that 
-    # module, and then run its 'initialize' routine to create an intial 
+    # For each element (anode 'an', separator 'sep', cathode 'ca') the 'class'
+    # variable from the inputs tells what kind of anode, separator, or cathode
+    # it is, and points to a '.py' file in this directory.  We import that
+    # module, and then run its 'initialize' routine to create an intial
     # solution vector and an object that stores needed parameters.
     # import single_particle_electrode as an_module_0
-    
-    an_module = importlib.import_module('electrode_models.' 
+
+    an_module = importlib.import_module('electrode_models.'
         + an_inputs['class'])
-    an = an_module.electrode(input_file, an_inputs, sep_inputs, ca_inputs, 
+    an = an_module.electrode(input_file, an_inputs, sep_inputs, ca_inputs,
             'anode', parameters, offset=0)
-    
-    sep_module = importlib.import_module('separator_models.' 
+
+    sep_module = importlib.import_module('separator_models.'
         + sep_inputs['class'])
-    sep = sep_module.separator(input_file, sep_inputs, parameters, 
+    sep = sep_module.separator(input_file, sep_inputs, parameters,
             offset=an.n_vars)
-    
+
     # Check to see if the anode object needs to adjust the separator properties:
     sep = an.adjust_separator(sep)
-    ca_module = importlib.import_module('electrode_models.' 
+    ca_module = importlib.import_module('electrode_models.'
         + ca_inputs['class'])
-    ca = ca_module.electrode(input_file, ca_inputs, sep_inputs, an_inputs, 
+    ca = ca_module.electrode(input_file, ca_inputs, sep_inputs, an_inputs,
         'cathode', parameters, offset= an.n_vars+sep.n_vars*sep.n_points)
 
-    # Check to see if the cathode object needs to adjust the separator 
-    # properties:    
+    # Check to see if the cathode object needs to adjust the separator
+    # properties:
     sep = ca.adjust_separator(sep)
 
     # Initialize the solution vector:
@@ -96,23 +96,23 @@ def bat_can(input, cores, print_flag):
     #===========================================================================
     #   RUN THE SIMULATION
     #===========================================================================
-    # The inputs tell us what type of experiment we will simulate.  Load the 
+    # The inputs tell us what type of experiment we will simulate.  Load the
     # module, then call its 'run' function:
 
     # If the user requests a specific initailization routine, run that first:
     if 'initialize' in parameters and parameters['initialize']['enable']:
         if parameters['initialize']['type'] == 'open-circuit':
-            model = importlib.import_module('.'+'CC_cycle', 
+            model = importlib.import_module('.'+'CC_cycle',
                 package='simulations')
-            
-            sim = {'i_ext': None, 'C-rate': 0.0, 'n_cycles': 0, 
-                'first-step': 'discharge', 'equilibrate': 
-                {'enable': True, 'time':  5}, 'phi-cutoff-lower': 2.0, 
+
+            sim = {'i_ext': None, 'C-rate': 0.0, 'n_cycles': 0,
+                'first-step': 'discharge', 'equilibrate':
+                {'enable': True, 'time':  5}, 'phi-cutoff-lower': 2.0,
                 'phi-cutoff-upper': 4.8, 'init':True}
-            
+
             solution = model.run(SV_0, an, sep, ca, algvars, parameters, sim)
-            
-            # Save final state as the initial state for all subsequent 
+
+            # Save final state as the initial state for all subsequent
             # simulation steps:
             SV_0 = model.final_state(solution)
 
@@ -121,53 +121,53 @@ def bat_can(input, cores, print_flag):
 
     def run_model(x_guess, final_flag = False):
         print('x = ',x_guess)
-    
+
         # Set the current fitting parameter values:
         for i, x in enumerate(x_guess):
             if fit_params[i]['type'] == 'elyte-transport':
                 sep.D_k[sep.elyte_obj.species_index(fit_params[i]['species'])] \
                     = x
             elif fit_params[i]['type'] == 'cathode-kinetics':
-                ca.surf_obj.set_multiplier(x, 
+                ca.surf_obj.set_multiplier(x,
                     fit_params[i]['reaction-index'])
             elif fit_params[i]['type'] == 'air-kinetics':
-                ca.gas_elyte_obj.set_multiplier(x, 
+                ca.gas_elyte_obj.set_multiplier(x,
                     fit_params[i]['reaction-index'])
             elif fit_params[i]['type'] == 'cathode-microstructure':
                 setattr(ca, fit_params[i]['parameter'], x)
 
         if final_flag:
             SSR_net = 0
             icolor = 0
-            fit_fig, fit_axs = plt.subplots(1, 1, sharex=True, 
+            fit_fig, fit_axs = plt.subplots(1, 1, sharex=True,
                 gridspec_kw = {'wspace':0, 'hspace':0})
-        
+
             fit_fig.set_size_inches((5.0, 2.25))
 
             # pool = mp.Pool(processes = int(cores))
             # SSR_net = pool.map(run, (list(parameters['simulations'])))
             for sim in parameters['simulations']:
-                SSR_net += run(sim, SV_0, algvars, parameters, an, ca, sep, 
+                SSR_net += run(sim, SV_0, algvars, parameters, an, ca, sep,
                                final_flag, fit_fig, fit_axs, icolor)
                 icolor += 1
         else:
             with Pool(processes = int(cores)) as pool:
-                SSR_net = pool.map(partial(run, SV_0=SV_0, algvars=algvars, 
-                            parameters=parameters, an=an, sep=sep, ca=ca), 
+                SSR_net = pool.map(partial(run, SV_0=SV_0, algvars=algvars,
+                            parameters=parameters, an=an, sep=sep, ca=ca),
                             list(parameters['simulations']))
 
         print('SSR = ', np.sum(SSR_net))
-        
+
         if final_flag:
             fit_axs.annotate(f"SSR = 1.0", xy=(0,0))#,
                 # xytext=(0.5, 0.5), textcoords='axes fraction')#, fontsize = 8)
-                
+
             if len(parameters['simulations']) == 1:
-                savename = (parameters['output'] +'_' 
+                savename = (parameters['output'] +'_'
                     + sim['outputs']['save-name'] )
             else:
                 savename = (parameters['output'])
-                
+
             fit_fig.savefig(savename + '/fit.pdf')
             plt.show()
 
@@ -183,73 +183,73 @@ def run_model(x_guess, final_flag = False):
     for sim in parameters['simulations']:
         sim['ref_data'] = pd.read_excel('data/' + sim['validation'])
 
-    
+
     if print_flag:
-        # Just print the fit of the starting guess - do not run the fitting 
+        # Just print the fit of the starting guess - do not run the fitting
         # routine:
         result = run_model(x_start, final_flag=True)
     else:
         # Fit to the reference data:
-        result = spo.minimize(run_model, x_start, bounds = x_bounds, 
+        result = spo.minimize(run_model, x_start, bounds = x_bounds,
                 options={'disp': True})
-    
+
         print(result)
 
         print("Best fit = ", result.x)
 
         run_model(result.x, final_flag=True)
 
     stop = timeit.default_timer()
-    print('Time: ', stop - start)  
+    print('Time: ', stop - start)
 
-def run(sim, SV_0, algvars, parameters, an, ca, sep, final_flag=False, 
+def run(sim, SV_0, algvars, parameters, an, ca, sep, final_flag=False,
             fit_fig=None, fit_axs=None, icolor=0):
-    
+
     sim_local = sim
     try:
-        model = importlib.import_module('.'+sim['type'], 
+        model = importlib.import_module('.'+sim['type'],
             package='simulations')
 
         sim['init'] = False
-        solution = model.run(SV_0, an, sep, ca, algvars, 
+        solution = model.run(SV_0, an, sep, ca, algvars,
             parameters, sim)
         # Read out the results:
         if final_flag:
             sim_local['outputs']['show-plots'] = False
-            results = model.output(solution, an, sep, ca, parameters, 
+            results = model.output(solution, an, sep, ca, parameters,
                 sim_local, plot_flag=True, return_flag=True, save_flag=True)
                 #     for sim in parameters['simulations']:
                 # model = importlib.import_module('.'+sim['type'], package='simulations')
             solution = model.plot(an, sep, ca, parameters, sim_local)
         else:
-            results = model.output(solution, an, sep, ca, parameters, 
+            results = model.output(solution, an, sep, ca, parameters,
                 sim, plot_flag=False, return_flag=True, save_flag=False)
         phi_sim = results['phi_ed'].to_numpy()[:,-1]
-            
+
         sim_data = np.array((results['capacity'].to_numpy(),
             phi_sim))
 
         # This function calculates the SSR for this simulation:
         # Scale to convert capacity units in validation data to mAh/cm2:
-        UnitsScale = 46.968 #todo #72
-        ssr_calc = fit.SSR(sim['ref_data'].to_numpy(), sim_data.T, 
+        UnitsScale = 1 #todo #72 46.968
+        ssr_calc = fit.SSR(sim['ref_data'].to_numpy(), sim_data.T,
                 units_scale = UnitsScale)
         print('SSR = ', ssr_calc)
-        
+
         ndata = len(parameters['simulations']) + 1
         cmap = plt.get_cmap('plasma')
-        
+
         color_ind = np.linspace(0,1,ndata)
         colors = list()
-        
+
         for i in np.arange(ndata):
             colors.append(cmap(color_ind[i]))
-  
+
         if final_flag:
-            fit_axs, fit_fig = fit.plot(sim['ref_data'].to_numpy(), 
-                sim_data.T, fit_axs, fit_fig, units_scale = UnitsScale, 
+            fit_axs, fit_fig = fit.plot(sim['ref_data'].to_numpy(),
+                sim_data.T, fit_axs, fit_fig, units_scale = UnitsScale,
                 color = colors[icolor])
-            
+
     except:
         # Assign a large penalty for failed parameter sets:
         ssr_calc = 1e23
@@ -265,12 +265,12 @@ def run(sim, SV_0, algvars, parameters, an, ca, sep, final_flag=False,
 if __name__ == '__main__':
     import argparse
 
-    # Currently, the only command line keyword enabled is --input, to specify 
+    # Currently, the only command line keyword enabled is --input, to specify
     # the input file location:
     parser = argparse.ArgumentParser()
     parser.add_argument('--input')
     parser.add_argument('--cores')
     parser.add_argument('--print', action='store_true')
     args = parser.parse_args()
-    
+
     bat_can(args.input, args.cores, args.print)

electrode_models/air_electrode.py

@@ -55,6 +55,8 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs,
 
         # Phase volume fractions
         eps_host = np.array([inputs['eps_host']])
+        eps_host = eps_host[0]
+
         if len(eps_host) == 1:
             self.eps_host = np.repeat(eps_host[0], self.n_points)
         elif len(eps_host) == self.n_points:
@@ -71,6 +73,8 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs,
         # radius, with no overlap.
         self.r_host = inputs['r_host']
         self.th_product = inputs['th_product']
+
+
         self.V_host = 4./3. * np.pi * (self.r_host)**3  # Volume of a single carbon / host particle [m3]
         self.A_host = 4. * np.pi * (self.r_host)**2 # Surface area of a single carbon / host particle [m2]
         # m2 of host-electrolyte interface / m3 of total volume [m^-1]
@@ -93,8 +97,7 @@ def __init__(self, input_file, inputs, sep_inputs, counter_inputs,
 
         # Microstructure-based transport scaling factor, based on Bruggeman
         # coefficient of -0.5:
-        self.elyte_microstructure = self.eps_elyte_init[0]**1.5 #This is used by the
-
+        self.elyte_microstructure = self.eps_elyte_init[0]**1.5 #This is used by the electrode_boundary_flux
         # SV_offset specifies the index of the first SV variable for the
         # electode (zero for anode, n_vars_anode + n_vars_sep for the cathode)
         self.SV_offset = offset
@@ -219,7 +222,7 @@ def residual(self, t, SV, SVdot, sep, counter, params):
         c_k_elyte = SV_loc[SVptr['C_k_elyte'][j]]
         eps_product = SV_loc[SVptr['eps_product'][j]]
         eps_elyte = 1. - eps_product - self.eps_host[j]
-
+        self.elyte_microstructure = eps_elyte**1.5
         # Read electrolyte fluxes at the separator boundary.  No matter the
         # electrode, the function returns a value where flux to the electrode
         # is considered positive. We multiply by `i_ext_flag` to get the
@@ -270,7 +273,7 @@ def residual(self, t, SV, SVdot, sep, counter, params):
             resid[SVptr['phi_ed'][j]] = i_io_in - i_io_out + i_el_in - i_el_out
 
             # Calculate available surface area (m2 interface per m3 electrode):
-            A_avail = self.A_init[j] - eps_product/self.th_product
+            A_avail = self.A_init[j] - eps_product/self.thickness_temp#self.th_product
             # Convert to m2 interface per m2 geometric area:
             A_surf_ratio = A_avail*self.dy
             # Multiplier to scale phase destruction rates.  As eps_product
@@ -363,7 +366,7 @@ def residual(self, t, SV, SVdot, sep, counter, params):
                     - params['potential'])
 
         # Calculate available surface area (m2 interface per m3 electrode):
-        A_avail = self.A_init[j] - eps_product/self.th_product
+        A_avail = self.A_init[j] - eps_product/self.thickness_temp#self.th_product
         # Convert to m2 interface per m2 geometric area:
         A_surf_ratio = A_avail*self.dy
 
@@ -413,6 +416,9 @@ def residual(self, t, SV, SVdot, sep, counter, params):
             + sdot_elyte_host * A_surf_ratio)
             * self.dyInv)/eps_elyte
 
+        eps_elyte = 1. - eps_product - self.eps_host[0]
+        self.elyte_microstructure = eps_elyte**1.5
+
         return resid
 
     def voltage_lim(self, SV, val):