Merge branch 'feature/2D_C3' into feature/temperature_effects

.devcontainer/Dockerfile

@@ -0,0 +1,16 @@
+FROM --platform=linux/amd64 mcr.microsoft.com/devcontainers/miniconda:0-3
+
+# Copy environment.yml (if found) to a temp location so we update the environment. Also
+# copy "noop.txt" so the COPY instruction does not fail if no environment.yml exists.
+COPY environment.yml* .devcontainer/noop.txt /tmp/conda-tmp/
+RUN if [ -f "/tmp/conda-tmp/environment.yml" ]; then umask 0002 && /opt/conda/bin/conda env update -n base -f /tmp/conda-tmp/environment.yml; fi \
+    && rm -rf /tmp/conda-tmp
+
+# [Optional] Uncomment to install a different version of Python than the default
+# RUN conda install -y python=3.6 \
+#     && pip install --no-cache-dir pipx \
+#     && pipx reinstall-all
+
+# [Optional] Uncomment this section to install additional OS packages.
+# RUN apt-get update && export DEBIAN_FRONTEND=noninteractive \
+#     && apt-get -y install --no-install-recommends <your-package-list-here>

.devcontainer/devcontainer.json

@@ -0,0 +1,24 @@
+// For format details, see https://aka.ms/devcontainer.json. For config options, see the
+// README at: https://github.com/devcontainers/templates/tree/main/src/miniconda
+{
+	"name": "Miniconda (Python 3)",
+	"build": { 
+		"context": "..",
+		"dockerfile": "Dockerfile"
+	}
+
+	// Features to add to the dev container. More info: https://containers.dev/features.
+	// "features": {},
+
+	// Use 'forwardPorts' to make a list of ports inside the container available locally.
+	// "forwardPorts": [],
+
+	// Use 'postCreateCommand' to run commands after the container is created.
+	// "postCreateCommand": "python --version",
+
+	// Configure tool-specific properties.
+	// "customizations": {},
+
+	// Uncomment to connect as root instead. More info: https://aka.ms/dev-containers-non-root.
+	// "remoteUser": "root"
+}

.devcontainer/noop.txt

@@ -0,0 +1,3 @@
+This file is copied into the container along with environment.yml* from the
+parent folder. This is done to prevent the Dockerfile COPY instruction from 
+failing if no environment.yml is found.

configs/params_LFP.cfg

@@ -5,7 +5,7 @@
 type = ACR
 discretization = 1e-9
 shape = C3
-thickness = 20e-9
+thickness = 40e-9
 
 [Material]
 muRfunc = LiFePO4

configs/params_LFP2D.cfg

@@ -0,0 +1,39 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = ACR2D
+discretization = 1e-9
+shape = C3
+discretization_ver = 0.75e-9
+# remember the it is half the thickness, check the normalizations
+thickness = 20e-9
+
+[Material]
+muRfunc = LiFePO42D
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+Omega_a = 1.8560e-20
+kappa_x = 5.0148e-10
+kappa_y = 10e-10
+Bx = 0.19e9
+By = 0.5e9
+rho_s = 1.3793e28
+D = 1e-14
+E_D = 0
+Dfunc = lattice
+dgammadc = 0e-30
+cwet = 0.98
+D_surf = 1e-16
+E_D_surf = 0
+
+[Reactions]
+surface_diffusion = true
+rxnType = BV
+k0 = 0.5
+E_A = 13000
+alpha = 0.5
+# Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0

configs/params_system.cfg

@@ -8,12 +8,12 @@
 profileType = CC
 # Battery (dis)charge c-rate (only used for CC), number of capacities / hr
 # (positive for discharge, negative for charge)
-Crate = 1
+Crate = 0.01
 # Optional nominal 1C current density for the cell, A/m^2
 # 1C_current_density = 12.705
 # Voltage cutoffs, V
 Vmax = 3.6
-Vmin = 2.0
+Vmin = 3.1
 # Battery applied voltage (only used for CV), V
 Vset = 0.12
 # Battery constant power (only used for CP), W/m^2
@@ -25,9 +25,10 @@ power = 1
 # This can have as many rows (segments) as desired.
 # Note: It's okay to leave commented lines within the segments list
 segments = [
-    (0.3, 0.4),
-    #            (0, 0.2),
-    (-0.5, 0.1),
+    (0.2,150),
+    # (5, 6),
+    (0, 30),
+    (3, 15),
     ]
 # Continuation directory. If false, begin a fresh simulation with the
 # specified input parameters here. Otherwise, this should be the
@@ -46,10 +47,11 @@ tend = 1.2e3
 # spacing between initial and final times with tsteps total outputs.
 tsteps = 200
 # Relative Tolerance
-relTol = 1e-6
+relTol = 1e-10
 # Absolute Tolerance
 absTol = 1e-6
 # Initial Temperature throughout electrode, K
+# Temperature, K
 T = 298
 # Nonisothermal: true for heat generation throughout electrode, false
 # for no heat generation
@@ -58,31 +60,31 @@ nonisothermal = true
 # (affects noise, particle size distribution). Set to true exactly
 # reproducible results -- useful for testing.
 # Options: true, false
-randomSeed = false
+randomSeed = true
 # Value of the random seed, must be an integer
 seed = 0
 # Data reporter: choice of mat (MATLAB), hdf5 (hdf5), or hdf5Fast (hdf5, without
 # printing internal variable concentrations) files. hdf5 files
 # are better for cycling, as they store less information and there is less
 # opening/rewriting of files. Default is mat
-dataReporter = hdf5
+dataReporter = mat
 # Series resistance, [Ohm m^2]
 Rser = 0.
 # Cathode, anode, and separator numer disc. in x direction (volumes in electrodes)
 # - Nvol_c must be >= 1
 # - If Nvol_c = 1 & Nvol_a = 0 & Nvol_s = 0, simulate a single volume with no
 #   separator and infinitely fast counter electrode
 # - If Nvol_a = 0, simulate a Li foil electrode
-Nvol_c = 10
-Nvol_s = 5
-Nvol_a = 10
+Nvol_c = 1
+Nvol_s = 0
+Nvol_a = 0
 # Number of particles per volume for cathode and anode
-Npart_c = 2
-Npart_a = 2
+Npart_c = 1
+Npart_a = 0
 
 [Electrodes]
 # The name of the parameter file describing the cathode particles
-cathode = params_LFP.cfg
+cathode = params_LFP2D.cfg
 # The name of the parameter file describing the anode particles
 # Used only if Nvol_a > 0
 anode = params_graphite_1param.cfg
@@ -98,8 +100,8 @@ Rfilm_foil = 0e-0
 #   sphere or cylinder -- radius
 # If using stddev = 0, set Npart = 1 for that electrode to avoid
 # wasting computational effort.
-mean_c = 100e-9
-stddev_c = 1e-9
+mean_c = 60e-9
+stddev_c = 0e-9
 mean_a = 100e-9
 stddev_a = 1e-9
 # Use a specific set of particle sizes for the distribution
@@ -109,7 +111,7 @@ specified_psd_c = False
 specified_psd_a = False
 # Initial electrode filling fractions
 # (for disch, anode starts full, cathode starts empty)
-cs0_c = 0.01
+cs0_c = 0.02
 cs0_a = 0.99
 
 [Conductivity]
@@ -152,7 +154,7 @@ entropy_heat_gen = False
 
 [Geometry]
 # Thicknesses, m
-L_c = 50e-6
+L_c = 30e-6
 L_a = 50e-6
 L_s = 25e-6
 # Volume loading percents of active material (volume fraction of solid

mpet/config/configuration.py

@@ -521,6 +521,10 @@ def _scale_electrode_parameters(self):
                 self[trode, 'lambda'] = self[trode, 'lambda'] / kT
             if self[trode, 'B'] is not None:
                 self[trode, 'B'] = self[trode, 'B'] / (kT * constants.N_A * self[trode, 'cs_ref'])
+            if self[trode, 'Bx'] is not None:
+                self[trode,'Bx'] = self[trode, 'Bx'] / (kT * constants.N_A * self[trode, 'cs_ref'])
+            if self[trode, 'By'] is not None:
+                self[trode,'By'] = self[trode, 'By'] / (kT * constants.N_A * self[trode, 'cs_ref'])
             for param in ['Omega_a', 'Omega_b', 'Omega_c', 'EvdW']:
                 value = self[trode, param]
                 if value is not None:
@@ -635,6 +639,9 @@ def _distr_part(self):
         self['psd_area'] = {}
         self['psd_vol'] = {}
         self['psd_vol_FracVol'] = {}
+        # if self[trode,'type'] in 'ACR2D':
+        self['psd_num_ver'] = {}
+        self['psd_len_ver'] = {}
 
         for trode in self['trodes']:
             solidType = self[trode, 'type']
@@ -645,6 +652,14 @@ def _distr_part(self):
             if not np.all(self['specified_psd'][trode]):
                 # If PSD is not specified, make a length-sampled particle size distribution
                 # Log-normally distributed
+                if self[trode,'type'] == 'ACR2D':
+                    np.random.seed(self.D_s['seed'])
+                    mean_t = self[trode,'thickness']
+                    stddev_t = mean_t*0
+                    var_t = stddev_t**2
+                    mu_t = np.log((mean_t**2) / np.sqrt(var_t + mean_t**2))
+                    sigma_t = np.sqrt(np.log(var_t/(mean_t**2) + 1))
+                    raw_t = np.random.lognormal(mu_t, sigma_t, size=(Nvol, Npart))
                 mean = self['mean'][trode]
                 stddev = self['stddev'][trode]
                 if np.allclose(stddev, 0., atol=1e-12):
@@ -667,6 +682,12 @@ def _distr_part(self):
             if solidType in ['ACR']:
                 psd_num = np.ceil(raw / solidDisc).astype(int)
                 psd_len = solidDisc * psd_num
+            elif solidType in ['ACR2D']:
+                solidDisc_ver = self[trode, 'discretization_ver']
+                psd_num = np.ceil(raw / solidDisc).astype(int)
+                psd_len = solidDisc * psd_num
+                psd_num_ver = np.ceil(raw_t/solidDisc_ver).astype(int)
+                psd_len_ver = solidDisc_ver * psd_num_ver
             elif solidType in ['CHR', 'diffn', 'CHR2', 'diffn2']:
                 psd_num = np.ceil(raw / solidDisc).astype(int) + 1
                 psd_len = solidDisc * (psd_num - 1)
@@ -704,6 +725,10 @@ def _distr_part(self):
             self['psd_area'][trode] = psd_area
             self['psd_vol'][trode] = psd_vol
             self['psd_vol_FracVol'][trode] = psd_frac_vol
+            # store values to config 2d
+            if self[trode,'type'] == 'ACR2D':
+                self['psd_num_ver'][trode] = psd_num_ver
+                self['psd_len_ver'][trode] = psd_len_ver
 
     def _G(self):
         """
@@ -757,15 +782,34 @@ def _indvPart(self):
                     kappa_ref = constants.k * constants.T_ref * cs_ref_part * plen**2  # J/m
                     gamma_S_ref = kappa_ref / plen  # J/m^2
                     # non-dimensional quantities
+                    if self[trode,'type'] == 'ACR2D':
+                        self[trode, 'indvPart']['N_ver'][i,j] = self['psd_num_ver'][trode][i,j]
                     if self[trode, 'kappa'] is not None:
                         self[trode, 'indvPart']['kappa'][i, j] = self[trode, 'kappa'] / kappa_ref
+                    if self[trode, 'kappa_x'] is not None:
+                        self[trode,'indvPart']['kappa_x'][i, j] = self[trode,'kappa_x'] / kappa_ref
+                    if self[trode, 'kappa_y'] is not None:
+                        pthick = self['psd_len_ver'][trode][i,j]
+                        self[trode,'indvPart']['kappa_y'][i, j] = self[trode,'kappa_y'] \
+                            / (kappa_ref * (pthick**2/plen**2))
                     if self[trode, 'dgammadc'] is not None:
                         nd_dgammadc = self[trode, 'dgammadc'] * cs_ref_part / gamma_S_ref
                         self[trode, 'indvPart']['beta_s'][i, j] = nd_dgammadc \
                             / self[trode, 'indvPart']['kappa'][i, j]
-                    self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] * self['t_ref'] / plen**2
+                    if self[trode,'type'] == 'ACR2D':
+                        pthick = self['psd_len_ver'][trode][i,j]
+                        self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] \
+                            * self['t_ref'] / pthick**2
+                        self[trode,'indvPart']['D_surf'][i, j] = self[trode,'D_surf'] \
+                            * self['t_ref'] / plen**2
+                    else:
+                        self[trode, 'indvPart']['D'][i, j] = self[trode, 'D'] \
+                            * self['t_ref'] / plen**2
                     self[trode, 'indvPart']['E_D'][i, j] = self[trode, 'E_D'] \
                         / (constants.k * constants.N_A * constants.T_ref)
+                    if self[trode, 'E_D_surf'] is not None:
+                        self[trode, 'indvPart']['E_D_surf'][i, j] = self[trode, 'E_D_surf'] \
+                            / (constants.k * constants.N_A * constants.T_ref)
                     self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \
                         / (constants.e * F_s_ref)
                     self[trode, 'indvPart']['E_A'][i, j] = self[trode, 'E_A'] \

mpet/config/constants.py

@@ -13,7 +13,7 @@
 #: General particle classification (1 var)
 two_var_types = ["diffn2", "CHR2", "homog2", "homog2_sdn"]
 #: General particle classification (2 var)
-one_var_types = ["ACR", "diffn", "CHR", "homog", "homog_sdn"]
+one_var_types = ["ACR", "diffn", "CHR", "homog", "homog_sdn", "ACR2D"]
 #: Reference concentration, mol/m^3 = 1M
 c_ref = 1000.
 #: Reaction rate epsilon for values close to zero
@@ -26,6 +26,7 @@
 #: subset of ``PARAMS_PER_TRODE``` that is defined for the separator as well
 PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp', 'k_h', 'cp', 'rhom']
 #: parameters that are defined for each particle, and their type
-PARAMS_PARTICLE = {'N': int, 'kappa': float, 'beta_s': float, 'D': float, 'k0': float,
+PARAMS_PARTICLE = {'N': int, 'N_ver':int, 'kappa': float, 'kappa_x': float, 'kappa_y': float,
+                   'beta_s': float, 'D': float,'D_surf': float, 'k0': float,
                    'Rfilm': float, 'delta_L': float, 'Omega_a': float, 'E_D': float,
-                   'E_A': float}
+                   'E_D_surf': float, 'E_A': float}

mpet/config/derived_values.py

@@ -241,7 +241,6 @@ def muR_ref(self, trode):
         muRfunc = props_am.muRfuncs(self.config, trode).muRfunc
         cs0bar = self.config['cs0'][trode]
         cs0 = np.array([cs0bar])
-
         solidType = self.config[trode, 'type']
         if solidType in constants.two_var_types:
             muR_ref = -muRfunc((cs0, cs0), (cs0bar, cs0bar), self.config['T'], 0.)[0][0]

mpet/config/schemas.py

@@ -145,6 +145,7 @@ def tobool(value):
                            'discretization': Use(float),
                            'shape': lambda x:
                                check_allowed_values(x, ["C3", "sphere", "cylinder", "homog_sdn"]),
+                           Optional('discretization_ver', default=None):Use(float),
                            Optional('thickness'): Use(float)},
              'Material': {Optional('muRfunc_filename', default=None): str,
                           'muRfunc': str,
@@ -155,19 +156,26 @@ def tobool(value):
                           Optional('Omega_b', default=None): Use(float),
                           Optional('Omega_c', default=None): Use(float),
                           Optional('kappa', default=None): Use(float),
+                          Optional('kappa_x', default=None): Use(float),
+                          Optional('kappa_y', default=None): Use(float),
                           Optional('B', default=None): Use(float),
+                          Optional('Bx', default=None): Use(float),
+                          Optional('By', default=None): Use(float),
                           Optional('EvdW', default=None): Use(float),
                           'rho_s': Use(float),
                           'D': Use(float),
+                          Optional('D_surf', default=0.): Use(float),
                           Optional('Dfunc_filename', default=None): str,
                           'Dfunc': str,
                           Optional('E_D', default=0.): Use(float),
+                          Optional('E_D_surf', default=0.): Use(float),
                           Optional('dgammadc', default=None): Use(float),
                           Optional('cwet', default=None): Use(float)},
              'Reactions': {Optional('rxnType_filename', default=None): str,
                            'rxnType': str,
                            'k0': Use(float),
                            Optional('E_A', default=0.): Use(float),
+                           Optional('surface_diffusion', default=None): Use(tobool),
                            'alpha': Use(float),
                            Optional('lambda', default=None): Use(float),
                            'Rfilm': Use(float)}}