Repository: TRI-AMDD/mpet (language: python)

Issue:
merging fix of temperature branch into new Delft branch



Discussion / hints:

Fix patch:
diff --git a/bin/MPET_automatic.py b/bin/MPET_automatic.py
new file mode 100644
index 00000000..546dff25
--- /dev/null
+++ b/bin/MPET_automatic.py
@@ -0,0 +1,33 @@
+import os
+import subprocess
+import numpy as np
+
+
+def run_MPET(cwd, config):
+    os.chdir(cwd)
+    subprocess.call(["python", cwd + r"\bin\mpetrun.py", config])
+
+
+cwd = os.getcwd()
+os.chdir(cwd)
+
+
+os.chdir(r".\configs\mark_diff_Crate")
+config_list = os.listdir()
+# print(config_list)
+config_I_want = np.array([])
+for config in config_list:
+    if config[:22] == 'params_system_NMC_Mark':
+        config_I_want = np.append(config_I_want,config)
+print(config_I_want)
+# conf = open('ensemble_parallel_configs.txt',"r")
+# lines = conf.readlines()
+os.chdir(cwd)
+# # print(cwd)
+for line in config_I_want:
+    cwd = os.getcwd()
+    run_MPET(cwd, (cwd + r'\\' + r'configs\\mark_diff_Crate\\' + line))
+
+
+
+            
diff --git a/bin/create_ensamble.py b/bin/create_ensamble.py
new file mode 100644
index 00000000..c9a1707e
--- /dev/null
+++ b/bin/create_ensamble.py
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+
+import configparser
+import sys
+import itertools
+import os
+
+# Values that need ensamble
+# ensamble = [
+#     [("Particles","mean_c"), ["100e-9","300e-9"]],
+#     [("Particles","stddev_c"), ["10e-9","100e-9"]],
+#     [("Sim Params","relTol"), ["1e-5","1e-7"]],
+#     [("Sim Params","absTol"), ["1e-5","1e-7"]],
+#     [("Electrodes", "cathode"), ['D=1e-18-k0=2.cfg','D=1e-18-k0=2.5.cfg','D=1e-18-k0=3.cfg',
+#                                  'D=5e-18-k0=2.cfg','D=5e-18-k0=2.5.cfg','D=5e-18-k0=3.cfg',
+#                                  'D=50e-18-k0=2.cfg','D=50e-18-k0=2.5.cfg','D=50e-18-k0=3.cfg',
+#                                  'D=500e-18-k0=2.cfg','D=50e-18-k0=2.5.cfg','D=500e-18-k0=3.cfg']]
+# ]
+
+# ensamble = [
+#     [("Material","D"), ["1e-18","5e-18","50e-18","500e-18"]],
+#     [("Reactions","k0"), ["2","2.5","3"]],
+# ]
+
+
+def ensemble_definitions():
+    # Values that need ensamble
+    ensamble = [[("Particles","mean_c"), ["100e-9","300e-9"]],
+                [("Particles","stddev_c"), ["10e-9","100e-9"]],
+                [("Sim Params","relTol"), ["1e-5","1e-7"]],
+                [("Sim Params","absTol"), ["1e-5","1e-7"]],
+                [("Electrodes", "cathode"), ['D=1e-18-k0=2.cfg','D=1e-18-k0=2.5.cfg','D=1e-18-k0=3.cfg',
+                                             'D=5e-18-k0=2.cfg','D=5e-18-k0=2.5.cfg','D=5e-18-k0=3.cfg',
+                                             'D=50e-18-k0=2.cfg','D=50e-18-k0=2.5.cfg','D=50e-18-k0=3.cfg',
+                                             'D=500e-18-k0=2.cfg','D=50e-18-k0=2.5.cfg','D=500e-18-k0=3.cfg']]
+                ]
+
+    # helpers
+    keys = [vals[0] for vals in ensamble]
+    val = [vals[1] for vals in ensamble]
+    return keys, val
+
+
+def create_ensemble(cff, keys=None, val=None):
+    with open('ensemble_parallel_configs.txt', "w") as ff:
+        if keys is None and val is None:
+            keys, val = ensemble_definitions()
+        cfg = configparser.ConfigParser()
+        cfg.optionxform = str
+        cfg.read(cff)
+        # Create all variations
+        combinations = list(itertools.product(*val))
+        for combination in combinations:
+            params = dict(zip(keys,combination))
+            new_cfg = cfg
+            nicename = []
+            for key, val in params.items():
+                new_cfg[key[0]][key[1]] = val
+                nicename.append(key[1] + "=" + val)
+
+            # Write config
+            cfg_dir = os.path.dirname(cff)
+            with open(cfg_dir + "/" + "-".join(nicename) + ".cfg", "w") as f:
+                new_cfg.write(f)
+                ff.write(str(cfg_dir + "/" + "-".join(nicename) + ".cfg\n"))
+    return
+
+
+if __name__ == '__main__':
+    # Read in file
+    if len(sys.argv) < 2:
+        print("need the config file [python create_enamble.py <baseconfig>]")
+        exit(1)
+    create_ensemble(sys.argv[1])
\ No newline at end of file
diff --git a/configs/params_LFP_CHR.cfg b/configs/params_LFP_CHR.cfg
new file mode 100644
index 00000000..f27c70a2
--- /dev/null
+++ b/configs/params_LFP_CHR.cfg
@@ -0,0 +1,31 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = CHR
+discretization = 1e-9
+shape = sphere
+thickness = 20e-9
+
+[Material]
+muRfunc = LiFePO4
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+Omega_a = 1.8560e-20
+kappa = 5.0148e-10
+B = 0.01e9
+rho_s = 1.3793e28
+D = 5e-14
+Dfunc = lattice
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = CIET
+k0 = 60
+E_A = 13000
+alpha = 0.5
+# Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0
diff --git a/configs/params_LFP_CHR_meta.cfg b/configs/params_LFP_CHR_meta.cfg
new file mode 100644
index 00000000..3c36c73c
--- /dev/null
+++ b/configs/params_LFP_CHR_meta.cfg
@@ -0,0 +1,31 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = CHR
+discretization = 1.1e-9
+shape = sphere
+thickness = 20e-9
+
+[Material]
+muRfunc = LiFePO4_meta
+noise = true
+noise_prefac = 1e-5
+numnoise = 200
+Omega_a = 2.931e-20
+kappa = 9.9193e-10
+B = 0.05e9
+rho_s = 1.3793e28
+D = 5e-16
+Dfunc = lattice
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = CIET
+k0 = 20
+E_A = 13000
+alpha = 0.5
+# Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0
diff --git a/configs/params_LFP_CHR_plate.cfg b/configs/params_LFP_CHR_plate.cfg
new file mode 100644
index 00000000..914173a0
--- /dev/null
+++ b/configs/params_LFP_CHR_plate.cfg
@@ -0,0 +1,31 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = CHR
+discretization = 1e-9
+shape = plate
+thickness = 100e-9
+
+[Material]
+muRfunc = LiFePO4
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+Omega_a = 1.8560e-20
+kappa = 5.0148e-10
+B = 0.5e9
+rho_s = 1.3793e28
+D = 5e-16
+Dfunc = lattice
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = BV
+k0 = 0.1
+E_A = 13000
+alpha = 0.5
+# Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0
diff --git a/configs/params_LFP_diff_classic.cfg b/configs/params_LFP_diff_classic.cfg
new file mode 100644
index 00000000..770d42f2
--- /dev/null
+++ b/configs/params_LFP_diff_classic.cfg
@@ -0,0 +1,35 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = ACR
+discretization = 1e-9
+shape = C3
+thickness = 50e-9
+
+[Material]
+muRfunc = LiFePO4
+noise = false
+noise_prefac = 1e-5
+numnoise = 200
+Omega_a = 1.8560e-20
+kappa = 5.0148e-10
+B = 0.1916e9
+rho_s = 1.3793e28
+D = 1e-16
+# the effective diff = D_surf/[(thick/0.75)/2]
+Dfunc = lattice
+dgammadc = 2e-29
+cwet = 0.98
+E_D = 1000
+
+[Reactions]
+surface_diffusion = true
+rxnType = BV
+k0 = 0.16
+E_A = 13000
+alpha = 0.5
+# Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0
+
diff --git a/configs/params_LFP_diff_meta.cfg b/configs/params_LFP_diff_meta.cfg
new file mode 100644
index 00000000..61ce2629
--- /dev/null
+++ b/configs/params_LFP_diff_meta.cfg
@@ -0,0 +1,36 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = ACR
+discretization = 1e-9
+shape = C3
+thickness = 100e-9
+
+[Material]
+muRfunc = LiFePO4_meta
+noise = false
+noise_prefac = 1e-5
+numnoise = 200
+Omega_a = 2.931e-20
+kappa = 9.9193e-10
+B = 0.1916e9
+rho_s = 1.3793e28
+D = 1e-15
+# the effective diff = D_surf/[(thick/0.75)/2]
+Dfunc = lattice
+dgammadc = 1e-29
+cwet = 0.98
+E_D = 1000
+
+[Reactions]
+surface_diffusion = true
+# rxnType = BV
+# k0 = 0.16
+rxnType = CIET
+k0 = 30
+E_A = 13000
+alpha = 0.5
+# # Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0
diff --git a/configs/params_LFP_meta.cfg b/configs/params_LFP_meta.cfg
new file mode 100644
index 00000000..891d157f
--- /dev/null
+++ b/configs/params_LFP_meta.cfg
@@ -0,0 +1,31 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = ACR
+discretization = 1e-9
+shape = C3
+thickness = 20e-9
+
+[Material]
+muRfunc = LiFePO4_meta
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+Omega_a = 2.931e-20
+kappa = 9.9193e-10
+B = 0.1916e9
+rho_s = 1.3793e28
+D = 5.3e-19
+Dfunc = lattice
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = BV
+k0 = 1.6e-1
+E_A = 13000
+alpha = 0.5
+# Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0
diff --git a/configs/params_NMC_Chen2020.cfg b/configs/params_NMC_Chen2020.cfg
new file mode 100644
index 00000000..711d9520
--- /dev/null
+++ b/configs/params_NMC_Chen2020.cfg
@@ -0,0 +1,24 @@
+# Default parameters for simulating NMC811 from Chen 2020
+
+[Particles]
+type = diffn
+discretization = 2.e-7
+shape = sphere
+thickness = 20e-9
+
+[Material]
+muRfunc = NMC811_Chen2020
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+rho_s = 3.116e28
+D = 1.8e-13
+Dfunc = constant
+E_D = 5000
+
+[Reactions]
+rxnType = BV_mod01
+k0 = 6.8
+E_A = 17800
+alpha = 0.5
+Rfilm = 0e-0
\ No newline at end of file
diff --git a/configs/params_NaAnode.cfg b/configs/params_NaAnode.cfg
new file mode 100644
index 00000000..6421684f
--- /dev/null
+++ b/configs/params_NaAnode.cfg
@@ -0,0 +1,28 @@
+# Cathode parameters from M. Doyle et al., J. Electrochem. Soc. 143, 1890 (1996).
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = CHR
+discretization = 3e-8
+shape = cylinder
+thickness = 12e-6
+
+[Material]
+muRfunc = NaAnode
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+B = 0.0
+kappa = 4.0e-6
+rho_s = 1.26e28
+D = 2.5e-15
+EvdW = 0.0
+Dfunc = amorph
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = BV_Na
+k0 = 50
+alpha = 0.5
+Rfilm = 0e-0
diff --git a/configs/params_system_NMC_Mark.cfg b/configs/params_system_NMC_Mark.cfg
new file mode 100644
index 00000000..843e1ee7
--- /dev/null
+++ b/configs/params_system_NMC_Mark.cfg
@@ -0,0 +1,210 @@
+# Note: In this file, parameters ending in _c, _s, and _a refer to the
+# cathode (positive electrode), separator, and anode (negative
+# electrode) respectively
+
+[Sim Params]
+# Constant voltage or current or segments of one of them
+# Options: CV, CC, CCsegments, CVsegments
+profileType = CC
+# Battery (dis)charge c-rate (only used for CC), number of capacities / hr
+# (positive for discharge, negative for charge)
+Crate = -0.62
+# Crate = -2.48
+# Optional nominal 1C current density for the cell, A/m^2
+# 1C_current_density = 12.705
+# Voltage cutoffs, V
+Vmax = 4.3
+Vmin = 2.7
+# Battery applied voltage (only used for CV), V
+Vset = 0.12
+# Battery constant power (only used for CP), W/m^2
+power = 1
+# CC/CV segments defining profile for profileType = CCsegments or CVsegments
+# Segments are input as a list of tuples, with the first entry of each
+# row being the setpoint (CC or CV) and the second being the time (in
+# minutes) for which that setpoint is held.
+# 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),
+    ]
+# Continuation directory. If false, begin a fresh simulation with the
+# specified input parameters here. Otherwise, this should be the
+# absolute path to the output directory of the simulation to continue.
+# If continuing, keep the same input parameters as the simulation to
+# continue, changing only `Sim Params` subheadings (other than Nvol
+# and Npart options).
+# Options: false, absolute directory path
+prevDir = false
+# Final time (only used for CV), [s]
+tend = 1.2e3
+# Number disc. in time
+# Note time stepping is adaptive, variable-order stepping, so this
+# affects only the interpolated output values, not simulation
+# accuracy. The output will have all simulation values at a linear
+# spacing between initial and final times with tsteps total outputs.
+tsteps = 200
+# Relative Tolerance
+relTol = 1e-8
+# Absolute Tolerance
+absTol = 1e-8
+# Initial Temperature throughout electrode, K
+T = 298
+# Nonisothermal: true for heat generation throughout electrode, false
+# for no heat generation
+nonisothermal = false
+# Random seed. Set to true to give a random seed in the simulation
+# (affects noise, particle size distribution). Set to true exactly
+# reproducible results -- useful for testing.
+# Options: true, 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 = 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 = 4
+Nvol_s = 10
+Nvol_a = 0
+# Number of particles per volume for cathode and anode
+Npart_c = 8
+Npart_a = 0
+
+[Electrodes]
+# The name of the parameter file describing the cathode particles
+cathode = params_NMC_Chen2020.cfg
+# The name of the parameter file describing the anode particles
+# Used only if Nvol_a > 0
+anode = params_graphite_1param.cfg
+# Rate constant of the Li foil electrode, A/m^2
+# Used only if Nvol_a = 0
+k0_foil = 1
+# Film resistance on the Li foil, Ohm m^2
+Rfilm_foil = 1e-3
+
+[Particles]
+# electrode particle size distribution info, m
+#   C3 -- size along [100] direction
+#   sphere or cylinder -- radius
+# If using stddev = 0, set Npart = 1 for that electrode to avoid
+# wasting computational effort.
+mean_c = 11.5e-6
+stddev_c = 4.5e-6
+mean_a = 100e-9
+stddev_a = 1e-9
+# Use a specific set of particle sizes for the distribution
+# If false, a randomly generated PSD is used. Otherwise the input should
+# be a 2D list of particle radii with 'Npart' rows and 'Nvol' columns.
+specified_psd_c = False
+specified_psd_a = False
+# Initial electrode filling fractions
+# (for disch, anode starts full, cathode starts empty)
+cs0_c = 0.87
+cs0_a = 0.99
+
+[Conductivity]
+# Simulate bulk cathode conductivity (Ohm's Law)?
+# Options: true, false
+simBulkCond_c = false
+simBulkCond_a = false
+# Dimensional conductivity (used if simBulkCond = true), S/m
+sigma_s_c = 1e-1
+sigma_s_a = 1e-1
+# Simulate particle connectivity losses (Ohm's Law)?
+# Options: true, false
+simPartCond_c = false
+simPartCond_a = false
+# Conductance between particles, S = 1/Ohm
+G_mean_c = 1e-14
+G_stddev_c = 0
+G_mean_a = 1e-14
+G_stddev_a = 0
+
+[Thermal Parameters]
+# Heat capacity of anode, cathode, electrolyte material in J/(kg*K).
+cp_c = 700
+cp_s = 700
+cp_a = 700
+# Mass density of anode, cathode, and electrolyte materialin kg/m^3.
+rhom_c = 2500
+rhom_s = 1100
+rhom_a = 2500
+# Heat transfer coefficient with the cell separator (W/(m^2*K))
+h_h = 2e4
+# Thermal conductivity in battery cell (only used for dilute electrolyte model), (W/(m*K)).
+# For Stefan-Maxwell concentrated electrolyte, input from props_elyte.py
+k_h_c = 2.1
+k_h_a = 1.7
+k_h_s = 0.16
+# Includes temperature dependence of entropic heat generation if true, does not include
+# if false
+entropy_heat_gen = False
+
+[Geometry]
+# Thicknesses, m
+L_c = 35e-6
+L_a = 50e-6
+L_s = 25e-6
+# Volume loading percents of active material (volume fraction of solid
+# that is active material)
+P_L_c = 0.82
+P_L_a = 0.69
+# Porosities (liquid volume fraction in each region)
+poros_c = 0.5
+poros_a = 0.1
+poros_s = 0.55
+# Bruggeman exponent (tortuosity = porosity^bruggExp)
+BruggExp_c = -0.5
+BruggExp_a = -0.5
+BruggExp_s = -0.5
+
+[Electrolyte]
+# Initial electrolyte conc., mol/m^3
+c0 = 1000
+# Cation/anion charge number (e.g. 2, -1 for CaCl_2)
+zp = 1
+zm = -1
+# Cation/anion dissociation number (e.g. 1, 2 for CaCl_2)
+nup = 1
+num = 1
+# Electrolyte model,
+# Options: dilute, SM
+#   dilute: assume dilute, binary electrolyte model; phi in
+#   electrolyte is an "inner potential" or an "quasi-electrostatic
+#   potential"
+#   SM: use Stefan-Maxwell model for electrolyte transport; phi in
+#   electrolyte is that measured with a Li metal reference electrode
+#   relative to some fixed position in solution.
+# WARNING: Using SM model with BV reaction models for electrode
+# particles assumes individual ion activities are given by their
+# concentrations for the reaction rate exchange current density.
+elyteModelType = SM
+# Stefan-Maxwell property set, see props_elyte.py file
+# Options:
+#   test1: parameter set for testing
+#   LiClO4_PC: electrolyte/solvent used in Fuller, Doyle, Newman 1994
+#              conductivity taken from dualfoil5.2.f
+#   valoen_bernardi: LiPF6 in carbonates as in Bernardi and Go 2011
+#SMset_filename = filename.py   # optional, to load SMset from custom file instead of props_elyte.py
+SMset = valoen_bernardi
+# Reference electrode (defining the electrolyte potential) information:
+# number of electrons transfered in the reaction, 1 for Li/Li+
+n = 1
+# Stoichiometric coefficient of cation, -1 for Li/Li+
+sp = -1
+# Dilute solution properties (used only if elyteModelType = "dilute")
+# Cation/anion diff, m^2/s
+# e.g. for LiPF6 in EC/DMC, Dp = 2.2e-10, Dm = 2.94e-10
+Dp = 2.2e-10
+Dm = 2.94e-10
diff --git a/configs/params_system_Na.cfg b/configs/params_system_Na.cfg
new file mode 100644
index 00000000..9338ad9c
--- /dev/null
+++ b/configs/params_system_Na.cfg
@@ -0,0 +1,215 @@
+# Note: In this file, parameters ending in _c, _s, and _a refer to the
+# cathode (positive electrode), separator, and anode (negative
+# electrode) respectively
+
+# ask: thickness of anode and cathode
+# mass loading
+# radius or diameter ? 
+
+
+
+[Sim Params]
+# Constant voltage or current or segments of one of them
+# Options: CV, CC, CCsegments, CVsegments
+profileType = CC
+# Battery (dis)charge c-rate (only used for CC), number of capacities / hr
+# (positive for discharge, negative for charge)
+Crate = 5
+# Optional nominal 1C current density for the cell, A/m^2
+# 1C_current_density = 12.705
+# Voltage cutoffs, V
+Vmax = 2
+Vmin = -0.25
+# Battery applied voltage (only used for CV), V
+Vset = 0.12
+# Battery constant power (only used for CP), W/m^2
+power = 1
+# CC/CV segments defining profile for profileType = CCsegments or CVsegments
+# Segments are input as a list of tuples, with the first entry of each
+# row being the setpoint (CC or CV) and the second being the time (in
+# minutes) for which that setpoint is held.
+# 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),
+    ]
+# Continuation directory. If false, begin a fresh simulation with the
+# specified input parameters here. Otherwise, this should be the
+# absolute path to the output directory of the simulation to continue.
+# If continuing, keep the same input parameters as the simulation to
+# continue, changing only `Sim Params` subheadings (other than Nvol
+# and Npart options).
+# Options: false, absolute directory path
+prevDir = false
+# Final time (only used for CV), [s]
+tend = 1.2e3
+# Number disc. in time
+# Note time stepping is adaptive, variable-order stepping, so this
+# affects only the interpolated output values, not simulation
+# accuracy. The output will have all simulation values at a linear
+# spacing between initial and final times with tsteps total outputs.
+tsteps = 200
+# Relative Tolerance
+relTol = 1e-12
+# Absolute Tolerance
+absTol = 1e-12
+# Initial Temperature throughout electrode, K
+T = 298
+# Nonisothermal: true for heat generation throughout electrode, false
+# for no heat generation
+nonisothermal = false
+# Random seed. Set to true to give a random seed in the simulation
+# (affects noise, particle size distribution). Set to true exactly
+# reproducible results -- useful for testing.
+# Options: true, 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 = 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 = 0
+Nvol_a = 0
+# Number of particles per volume for cathode and anode
+Npart_c = 10
+Npart_a = 0
+
+[Electrodes]
+# The name of the parameter file describing the cathode particles
+cathode = params_LiC6_LIONSIMBA.cfg
+# The name of the parameter file describing the anode particles
+# Used only if Nvol_a > 0
+anode = params_graphite_1param.cfg
+# Rate constant of the Li foil electrode, A/m^2
+# Used only if Nvol_a = 0
+k0_foil = 1e0
+# Film resistance on the Li foil, Ohm m^2
+Rfilm_foil = 0e-0
+
+[Particles]
+# electrode particle size distribution info, m
+#   C3 -- size along [100] direction
+#   sphere or cylinder -- radius
+# If using stddev = 0, set Npart = 1 for that electrode to avoid
+# wasting computational effort.
+mean_c = 2.5e-6
+stddev_c = 10e-7
+mean_a = 100e-9
+stddev_a = 1e-9
+# Use a specific set of particle sizes for the distribution
+# If false, a randomly generated PSD is used. Otherwise the input should
+# be a 2D list of particle radii with 'Npart' rows and 'Nvol' columns.
+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_a = 0.99
+
+[Conductivity]
+# Simulate bulk cathode conductivity (Ohm's Law)?
+# Options: true, false
+simBulkCond_c = false
+simBulkCond_a = false
+# Dimensional conductivity (used if simBulkCond = true), S/m
+sigma_s_c = 1e-1
+sigma_s_a = 1e-1
+# Simulate particle connectivity losses (Ohm's Law)?
+# Options: true, false
+simPartCond_c = false
+simPartCond_a = false
+# Conductance between particles, S = 1/Ohm
+G_mean_c = 1e-14
+G_stddev_c = 0
+G_mean_a = 1e-14
+G_stddev_a = 0
+
+[Thermal Parameters]
+# Heat capacity of anode, cathode, electrolyte material in J/(kg*K).
+cp_c = 700
+cp_s = 700
+cp_a = 700
+# Mass density of anode, cathode, and electrolyte materialin kg/m^3.
+rhom_c = 2500
+rhom_s = 1100
+rhom_a = 2500
+# Heat transfer coefficient with the cell separator (W/(m^2*K))
+h_h = 2e4
+# Thermal conductivity in battery cell (only used for dilute electrolyte model), (W/(m*K)).
+# For Stefan-Maxwell concentrated electrolyte, input from props_elyte.py
+k_h_c = 2.1
+k_h_a = 1.7
+k_h_s = 0.16
+# Includes temperature dependence of entropic heat generation if true, does not include
+# if false
+entropy_heat_gen = False
+
+[Geometry]
+# Thicknesses, m
+L_c = 80e-6
+L_a = 50e-6
+L_s = 25e-6
+# Volume loading percents of active material (volume fraction of solid
+# that is active material)
+P_L_c = 0.75
+P_L_a = 0.69
+# Porosities (liquid volume fraction in each region)
+poros_c = 0.6
+poros_a = 0.4
+poros_s = 1.0
+# Bruggeman exponent (tortuosity = porosity^bruggExp)
+BruggExp_c = -0.4
+BruggExp_a = -0.5
+BruggExp_s = -0.5
+
+[Electrolyte]
+# Initial electrolyte conc., mol/m^3
+c0 = 1000
+# Cation/anion charge number (e.g. 2, -1 for CaCl_2)
+zp = 1
+zm = -1
+# Cation/anion dissociation number (e.g. 1, 2 for CaCl_2)
+nup = 1
+num = 1
+# Electrolyte model,
+# Options: dilute, SM
+#   dilute: assume dilute, binary electrolyte model; phi in
+#   electrolyte is an "inner potential" or an "quasi-electrostatic
+#   potential"
+#   SM: use Stefan-Maxwell model for electrolyte transport; phi in
+#   electrolyte is that measured with a Li metal reference electrode
+#   relative to some fixed position in solution.
+# WARNING: Using SM model with BV reaction models for electrode
+# particles assumes individual ion activities are given by their
+# concentrations for the reaction rate exchange current density.
+elyteModelType = dilute
+# Stefan-Maxwell property set, see props_elyte.py file
+# Options:
+#   test1: parameter set for testing
+#   LiClO4_PC: electrolyte/solvent used in Fuller, Doyle, Newman 1994
+#              conductivity taken from dualfoil5.2.f
+#   valoen_bernardi: LiPF6 in carbonates as in Bernardi and Go 2011
+#SMset_filename = filename.py   # optional, to load SMset from custom file instead of props_elyte.py
+SMset = valoen_bernardi
+# Reference electrode (defining the electrolyte potential) information:
+# number of electrons transfered in the reaction, 1 for Li/Li+
+n = 1
+# Stoichiometric coefficient of cation, -1 for Li/Li+
+sp = -1
+# Dilute solution properties (used only if elyteModelType = "dilute")
+# Cation/anion diff, m^2/s
+# e.g. for LiPF6 in EC/DMC, Dp = 2.2e-10, Dm = 2.94e-10
+Dp = 2e-10
+Dm = 3e-10
diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py
index 221312dc..e723b6cf 100644
--- a/mpet/config/configuration.py
+++ b/mpet/config/configuration.py
@@ -711,6 +711,9 @@ def _distr_part(self):
             elif solidShape == 'cylinder':
                 psd_area = 2 * np.pi * psd_len * self[trode, 'thickness']
                 psd_vol = np.pi * psd_len**2 * self[trode, 'thickness']
+            elif solidShape == 'plate':
+                psd_area = 2 * (psd_len/psd_len) * self[trode, 'thickness']**2
+                psd_vol = psd_len * self[trode, 'thickness']**2
             else:
                 raise NotImplementedError(f'Unknown solid shape: {solidShape}')
 
@@ -839,8 +842,9 @@ def _verify_config(self):
             solidType = self[trode, 'type']
             if solidType in ["ACR", "homog_sdn"] and solidShape != "C3":
                 raise Exception("ACR and homog_sdn req. C3 shape")
-            if (solidType in ["CHR", "diffn"] and solidShape not in ["sphere", "cylinder"]):
-                raise NotImplementedError("CHR and diffn req. sphere or cylinder")
+            if (solidType in ["CHR", "diffn"] and solidShape not in ["sphere", "cylinder",
+                                                                     "plate"]):
+                raise NotImplementedError("CHR and diffn req. sphere, cylinder or plate shape")
 
     @staticmethod
     def size2regsln(size):
diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py
index b3f95609..89fb94a9 100644
--- a/mpet/config/schemas.py
+++ b/mpet/config/schemas.py
@@ -144,8 +144,8 @@ def tobool(value):
                                                                   + constants.two_var_types),
                            'discretization': Use(float),
                            'shape': lambda x:
-                               check_allowed_values(x, ["C3", "sphere", "cylinder", "homog_sdn"]),
-                           Optional('discretization_ver', default=None):Use(float),
+                               check_allowed_values(x, ["C3", "sphere", "cylinder",
+                                                        "plate", "homog_sdn"]),
                            Optional('thickness'): Use(float)},
              'Material': {Optional('muRfunc_filename', default=None): str,
                           'muRfunc': str,
@@ -172,6 +172,7 @@ def tobool(value):
                           Optional('dgammadc', default=None): Use(float),
                           Optional('cwet', default=None): Use(float)},
              'Reactions': {Optional('rxnType_filename', default=None): str,
+                           Optional('surface_diffusion', default=None): Use(tobool),
                            'rxnType': str,
                            'k0': Use(float),
                            Optional('E_A', default=0.): Use(float),
diff --git a/mpet/electrode/diffusion/amorph.py b/mpet/electrode/diffusion/amorph.py
new file mode 100644
index 00000000..3f6990d5
--- /dev/null
+++ b/mpet/electrode/diffusion/amorph.py
@@ -0,0 +1,2 @@
+def amorph(y):
+    return 1+y
diff --git a/mpet/electrode/materials/LiFePO4_meta.py b/mpet/electrode/materials/LiFePO4_meta.py
new file mode 100644
index 00000000..eb53f4da
--- /dev/null
+++ b/mpet/electrode/materials/LiFePO4_meta.py
@@ -0,0 +1,20 @@
+import numpy as np
+
+
+def LiFePO4_meta(self, y, ybar, T, muR_ref):
+    """ Bai, Cogswell, Bazant 2011 """
+    muRtheta = -self.eokT*3.422
+    muRhomog = 2*self.reg_sln(y, T, self.get_trode_param("Omega_a")/2)
+    b = 0.3
+    pi = 3.14159
+    a = -0.1
+    phi = -700
+    muR_meta_1 = a*20*pi*np.cos(pi*y-b)*(np.sin(pi*y-b))**19
+    muR_meta_2 = phi*6*((y*(1-y))**5)*(1-2*y)
+    # a = -0.1
+    # muR_meta = a*4*pi*np.cos(pi*y-b)*(np.sin(pi*y-b))**3
+    muRnonHomog = self.general_non_homog(y, ybar)
+    muR = muRhomog + muRnonHomog + muR_meta_1 + muR_meta_2
+    actR = np.exp(muR/T)
+    muR += muRtheta + muR_ref
+    return muR, actR
diff --git a/mpet/electrode/materials/NMC811_Chen2020.py b/mpet/electrode/materials/NMC811_Chen2020.py
new file mode 100644
index 00000000..851b9a24
--- /dev/null
+++ b/mpet/electrode/materials/NMC811_Chen2020.py
@@ -0,0 +1,13 @@
+import numpy as np
+
+
+def NMC811_Chen2020(self, y, ybar, T, muR_ref):
+    """
+    This function was obtained from Chen et al. 2020.
+    """
+    OCV = (-0.8090*y + 4.4875 - 0.0428*np.tanh(18.5138*(y - 0.5542)) 
+           - 17.7326*np.tanh(15.7890*(y - 0.3117)) 
+           + 17.5842*np.tanh(15.9308*(y - 0.3120)))
+    muR = self.get_muR_from_OCV(OCV, muR_ref)
+    actR = None
+    return muR, actR
diff --git a/mpet/electrode/materials/NaAnode.py b/mpet/electrode/materials/NaAnode.py
new file mode 100644
index 00000000..89a449b0
--- /dev/null
+++ b/mpet/electrode/materials/NaAnode.py
@@ -0,0 +1,28 @@
+import numpy as np
+
+#     muRtheta = -self.eokT*1.3356
+#     OCV = (-13.522*y**6 + 42.451*y**5 - 46.305*y**4
+#            + 15.214*y**3 + 7.0034*y**2 - 6.0744*y)
+
+
+def NaAnode(self, y, ybar, T, muR_ref):
+
+    muRtheta = -self.eokT*1.435080161287005
+    muRtheta = -self.eokT*1.45
+
+    OCV = (-7.6991314610297525*y
+           + 16.406426373526156*y**2 + -13.120937152671123*y**3
+           + -8.140487082887885*y**4 + 22.754536573994052*y**5
+           + -11.76391247146248*y**6)
+    
+#     muRtheta = -self.eokT*1.2356
+#     OCV = (-13.522*y**6 + 42.451*y**5 - 46.305*y**4
+#            + 15.214*y**3 + 7.0034*y**2 - 6.0744*y)
+
+    muRhomog = -self.eokT*OCV
+    muRnonHomog = self.general_non_homog(y, ybar)
+    muR = muRhomog + muRnonHomog
+    actR = np.exp(muR/T)
+    muR += muRtheta + muR_ref
+
+    return muR, actR
diff --git a/mpet/electrode/reactions/BV_Na.py b/mpet/electrode/reactions/BV_Na.py
new file mode 100644
index 00000000..19e1e2a4
--- /dev/null
+++ b/mpet/electrode/reactions/BV_Na.py
@@ -0,0 +1,12 @@
+import numpy as np
+
+
+def BV_Na(eta, c_sld, c_lyte, k0, E_A, T, act_R=None,
+          act_lyte=None, lmbda=None, alpha=None):
+    if act_R is None:
+        act_R = c_sld/(1-c_sld)
+    gamma_ts = (1./(1-c_sld)**0.5)
+    ecd = (k0 * act_lyte**(1-alpha)
+           * act_R**(alpha) / gamma_ts)
+    Rate = ecd * np.exp(-E_A/T + E_A/1) * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T))
+    return Rate
diff --git a/mpet/geometry.py b/mpet/geometry.py
index ed81f2e4..f02ba980 100644
--- a/mpet/geometry.py
+++ b/mpet/geometry.py
@@ -34,6 +34,16 @@ def get_unit_solid_discr(Shape, N):
         vol_vec[-1] = np.pi * h * (Rs * dr - dr**2 / 4.)
         Vp = np.pi * Rs**2 * h
         volfrac_vec = vol_vec / Vp
+    elif Shape == "plate":
+        L = 1.
+        h = 1.
+        dx = h / (N - 1)
+        r_vec = np.linspace(0, h, N)
+        vol_vec = L**2 * dx * np.ones(N)
+        vol_vec[0] = dx/2.
+        vol_vec[-1] = dx/2.
+        Vp = L**2 * h
+        volfrac_vec = vol_vec / Vp
     else:
         raise NotImplementedError("Fix shape volumes!")
     return r_vec, volfrac_vec
@@ -68,8 +78,13 @@ def calc_curv(c, dr, r_vec, Rs, beta_s, particleShape):
         curv[N-1] = (
             (1./Rs)*beta_s
             + (2*c[-2] - 2*c[-1] + 2*dr*beta_s)/dr**2)
+    elif particleShape == "plate":
+        dx = dr
+        curv[0] = (2*c[1] - 2*c[0]) / dx**2
+        curv[1:N-1] = np.diff(c,2)/dx**2
+        curv[N-1] = (2*c[-2] - 2*c[-1] + 2*dr*beta_s)/dx**2
     else:
-        raise NotImplementedError("calc_curv_c only for sphere and cylinder")
+        raise NotImplementedError("calc_curv_c only for sphere, cylinder or plate")
     return curv
 
 
diff --git a/mpet/main.py b/mpet/main.py
index 3ce8c68d..9a4542e3 100644
--- a/mpet/main.py
+++ b/mpet/main.py
@@ -94,6 +94,9 @@ def main(paramfile, keepArchive=True):
 
     # Directories we'll store output in.
     outdir_name = time.strftime("%Y%m%d_%H%M%S", time.localtime())
+    config_file = os.path.basename(paramfile)
+    config_base = os.path.splitext(config_file)[0]
+    outdir_name = "_".join((time.strftime("%Y%m%d_%H%M%S", time.localtime()), config_base))
     outdir_path = os.path.join(os.getcwd(), "history")
     outdir = os.path.join(outdir_path, outdir_name)
     # Make sure there's a place to store the output
diff --git a/mpet/mod_battery.py b/mpet/mod_battery.py
new file mode 100644
index 00000000..53e9388c
--- /dev/null
+++ b/mpet/mod_battery.py
@@ -0,0 +1,40 @@
+# import daetools.pyDAE as dae
+# from pyUnits import s
+
+# import numpy as np
+
+# import mpet.extern_funcs as extern_funcs
+# import mpet.geometry as geom
+# import mpet.mod_electrodes as mod_electrodes
+# import mpet.ports as ports
+# import mpet.utils as utils
+# from mpet.config import constants
+# from mpet.daeVariableTypes import mole_frac_t, elec_pot_t, conc_t, temp_t
+
+# # Dictionary of end conditions
+# endConditions = {
+#     1:"Vmax reached",
+#     2:"Vmin reached"}
+
+# class ModCell(dae.daeModel):
+#     def __init__(self, config, Name, Parent=None, Description=""):
+#         dae.daeModel.__init__(self, Name, Parent, Description)
+
+#         self.config = config
+#         Nlayers = config["Nlayers"]
+
+#         self.DmnCell = dae.daeDomain(
+#             "DmnCell", self, dae.unit(),
+#             "Simulated layers in battery")
+               
+#         # Variables
+#         self.phi_bound = {}
+#         self.T_bound = {}
+
+#         # ports
+#         self.portsOutT = {}
+
+#          for lInd in range(Nlayers):
+#            self.portsOutLyte[lInd] = ports.portFromCell(
+#                     "portCell{lInd}".format(lInd=lInd), dae.eOutletPort,
+#                     self, "Electrolyte port to particles")
\ No newline at end of file
diff --git a/mpet/mod_cell.py b/mpet/mod_cell.py
index bd5f487b..23b26f93 100644
--- a/mpet/mod_cell.py
+++ b/mpet/mod_cell.py
@@ -215,7 +215,7 @@ def DeclareEquations(self):
                              * self.particles[trode][vInd,pInd].dcbardt())
                 eq.Residual = self.R_Vp[trode](vInd) - RHS
 
-        # Define dimensionless R_Vp for each electrode volume
+        # Define dimensionless Q_Vp for each electrode volume
         for trode in trodes:
             for vInd in range(Nvol[trode]):
                 eq = self.CreateEquation(
@@ -589,8 +589,8 @@ def get_lyte_internal_fluxes(c_lyte, phi_lyte, T_lyte, disc, config, Nvol):
                             + (1./(num*zm)*(1-tp0(c_edges_int, T_edges_int))*i_edges_int))
     q_edges_int = -k_h*np.diff(T_lyte)/dxd1
     # replace boundary conditions since they are ghost points with convective BCs
-    q_edges_int[0] = config["h_h"]*(1 - T_lyte[1])
-    q_edges_int[-1] = config["h_h"]*(T_lyte[-2] - 1)
+    q_edges_int[0] = config["h_h"]*(config["T"] - T_lyte[1])
+    q_edges_int[-1] = config["h_h"]*(T_lyte[-2] - config["T"])
     return Nm_edges_int, i_edges_int, q_edges_int
 
 
diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py
index f0cf445f..7e4a14b1 100644
--- a/mpet/mod_electrodes.py
+++ b/mpet/mod_electrodes.py
@@ -636,6 +636,11 @@ def DeclareEquations(self):
             eq.Residual -= self.c.dt(k) * volfrac_vec[k]
 
         c = np.empty(N, dtype=object)
+        c_surf = np.empty(N, dtype=object)
+        # if self.get_trode_param("surface_diffusion"):
+        #     ctmp = utils.get_var_vec(self.c,N)
+        #     c = np.hstack((self.c_left_GP(),ctmp,self.c_right_GP()))
+        # else:
         c[:] = [self.c(k) for k in range(N)]
         if self.get_trode_param("type") in ["ACR", "diffn", "CHR"]:
             # Equations for 1D particles of 1 field varible
@@ -685,6 +690,12 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise):
             c_surf = c
             muR_surf, actR_surf = calc_muR(
                 c_surf, self.cbar(), self.T_lyte(), self.config, self.trode, self.ind)
+            if noise is None:
+                muR_surf = muR_surf
+            else:
+                muR_surf = muR_surf + noise(dae.Time().Value)
+                actR_surf = np.exp(muR_surf/self.T_lyte())
+
         elif self.get_trode_param("type") in ["diffn", "CHR"]:
             muR, actR = calc_muR(c, self.cbar(), self.T_lyte(),
                                  self.config, self.trode, self.ind)
@@ -735,21 +746,44 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise):
                 area_vec = 4*np.pi*edges**2
             elif self.get_trode_param("shape") == "cylinder":
                 area_vec = 2*np.pi*edges  # per unit height
+            elif self.get_trode_param("shape") == "plate":
+                area_vec = np.ones(np.shape(edges))
             RHS = -np.diff(Flux_vec * area_vec)
 
         dcdt_vec = np.empty(N, dtype=object)
         dcdt_vec[0:N] = [self.c.dt(k) for k in range(N)]
         LHS_vec = MX(Mmat, dcdt_vec)
-        for k in range(N):
-            eq = self.CreateEquation("dcsdt_discr{k}".format(k=k))
-            eq.Residual = LHS_vec[k] - RHS[k]
+        if self.get_trode_param("surface_diffusion"):
+            surf_diff_vec = self.calc_surf_diff(c_surf, muR_surf)
+            for k in range(N):
+                eq = self.CreateEquation("dcsdt_discr{k}".format(k=k))
+                eq.Residual = LHS_vec[k] - RHS[k] - surf_diff_vec[k]
+        else:
+            for k in range(N):
+                eq = self.CreateEquation("dcsdt_discr{k}".format(k=k))
+                eq.Residual = LHS_vec[k] - RHS[k]
 
         if self.get_trode_param("type") in ["ACR"]:
             return eta[-1], c_surf[-1]
         else:
             return eta, c_surf
 
-
+    def calc_surf_diff(self,c_surf, muR_surf):
+        # the surface diffusion keeps the volume centered method of the ACR
+        N = np.size(c_surf)
+        dxs = 1./N
+        D_eff = (self.get_trode_param("D")/(self.T_lyte())
+                 * np.exp(-self.get_trode_param("E_D")/(self.T_lyte())
+                          + self.get_trode_param("E_D")/1))
+        c_edges = utils.mean_linear(c_surf)
+        surf_flux = np.empty(N+1, dtype=object)
+        surf_flux[1:-1] = -D_eff*c_edges*(1-c_edges)*np.diff(muR_surf)/(dxs)
+        surf_flux[0] = 0.
+        surf_flux[-1] = 0.
+        surf_diff = -(np.diff(surf_flux))/(dxs)
+        return surf_diff
+
+   
 def calc_eta(muR, muO):
     return muR - muO
 
@@ -758,7 +792,7 @@ def get_Mmat(shape, N):
     r_vec, volfrac_vec = geo.get_unit_solid_discr(shape, N)
     if shape == "ACR":
         Mmat = sprs.eye(N, N, format="csr")
-    elif shape in ["sphere", "cylinder","C3"]:
+    elif shape in ["sphere", "cylinder", 'plate']:
         Rs = 1.
         # For discretization background, see Zeng & Bazant 2013
         # Mass matrix is common for each shape, diffn or CHR
@@ -766,8 +800,10 @@ def get_Mmat(shape, N):
             Vp = 4./3. * np.pi * Rs**3
         elif shape == "cylinder":
             Vp = np.pi * Rs**2  # per unit height
-        elif shape == "C3":
-            Vp = Rs
+        elif shape == "plate":
+            L = 1.
+            h = 1.
+            Vp = L**2 * h
         vol_vec = Vp * volfrac_vec
         M1 = sprs.diags([1./8, 3./4, 1./8], [-1, 0, 1],
                         shape=(N, N), format="csr")
diff --git a/mpet/plot/plot_data.py b/mpet/plot/plot_data.py
index 65cc3b77..6941d090 100644
--- a/mpet/plot/plot_data.py
+++ b/mpet/plot/plot_data.py
@@ -699,7 +699,7 @@ def animate(tind):
         # Set up colors.
         # Define if you want smooth or discrete color changes
         # Option: "smooth" or "discrete"
-        color_changes = "discrete"
+        color_changes = "smooth"
 #        color_changes = "smooth"
         # Discrete color changes:
         if color_changes == "discrete":
@@ -721,8 +721,9 @@ def animate(tind):
         # Smooth colormap changes:
         if color_changes == "smooth":
             # generated with colormap.org
-            cmaps = np.load("colormaps_custom.npz")
-            cmap_data = cmaps["GnYlRd_3"]
+            cmap_location = os.path.dirname(os.path.abspath(__file__)) + r'\colormaps_custom.npz'
+            cmaps = np.load(cmap_location)
+            cmap_data = cmaps['GnYlRd_1']
             cmap = mpl.colors.ListedColormap(cmap_data/255.)
 
         size_frac_min = 0.10
@@ -775,6 +776,8 @@ def animate(tind):
                 bottom_left = center - size / 2
                 rects[vInd,pInd] = plt.Rectangle(
                     bottom_left, size, size, color=color)
+                # rects[vInd,pInd] = plt.Circle(
+                #     center, radius=size/4, color=color)
             # Create a group of rectange "patches" from the rects array
             collection[indx] = mcollect.PatchCollection(rects.reshape(-1))
             # Put them on the axes
diff --git a/mpet/ports.py b/mpet/ports.py
index 0c1686c6..9decd8a7 100644
--- a/mpet/ports.py
+++ b/mpet/ports.py
@@ -3,6 +3,19 @@
 
 from mpet.daeVariableTypes import mole_frac_t, elec_pot_t
 
+# class portFromCell(dae.daePort):
+#     def __init__(self, Name, PortType, Model, Description=""):
+#         dae.daePort.__init__(self, Name, PortType, Model, Description)
+#         self.c_lyte = dae.daeVariable(
+#             "c_lyte", mole_frac_t, self,
+#             "Concentration in the electrolyte")
+#         self.T_lyte = dae.daeVariable(
+#             "T_lyte", mole_frac_t, self,
+#             "Temperature in the electrolyte")
+#         self.phi_lyte = dae.daeVariable(
+#             "phi_lyte", elec_pot_t, self,
+#             "Electric potential in the electrolyte")
+
 
 class portFromElyte(dae.daePort):
     def __init__(self, Name, PortType, Model, Description=""):
diff --git a/mpet/props_am.py b/mpet/props_am.py
index f7b1cb38..2ad61fc7 100644
--- a/mpet/props_am.py
+++ b/mpet/props_am.py
@@ -161,6 +161,19 @@ def non_homog_round_wetting(self, y, ybar, B, kappa, beta_s, shape, r_vec):
         muR_nh = B*(y - ybar) - kappa*curv
         return muR_nh
 
+    def non_homog_rect_natural(self, y, ybar, B, kappa, beta_s):
+        """ Helper function """
+        # the taylor expansion at the edges is used
+        N = len(y)
+        ytmp = np.empty(N+2, dtype=object)
+        dxs = 1./N
+        ytmp[1:-1] = y
+        ytmp[0] = y[0] + dxs*beta_s*(y[0]*(1-y[0]))*6
+        ytmp[-1] = y[-1] + dxs*beta_s*(y[-1]*(1-y[-1]))*6
+        curv = np.diff(ytmp, 2)/(dxs**2)
+        muR_nh = -kappa*curv + B*(y - ybar)
+        return muR_nh
+
     def general_non_homog(self, y, ybar):
         """ Helper function """
         ptype = self.get_trode_param("type")
@@ -178,14 +191,18 @@ def general_non_homog(self, y, ybar):
             shape = self.get_trode_param("shape")
             kappa = self.get_trode_param("kappa")
             B = self.get_trode_param("B")
-            if shape == "C3":
+            if shape in ["C3"]:
                 if mod1var:
                     cwet = self.get_trode_param("cwet")
-                    muR_nh = self.non_homog_rect_fixed_csurf(
-                        y, ybar, B, kappa, cwet)
+                    if self.get_trode_param("surface_diffusion"):
+                        beta_s = self.get_trode_param("beta_s")
+                        muR_nh = self.non_homog_rect_natural(y, ybar, B, kappa, beta_s)
+                    else:
+                        muR_nh = self.non_homog_rect_fixed_csurf(
+                            y, ybar, B, kappa, cwet)
                 elif mod2var:
                     raise NotImplementedError("no 2param C3 model known")
-            elif shape in ["cylinder", "sphere"]:
+            elif shape in ["cylinder", "sphere", "plate"]:
                 beta_s = self.get_trode_param("beta_s")
                 r_vec = geo.get_unit_solid_discr(shape, N)[0]
                 if mod1var:
diff --git a/mpet/sim.py b/mpet/sim.py
index ff854cca..5a809a2e 100644
--- a/mpet/sim.py
+++ b/mpet/sim.py
@@ -67,9 +67,13 @@ def SetUpVariables(self):
         if not config["prevDir"] or config["prevDir"] == "false":
             # Solids
             for tr in config["trodes"]:
-                cs0 = config['cs0'][tr]
+                cs0_ic = config['cs0'][tr]
                 # Guess initial filling fractions
-                self.m.ffrac[tr].SetInitialGuess(cs0)
+                self.m.ffrac[tr].SetInitialGuess(cs0_ic)
+                epsrnd = 0
+                rndavg = epsrnd*(np.random.rand(Nvol[tr],Npart[tr]) - 0.5)
+                rndavg -= np.mean(rndavg)
+                cs0_mat = cs0_ic*np.ones((Nvol[tr],Npart[tr]), dtype=object) + rndavg
                 for i in range(Nvol[tr]):
                     # Guess initial volumetric reaction rates
                     self.m.R_Vp[tr].SetInitialGuess(i, 0.0)
@@ -84,6 +88,7 @@ def SetUpVariables(self):
                     for j in range(Npart[tr]):
                         Nij = config["psd_num"][tr][i,j]
                         part = self.m.particles[tr][i,j]
+                        cs0 = cs0_mat[i,j]
                         # Guess initial value for the average solid
                         # concentrations and set initial value for
                         # solid concentrations


Test patch:
diff --git a/configs/params_system_LFP_test.cfg b/configs/params_system_LFP_test.cfg
new file mode 100644
index 00000000..cac0e824
--- /dev/null
+++ b/configs/params_system_LFP_test.cfg
@@ -0,0 +1,211 @@
+# Note: In this file, parameters ending in _c, _s, and _a refer to the
+# cathode (positive electrode), separator, and anode (negative
+# electrode) respectively
+
+[Sim Params]
+# Constant voltage or current or segments of one of them
+# Options: CV, CC, CCsegments, CVsegments
+profileType = CCsegments
+# Battery (dis)charge c-rate (only used for CC), number of capacities / hr
+# (positive for discharge, negative for charge)
+Crate = -3
+# Optional nominal 1C current density for the cell, A/m^2
+# 1C_current_density = 12.705
+# Voltage cutoffs, V
+Vmax = 3.8
+Vmin = 3.1
+# Battery applied voltage (only used for CV), V
+Vset = 0.12
+# Battery constant power (only used for CP), W/m^2
+power = 1
+# CC/CV segments defining profile for profileType = CCsegments or CVsegments
+# Segments are input as a list of tuples, with the first entry of each
+# row being the setpoint (CC or CV) and the second being the time (in
+# minutes) for which that setpoint is held.
+# This can have as many rows (segments) as desired.
+# Note: It's okay to leave commented lines within the segments list
+segments = [
+    # (-0.2, 150),
+    (-5,6),
+    (0, 30),
+    (-3, 10),
+    ]
+# Continuation directory. If false, begin a fresh simulation with the
+# specified input parameters here. Otherwise, this should be the
+# absolute path to the output directory of the simulation to continue.
+# If continuing, keep the same input parameters as the simulation to
+# continue, changing only `Sim Params` subheadings (other than Nvol
+# and Npart options).
+# Options: false, absolute directory path
+prevDir = false
+# prevDir = C:\Users\pierfrancescoo\Documents\Phase-field\common-mpet\mpet-dev\history\20230323_191451_params_system_LFP_test
+# Final time (only used for CV), [s]
+tend = 1.2e3
+# Number disc. in time
+# Note time stepping is adaptive, variable-order stepping, so this
+# affects only the interpolated output values, not simulation
+# accuracy. The output will have all simulation values at a linear
+# spacing between initial and final times with tsteps total outputs.
+tsteps = 1000
+# Relative Tolerance
+relTol = 1e-6
+# Absolute Tolerance
+absTol = 1e-6
+# Initial Temperature throughout electrode, K
+T = 298
+# Nonisothermal: true for heat generation throughout electrode, false
+# for no heat generation
+nonisothermal = false
+# Random seed. Set to true to give a random seed in the simulation
+# (affects noise, particle size distribution). Set to true exactly
+# reproducible results -- useful for testing.
+# Options: true, 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 = 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 = 8
+Nvol_s = 0
+Nvol_a = 0
+# Number of particles per volume for cathode and anode
+Npart_c = 8
+Npart_a = 0
+
+[Electrodes]
+# The name of the parameter file describing the cathode particles
+cathode = params_LFP_CHR.cfg
+# The name of the parameter file describing the anode particles
+# Used only if Nvol_a > 0
+anode = params_graphite_1param.cfg
+# Rate constant of the Li foil electrode, A/m^2
+# Used only if Nvol_a = 0
+k0_foil = 10e0
+# Film resistance on the Li foil, Ohm m^2
+Rfilm_foil = 0e-0
+
+[Particles]
+# electrode particle size distribution info, m
+#   C3 -- size along [100] direction
+#   sphere or cylinder -- radius
+# If using stddev = 0, set Npart = 1 for that electrode to avoid
+# wasting computational effort.
+mean_c = 200e-9
+stddev_c = 50e-9
+mean_a = 100e-9
+stddev_a = 1e-9
+# Use a specific set of particle sizes for the distribution
+# If false, a randomly generated PSD is used. Otherwise the input should
+# be a 2D list of particle radii with 'Npart' rows and 'Nvol' columns.
+specified_psd_c = False
+specified_psd_a = False
+# Initial electrode filling fractions
+# (for disch, anode starts full, cathode starts empty)
+cs0_c = 0.98
+cs0_a = 0.99
+
+[Conductivity]
+# Simulate bulk cathode conductivity (Ohm's Law)?
+# Options: true, false
+simBulkCond_c = true
+simBulkCond_a = false
+# Dimensional conductivity (used if simBulkCond = true), S/m
+sigma_s_c = 0.5
+sigma_s_a = 1e-1
+# Simulate particle connectivity losses (Ohm's Law)?
+# Options: true, false
+simPartCond_c = false
+simPartCond_a = false
+# Conductance between particles, S = 1/Ohm
+G_mean_c = 1e-14
+G_stddev_c = 0
+G_mean_a = 1e-14
+G_stddev_a = 0
+
+[Thermal Parameters]
+# Heat capacity of anode, cathode, electrolyte material in J/(kg*K).
+cp_c = 700
+cp_s = 700
+cp_a = 700
+# Mass density of anode, cathode, and electrolyte materialin kg/m^3.
+rhom_c = 2500
+rhom_s = 1100
+rhom_a = 2500
+# Heat transfer coefficient with the cell separator (W/(m^2*K))
+h_h = 2e4
+# Thermal conductivity in battery cell (only used for dilute electrolyte model), (W/(m*K)).
+# For Stefan-Maxwell concentrated electrolyte, input from props_elyte.py
+k_h_c = 2.1
+k_h_a = 1.7
+k_h_s = 0.16
+# Includes temperature dependence of entropic heat generation if true, does not include
+# if false
+entropy_heat_gen = False
+
+[Geometry]
+# Thicknesses, m
+L_c = 50e-6
+L_a = 50e-6
+L_s = 25e-6
+# Volume loading percents of active material (volume fraction of solid
+# that is active material)
+P_L_c = 0.69
+P_L_a = 0.69
+# Porosities (liquid volume fraction in each region)
+poros_c = 0.5
+poros_a = 0.4
+poros_s = 1.0
+# Bruggeman exponent (tortuosity = porosity^bruggExp)
+BruggExp_c = -0.5
+BruggExp_a = -0.5
+BruggExp_s = -0.5
+
+[Electrolyte]
+# Initial electrolyte conc., mol/m^3
+c0 = 1000
+# Cation/anion charge number (e.g. 2, -1 for CaCl_2)
+zp = 1
+zm = -1
+# Cation/anion dissociation number (e.g. 1, 2 for CaCl_2)
+nup = 1
+num = 1
+# Electrolyte model,
+# Options: dilute, SM
+#   dilute: assume dilute, binary electrolyte model; phi in
+#   electrolyte is an "inner potential" or an "quasi-electrostatic
+#   potential"
+#   SM: use Stefan-Maxwell model for electrolyte transport; phi in
+#   electrolyte is that measured with a Li metal reference electrode
+#   relative to some fixed position in solution.
+# WARNING: Using SM model with BV reaction models for electrode
+# particles assumes individual ion activities are given by their
+# concentrations for the reaction rate exchange current density.
+elyteModelType = SM
+# Stefan-Maxwell property set, see props_elyte.py file
+# Options:
+#   test1: parameter set for testing
+#   LiClO4_PC: electrolyte/solvent used in Fuller, Doyle, Newman 1994
+#              conductivity taken from dualfoil5.2.f
+#   valoen_bernardi: LiPF6 in carbonates as in Bernardi and Go 2011
+#SMset_filename = filename.py   # optional, to load SMset from custom file instead of props_elyte.py
+SMset = valoen_bernardi
+# Reference electrode (defining the electrolyte potential) information:
+# number of electrons transfered in the reaction, 1 for Li/Li+
+n = 1
+# Stoichiometric coefficient of cation, -1 for Li/Li+
+sp = -1
+# Dilute solution properties (used only if elyteModelType = "dilute")
+# Cation/anion diff, m^2/s
+# e.g. for LiPF6 in EC/DMC, Dp = 2.2e-10, Dm = 2.94e-10
+Dp = 2.2e-10
+Dm = 2.94e-10