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

Issue:
Example config/lionsimba
Config files that reproduce the isothermal LIONSIMBA benchmark problem. The results were checked against LIONSIMBA simulations with 10 volumes for anode, cathode, and separator. Note that even LIONSIMBA does not quite match Fig 5c of the paper. I suspect that the results in the paper might have been computed at a different resolution.


Discussion / hints:

Fix patch:
diff --git a/.flake8 b/.flake8
index b5d9ab82..5acc5100 100644
--- a/.flake8
+++ b/.flake8
@@ -1,8 +1,8 @@
 [flake8]
 # W503 line break before binary operator
 # E226 missing whitespace around arithmetic operator
-# E123 closing bracket does not match indentation of opening bracket’s line
+# E123 closing bracket does not match indentation of opening bracket's line
 # E231 missing whitespace after ',', ';', or ':' (should only be ignored for array indexing)
 ignore = W503,E226,E123,E231
-exclude = mpetLogTimeSpace.py,mpetMaterialsMHC.py,mod_cell_M.py
+exclude = mpetLogTimeSpace.py,mpetMaterialsMHC.py,mod_cell_M.py,bin/workdir
 max-line-length = 99
diff --git a/.github/workflows/mpet-code-style.yml b/.github/workflows/mpet-code-style.yml
new file mode 100644
index 00000000..c63a21ed
--- /dev/null
+++ b/.github/workflows/mpet-code-style.yml
@@ -0,0 +1,31 @@
+name: MPET code style check
+
+on: [push, workflow_dispatch]
+
+jobs:
+  codestyle:
+
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Set up python
+        uses: actions/setup-python@v2
+        with:
+          python-version: 3.9
+          architecture: x64
+      
+      - name: Checkout MPET
+        uses: actions/checkout@v2
+        with:
+          fetch-depth: 1
+          path: mpet
+
+      - name: Install additional dependencies
+        run: |
+          pip install flake8
+
+      - name: Check code style
+        run: |
+          cd mpet
+          flake8
+
diff --git a/.gitignore b/.gitignore
index 51f4e752..e775f220 100644
--- a/.gitignore
+++ b/.gitignore
@@ -36,4 +36,5 @@ mpetMaterialsMHC.py
 venv/
 .remarkrc
 .vim
-*/.flake8
+.idea
+*.egg-info
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9571687d..f703dd0d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,16 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
+## [0.1.5] - 2021-01-19
+### Added
+- Improved regression test suite based on pytest.
+- Docker config file for running MPET in a container.
+- Continuous integration implemented for Github and Bitbucket source code repositories.
+
+### Fixed
+- Port variables and redundant time variables are no longer saved as output, reducing file size.
+
+
 ## [0.1.4] - 2020-07-27
 ### Added
 - Adds a separate package for mpet.plot.
diff --git a/README.md b/README.md
index 7c0695ae..457cc389 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,5 @@
+Development: [![Coverage Status](https://coveralls.io/repos/github/TRI-AMDD/mpet-dev/badge.svg?branch=development)](https://coveralls.io/github/TRI-AMDD/mpet-dev?branch=development)
+Master: [![Coverage Status](https://coveralls.io/repos/github/TRI-AMDD/mpet-dev/badge.svg?branch=master)](https://coveralls.io/github/TRI-AMDD/mpet-dev?branch=master)
 # MPET -- Multiphase Porous Electrode Theory
 
 This software is designed to run simulations of batteries with porous electrodes using porous electrode theory, which is a volume-averaged, multiscale approach to capture the coupled behavior of electrolyte and active material within electrodes. As a result, with physical parameter inputs and run protocols (specified current or voltage profiles), it makes predictions about the internal dynamics within a battery (electrolyte concentration and potential, solid phase concentrations, reaction rates, etc.) and also macroscopic, easily measurable electrochemical quantities such as total current and voltage. In this way, it is similar to the [`dualfoil`](http://www.cchem.berkeley.edu/jsngrp/fortran.html) code released by Newman and coworkers from Berkeley. This software has much of the functionality contained in `dualfoil` (it is currently missing, e.g., temperature dependence). However, beyond the standard porous electrode theory simulations, this software can also simulate electrodes in which the active materials phase separate using non-equilibrium thermodynamics within a phase field modeling framework. Such behavior is common in widely used electrode materials, including graphite and LiFePO4.
@@ -13,11 +15,14 @@ Smith, R. B., and Bazant M. Z., Multiphase Porous Electrode Theory, [Journal of
 
 ## Installation
 
-1.  Download the [latest release of MPET](https://bitbucket.org/bazantgroup/mpet/downloads/?tab=tags), or clone a copy of this source code repository.
-2.  Enter the mpet folder, and use the setup.py script to install the mpet Python package.
+1.  Install the prerequisites above.
+2.  Download the [latest release of MPET](https://bitbucket.org/bazantgroup/mpet/downloads/?tab=tags), or clone a copy of this source code repository.
+3.  Enter the mpet folder, and use the setup.py script to install the mpet Python package:
     - We recommend using the pip package manager: `pip install .`
     - The legacy approach also works: `python setup.py install`
 
+MPET is also available on [PyPI](https://pypi.org/project/mpet/), the Python Package Index, and can be installed with `pip install mpet`.
+
 ## Simulation
 
 1.  Copy the overall system parameters file, `configs/params_system.cfg`, to your working directory.
diff --git a/bin/mpetplot.py b/bin/mpetplot.py
index 2c8bc186..26112636 100755
--- a/bin/mpetplot.py
+++ b/bin/mpetplot.py
@@ -11,8 +11,8 @@
 import mpet.plot.outmat2txt as outmat2txt
 import mpet.plot.plot_data as plot_data
 
-#Ordered dictionary of plot types
-plotTypes=OrderedDict([
+# Ordered dictionary of plot types
+plotTypes = OrderedDict([
     ('v','voltage vs filling fraction'),
     ('vt','voltage vs time'),
     ('curr','current vs time'),
@@ -38,15 +38,16 @@
     ('text','convert the output to plain text (csv)')
 ])
 
-#Generate help description for plotTypes
-plotTypesHelp="""The type of plot to be generated:"""
+# Generate help description for plotTypes
+plotTypesHelp = """The type of plot to be generated:"""
 for key,value in plotTypes.items():
     plotTypesHelp = plotTypesHelp + '\n' + key
-    if len(key)<7:
+    if len(key) < 7:
         plotTypesHelp = plotTypesHelp + '\t'
     plotTypesHelp = plotTypesHelp + '\t' + value
 
-parser = argparse.ArgumentParser(description='Process and plot results generated by mpetrun.py.', formatter_class=RawTextHelpFormatter)
+parser = argparse.ArgumentParser(description='Process and plot results generated by mpetrun.py.',
+                                 formatter_class=RawTextHelpFormatter)
 parser.add_argument('directory', help='location of the mpet results')
 parser.add_argument('plotType', metavar='plotType', help=plotTypesHelp, choices=plotTypes.keys())
 parser.add_argument('save',nargs='?',choices=['save'],help='Optionally save the output')
diff --git a/bin/mpetrun.py b/bin/mpetrun.py
index d6c6757f..da5b18a6 100755
--- a/bin/mpetrun.py
+++ b/bin/mpetrun.py
@@ -7,7 +7,7 @@
 from mpet.version import __version__
 import mpet.main as main
 
-desc="""MPET - Multiphase Porous Electrode Theory
+desc = """MPET - Multiphase Porous Electrode Theory
 This software is designed to run simulations of batteries with porous electrodes
 using porous electrode theory, which is a volume-averaged, multiscale approach
 to capture the coupled behavior of electrolyte and active material within
diff --git a/bitbucket-pipelines.yml b/bitbucket-pipelines.yml
new file mode 100644
index 00000000..52d7a4f6
--- /dev/null
+++ b/bitbucket-pipelines.yml
@@ -0,0 +1,84 @@
+image: vikko00/daetools:1.9.0
+options:
+ docker: true
+
+pipelines:
+  pull-requests:
+    '**':
+      - step:
+          caches:
+            - docker
+          script:
+            # update container
+            - apt-get update
+            - apt-get upgrade -y
+            
+            # Install additional dependencies using mpet's setup.py
+            - pip install .[test]
+            
+            # Set up test for modified branch
+            - cd bin
+            - rm -rf workdir
+            - mkdir workdir
+            - cd workdir
+
+            - cp ../run_tests.py .
+            - ln -s ../../mpet .
+            - ln -s ../../tests .
+
+            # run tests for modified branch and get coverage
+            - coverage run --source=../../mpet/ run_tests.py --test_dir ./tests --output_dir ../../bin/workdir/modified > /dev/null
+
+            # upload Coveralls
+            - export COVERALLS_REPO_TOKEN=$COVERALLS_REPO_TOKEN
+            - coveralls || true #Dont fret if if fails
+            
+            # run tests for base branch
+            - mkdir base-mpet
+            - ln -sf base-mpet mpet
+            - git --work-tree=base-mpet/ checkout remotes/origin/$BITBUCKET_PR_DESTINATION_BRANCH -- mpet
+            - python run_tests.py --test_dir ./tests --output_dir ../../bin/workdir/stable > /dev/null
+
+            # Checks test results
+            - cd tests
+            - pytest --baseDir=ref_outputs           --modDir=../bin/workdir/modified compare_tests.py
+            - pytest --baseDir=../bin/workdir/stable --modDir=../bin/workdir/modified compare_tests.py
+
+            #Check Timings
+            - pytest --baseDir=../bin/workdir/stable --modDir=../bin/workdir/modified timings_tests.py || true #Ok to fail for now
+  branches:
+    '{master,development}':
+      - step:
+          caches:
+            - docker
+          script:
+            # update container
+            - apt-get update
+            - apt-get upgrade -y
+
+            # Install additional dependencies using mpet's setup.py
+            - pip install .[test]
+            
+            # Set up test for modified branch
+            - cd bin
+            - rm -rf workdir
+            - mkdir workdir
+            - cd workdir
+
+            - cp ../run_tests.py .
+            - ln -s ../../mpet .
+            - ln -s ../../tests .
+
+            # run tests for modified branch and get coverage
+            - coverage run --source=../../mpet/ run_tests.py --test_dir ./tests --output_dir ../../bin/workdir/modified > /dev/null
+
+            # upload Coveralls
+            - export COVERALLS_REPO_TOKEN=$COVERALLS_REPO_TOKEN
+            - coveralls || true #Dont fret if if fails
+            
+            # Checks test results
+            - cd tests
+            - pytest --baseDir=ref_outputs           --modDir=../bin/workdir/modified compare_tests.py
+            
+            #Check Timings
+            - pytest --baseDir=../bin/workdir/stable --modDir=../bin/workdir/modified timings_tests.py || true #Ok to fail for now
diff --git a/configs/params_LiC6_LIONSIMBA.cfg b/configs/params_LiC6_LIONSIMBA.cfg
new file mode 100644
index 00000000..942f1056
--- /dev/null
+++ b/configs/params_LiC6_LIONSIMBA.cfg
@@ -0,0 +1,29 @@
+#Anode parameters for isothermal benchmark from M. Torchio et al., J. Electrochem. Soc. 163, A1192 (2016).
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = diffn
+discretization = 2.e-7
+shape = sphere
+thickness = 20e-9
+
+[Material]
+muRfunc = LiC6_LIONSIMBA
+logPad = false
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+kappa = 4.0e-7
+B = 0.0
+rho_s = 1.839e28
+D = 3.9e-14
+Dfunc = constant
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = BV_mod01
+k0 = 4.690
+alpha = 0.5
+lambda = 6.26e-20
+Rfilm = 0e-0
diff --git a/configs/params_LiCoO2_LIONSIMBA.cfg b/configs/params_LiCoO2_LIONSIMBA.cfg
new file mode 100644
index 00000000..e0bfc785
--- /dev/null
+++ b/configs/params_LiCoO2_LIONSIMBA.cfg
@@ -0,0 +1,29 @@
+#Cathode parameters for isothermal benchmark from M. Torchio et al., J. Electrochem. Soc. 163, A1192 (2016).
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = diffn
+discretization = 2.e-7
+shape = sphere
+thickness = 20e-9
+
+[Material]
+muRfunc = LiCoO2_LIONSIMBA
+logPad = false
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+kappa = 5.0148e-10
+B = 0.1916e9
+rho_s = 3.1036e28
+D = 1e-14
+Dfunc = constant
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = BV_mod01
+k0 = 3.671
+alpha = 0.5
+lambda = 6.26e-20
+Rfilm = 0e-0
diff --git a/configs/params_system.cfg b/configs/params_system.cfg
index 0d81b608..f5b1dcc4 100644
--- a/configs/params_system.cfg
+++ b/configs/params_system.cfg
@@ -9,9 +9,11 @@ profileType = CC
 # Battery (dis)charge c-rate (only used for CC), number of capacities / hr
 # (positive for discharge, negative for charge)
 Crate = 1
+#Optional nominal 1C current density for the cell, A/m^2
+# 1C_current_density = 12.705
 # Voltage cutoffs, V
 Vmax = 3.6
-Vmin = 2.2
+Vmin = 2.0
 # Battery applied voltage (only used for CV), V
 Vset = 0.12
 # CC/CV segments defining profile for profileType = CCsegments or CVsegments
@@ -62,12 +64,12 @@ Rser = 0.
 # - 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 = 5
+Nvol_c = 10
 Nvol_s = 5
-Nvol_a = 0
+Nvol_a = 10
 # Number of particles per volume for cathode and anode
-Npart_c = 5
-Npart_a = 1
+Npart_c = 2
+Npart_a = 2
 
 [Electrodes]
 # The name of the parameter file describing the cathode particles
@@ -118,7 +120,7 @@ G_stddev_a = 0
 # Thicknesses, m
 L_c = 30e-6
 L_a = 50e-6
-L_s = 100e-6
+L_s = 25e-6
 # Volume loading percents of active material (volume fraction of solid
 # that is active material)
 P_L_c = 0.69
@@ -133,21 +135,8 @@ BruggExp_a = -0.5
 BruggExp_s = -0.5
 
 [Electrolyte]
-#Fraction of mobile Li in equilibrium [nondim] (source Raijmakers 2020)
-delta = 0.999
-#Maximal concentration of host sites [mol m^-3] (source Raijmakers 2020)
-cmax = 500
-#recombination rate constant [m^3 mol^-1 s^-1]
-kr = 0
-
-#mobile electrolyte diffusion constant, m^2/s
-D_li = 1.73e-16
-D_n = 5.69e-16
-
-# Initial electrolyte conc., mol/m^3  = cmax*delta
-# c0 and c0fa below are not needed
-c0 = 640
-c0_fac = 2
+# Initial electrolyte conc., mol/m^3
+c0 = 1000
 # Cation/anion charge number (e.g. 2, -1 for CaCl_2)
 zp = 1
 zm = -1
@@ -162,22 +151,17 @@ num = 1
 #   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.
-#   Solid: Takes the lattice gas model for solid electrolytes.
 # 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 = Solid
-# Solid solution parameter for Solid Electrolyte 
-# a_slyte is the alpha used in the lattice gas model
-a_slyte = 0
+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
-#   Solid_elyte_func: Takes the lattice gas model for solid electrolytes.
-SMset = Solid_elyte_func
+SMset = valoen_bernardi
 # Reference electrode (defining the electrolyte potential) information:
 # number of electrons transfered in the reaction, 1 for Li/Li+
 n = 1
@@ -186,5 +170,5 @@ 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 = 1.18e-10
+Dp = 2.2e-10
 Dm = 2.94e-10
diff --git a/configs/params_system_LIONSIMBA.cfg b/configs/params_system_LIONSIMBA.cfg
new file mode 100644
index 00000000..9479ad3c
--- /dev/null
+++ b/configs/params_system_LIONSIMBA.cfg
@@ -0,0 +1,125 @@
+#Cell parameters for isothermal benchmark from M. Torchio et al., J. Electrochem. Soc. 163, A1192 (2016).
+# See params_system.cfg for parameter explanations.
+
+[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 = 1
+#Optional nominal 1C current density for the cell, A/m^2
+1C_current_density = 30
+# Voltage cutoffs, V
+Vmax = 5
+Vmin = 2.5
+# Battery applied voltage (only used for CV), V
+Vset = 0.12
+# CC/CV segments defining profile for profileType = CCsegments or CVsegments
+segments = [
+            (0.3, 0.4),
+            (-0.5, 0.1),
+            ]
+# Continuation directory.
+# Options: false, absolute directory path
+prevDir = false
+# Final time (only used for CV), [s]
+tend = 1.2e3
+# Number disc. in time
+tsteps = 200
+# Numerical Tolerances
+relTol = 1e-6
+absTol = 1e-6
+# Temp, K
+# WARNING -- temperature dependence not fully implemented. Use 298 K.
+T = 298
+# Random seed. Set to true to give a random seed in the simulation
+randomSeed = false
+# Value of the random seed, must be an integer
+seed = 0
+# Series resistance, [Ohm m^2]
+Rser = 0.
+# Cathode, anode, and separator numer disc. in x direction (volumes in electrodes)
+Nvol_c = 10
+Nvol_s = 10
+Nvol_a = 10
+# Number of particles per volume for cathode and anode
+Npart_c = 1
+Npart_a = 1
+
+[Electrodes]
+cathode = params_LiCoO2_LIONSIMBA.cfg
+anode = params_LiC6_LIONSIMBA.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
+mean_c = 2e-6
+stddev_c = 0
+mean_a = 2e-6
+stddev_a = 0
+# Initial electrode filling fractions
+cs0_c = 0.4995
+cs0_a = 0.8551
+
+[Conductivity]
+# Simulate bulk cathode conductivity (Ohm's Law)?
+simBulkCond_c = true
+simBulkCond_a = true
+# Dimensional conductivity (used if simBulkCond = true), S/m
+sigma_s_c = 412.43
+sigma_s_a = 685.77
+# Simulate particle connectivity losses (Ohm's Law)?
+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
+
+[Geometry]
+# Thicknesses, m
+L_c = 8e-5
+L_a = 8.8e-5
+L_s = 2.5e-5
+# Volume loading percents of active material (volume fraction of solid
+# that is active material)
+P_L_c = 0.9593
+P_L_a = 0.9367
+# Porosities (liquid volume fraction in each region)
+poros_c = 0.385
+poros_a = 0.485
+poros_s = 0.724
+# Bruggeman exponent (tortuosity = porosity^bruggExp)
+BruggExp_c = -3
+BruggExp_a = -3
+BruggExp_s = -3
+
+[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
+elyteModelType = SM
+# Stefan-Maxwell property set, see props_elyte.py file
+SMset = LIONSIMBA_isothermal
+# 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
+Dp = 7.5e-10
+Dm = 7.5e-10
diff --git a/configs/params_system_solid.cfg b/configs/params_system_solid.cfg
new file mode 100644
index 00000000..0b1a82a3
--- /dev/null
+++ b/configs/params_system_solid.cfg
@@ -0,0 +1,196 @@
+# 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 = 1
+# Voltage cutoffs, V
+Vmax = 10
+Vmin = -10
+# Fraction of full (dis)charge to simulate (only used for CC)
+# E.g. to go from cathode filling fraction 0.1 to 0.5 with cathode
+# limiting capacity and Crate>0, set cs0_c = 0.1, capFrac=0.4
+capFrac = 0.98
+# Battery applied voltage (only used for CV), V
+Vset = 0.12
+# 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-6
+# Absolute Tolerance
+absTol = 1e-6
+# Temp, K
+# WARNING -- temperature dependence not fully implemented. Use 298 K.
+T = 298
+# 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
+# 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 = 5
+Nvol_s = 3
+Nvol_a = 0
+# Number of particles per volume for cathode and anode
+Npart_c = 2
+Npart_a = 1
+
+[Electrodes]
+# The name of the parameter file describing the cathode particles
+cathode = params_LFP.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 = 4e0
+# 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 = 100e-9
+stddev_c = 1e-9
+mean_a = 100e-9
+stddev_a = 1e-9
+# 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
+
+[Geometry]
+# Thicknesses, m
+L_c = 50e-6
+L_a = 50e-6
+L_s = 30e-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.4
+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]
+#Fraction of mobile Li in equilibrium [nondim] (source Raijmakers 2020)
+delta = 0.64
+#Maximal concentration of host sites [mol m^-3] (source Raijmakers 2020)
+cmax = 61141
+#recombination rate constant [m^3 mol^-1 s^-1]
+kr = 8e-7
+
+#mobile electrolyte diffusion constant, m^2/s
+D_li = 1.73e-16
+D_n = 5.69e-16
+
+# Initial electrolyte conc., mol/m^3  = cmax*delta
+# c0 and c0fa below are not needed
+c0 = 640
+c0fac = 2
+# 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.
+#   Solid: Takes the lattice gas model for solid electrolytes.
+# 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 = solid
+# Solid solution parameter for Solid Electrolyte 
+# c_dep to increase the effect of solid (not used)
+c_dep = 1
+# a_slyte is the alpha used in the lattice gas model
+a_slyte = 0
+# 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
+#   solid_elyte_func: Takes the lattice gas model for solid electrolytes.
+SMset = solid_elyte
+# 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 = 1.18e-10
+Dm = 2.94e-10
diff --git a/dockerfiles/Dockerfile-daetools b/dockerfiles/Dockerfile-daetools
new file mode 100644
index 00000000..5fe3e553
--- /dev/null
+++ b/dockerfiles/Dockerfile-daetools
@@ -0,0 +1,13 @@
+FROM python:3.7.8-stretch
+
+#Install necessary packages for daetools
+
+RUN apt-get update
+
+RUN apt-get install --yes libgl1-mesa-glx libgfortran3
+
+#Download and install daetools 1.9.0
+RUN wget https://sourceforge.net/projects/daetools/files/daetools/1.9.0/daetools-1.9.0-gnu_linux-x86_64.tar.gz && \ 
+    tar -xf daetools-1.9.0-gnu_linux-x86_64.tar.gz && \
+    pip install ./daetools-1.9.0-gnu_linux && \
+    rm -r ./daetools-1.9.0-gnu_linux
diff --git a/dockerfiles/Dockerfile-mpet-deps b/dockerfiles/Dockerfile-mpet-deps
new file mode 100644
index 00000000..632aa813
--- /dev/null
+++ b/dockerfiles/Dockerfile-mpet-deps
@@ -0,0 +1,6 @@
+FROM vikko00/daetools:1.9.0
+
+#Install mpet dependencies
+
+RUN pip install numpy scipy matplotlib pyQt5 h5py
+RUN pip install coverage coveralls pytest
diff --git a/mpet/__init__.py b/mpet/__init__.py
index 71525559..9ef2efff 100644
--- a/mpet/__init__.py
+++ b/mpet/__init__.py
@@ -1,1 +1,1 @@
-from .version import __version__
\ No newline at end of file
+from .version import __version__  # noqa - don't complain about unused import
diff --git a/mpet/daeVariableTypes.py b/mpet/daeVariableTypes.py
index 3402452f..be994d10 100644
--- a/mpet/daeVariableTypes.py
+++ b/mpet/daeVariableTypes.py
@@ -9,4 +9,4 @@
     upperBound=1e20, initialGuess=1.00, absTolerance=1.e-6)
 elec_pot_t = dae.daeVariableType(
     name="elec_pot_t", units=dae.unit(), lowerBound=-1e20,
-    upperBound=1e20, initialGuess=0, absTolerance=1.e-6)
\ No newline at end of file
+    upperBound=1e20, initialGuess=0, absTolerance=1.e-6)
diff --git a/mpet/data_reporting.py b/mpet/data_reporting.py
index c5b7d1ed..188b4032 100644
--- a/mpet/data_reporting.py
+++ b/mpet/data_reporting.py
@@ -12,6 +12,10 @@ class MyMATDataReporter(daeMatlabMATFileDataReporter):
     """See source code for pyDataReporting.daeMatlabMATFileDataReporter"""
     def WriteDataToFile(self):
         mdict = {}
+
+        # Get the simulation times from the phi_applied variable
+        mdict['phi_applied_times'] = self.Process.dictVariables['mpet.phi_applied'].TimeValues
+
         for var in self.Process.Variables:
             # Remove the model name part of the output key for
             # brevity.
@@ -20,7 +24,7 @@ def WriteDataToFile(self):
             # file to be read by, e.g., MATLAB.
             dkeybase = dkeybase.replace(".", "_")
             mdict[dkeybase] = var.Values
-            mdict[dkeybase + '_times'] = var.TimeValues
+
         sio.savemat(self.ConnectionString,
                     mdict, appendmat=False, format='5',
                     long_field_names=False, do_compression=False,
diff --git a/mpet/electrode/reactions.py b/mpet/electrode/reactions.py
index ff5a26ba..f26c21f7 100644
--- a/mpet/electrode/reactions.py
+++ b/mpet/electrode/reactions.py
@@ -4,50 +4,54 @@
 
 eps = -1e-12
 
-#Reaction rate functions 
-#Should have the same name as the rxnType config field
+
+# Reaction rate functions
+# Should have the same name as the rxnType config field
 def BV(eta, c_sld, c_lyte, k0, T, act_R=None,
-                  act_lyte=None, lmbda=None, alpha=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))
-    ecd = (k0 * act_lyte**(1-alpha)
-               * act_R**(alpha) / gamma_ts)
+    ecd = (k0 * act_lyte**(1-alpha) * act_R**(alpha) / gamma_ts)
     Rate = ecd * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T))
     return Rate
 
+
 def BV_gMod01(eta, c_sld, c_lyte, k0, T, act_R=None,
-                  act_lyte=None, lmbda=None, alpha=None):
+              act_lyte=None, lmbda=None, alpha=None):
     if act_R is None:
         act_R = c_sld/(1-c_sld)
     gamma_ts = (1./(c_sld*(1-c_sld)))
-    ecd = (k0 * act_lyte**(1-alpha)
-               * act_R**(alpha) / gamma_ts)
+    ecd = (k0 * act_lyte**(1-alpha) * act_R**(alpha) / gamma_ts)
     Rate = ecd * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T))
     return Rate
 
+
 def BV_raw(eta, c_sld, c_lyte, k0, T, act_R=None,
-                  act_lyte=None, lmbda=None, alpha=None):
-    ecd=k0
+           act_lyte=None, lmbda=None, alpha=None):
+    ecd = k0
     Rate = ecd * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T))
     return Rate
 
+
 def BV_mod01(eta, c_sld, c_lyte, k0, T, act_R=None,
-                  act_lyte=None, lmbda=None, alpha=None):
+             act_lyte=None, lmbda=None, alpha=None):
     ecd = (k0 * c_lyte**(1-alpha)
            * (1.0 - c_sld)**(1 - alpha) * c_sld**alpha)
     Rate = ecd * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T))
     return Rate
 
+
 def BV_mod02(eta, c_sld, c_lyte, k0, T, act_R=None,
-                  act_lyte=None, lmbda=None, alpha=None):
+             act_lyte=None, lmbda=None, alpha=None):
     ecd = (k0 * c_lyte**(1-alpha)
            * (0.5 - c_sld)**(1 - alpha) * c_sld**alpha)
     Rate = ecd * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T))
     return Rate
 
+
 def Marcus(eta, c_sld, c_lyte, k0, T, act_R=None,
-                  act_lyte=None, lmbda=None, alpha=None):
+           act_lyte=None, lmbda=None, alpha=None):
     if isinstance(c_sld, np.ndarray):
         c_sld = np.array([
             dae.Max(eps, c_sld[i]) for i in range(len(c_sld))])
@@ -57,11 +61,11 @@ def Marcus(eta, c_sld, c_lyte, k0, T, act_R=None,
     # We'll assume c_e = 1 (at the standard state for electrons)
 #        ecd = ( k0 * np.exp(-lmbda/(4.*T)) *
 #        ecd = ( k0 *
-    ecd = (k0 * (1-c_sld) *
-           c_lyte**((3-2*alpha)/4.) *
-           c_sld**((1+2*alpha)/4.))
-    Rate = (ecd * np.exp(-eta**2/(4.*T*lmbda)) *
-            (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T)))
+    ecd = (k0 * (1-c_sld)
+           * c_lyte**((3-2*alpha)/4.)
+           * c_sld**((1+2*alpha)/4.))
+    Rate = (ecd * np.exp(-eta**2/(4.*T*lmbda))
+            * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T)))
     return Rate
 
 
@@ -78,10 +82,10 @@ def MHC_kfunc(eta, lmbda):
 
 
 def MHC(eta, c_sld, c_lyte, k0, T, act_R=None,
-                  act_lyte=None, lmbda=None, alpha=None):
+        act_lyte=None, lmbda=None, alpha=None):
     # See Zeng, Smith, Bai, Bazant 2014
     # Convert to "MHC overpotential"
-    k0=k0/MHC_kfunc(0., lmbda)
+    k0 = k0/MHC_kfunc(0., lmbda)
     eta_f = eta + T*np.log(c_lyte/c_sld)
     gamma_ts = 1./(1. - c_sld)
     alpha = 0.5
@@ -96,4 +100,4 @@ def MHC(eta, c_sld, c_lyte, k0, T, act_R=None,
         krd = k0*MHC_kfunc(-eta_f, lmbda)
         kox = k0*MHC_kfunc(eta_f, lmbda)
         Rate = ecd_extras*(krd*c_lyte - kox*c_sld)
-    return Rate
\ No newline at end of file
+    return Rate
diff --git a/mpet/extern_funcs.py b/mpet/extern_funcs.py
index b1df52d8..c220279f 100644
--- a/mpet/extern_funcs.py
+++ b/mpet/extern_funcs.py
@@ -52,8 +52,8 @@ def Calculate(self, values):
             return adouble(0)
         # Store the previous time value to prevent excessive
         # interpolation.
-        if (len(self.previous_output) > 0 and
-                self.previous_output[0] == time.Value):
+        if (len(self.previous_output) > 0
+                and self.previous_output[0] == time.Value):
             return adouble(float(self.previous_output[1][self.position]))
         noise_vec = self.interp(time.Value)
         # it is a list now not a tuple
diff --git a/mpet/geometry.py b/mpet/geometry.py
index b1233577..2b3fc98d 100644
--- a/mpet/geometry.py
+++ b/mpet/geometry.py
@@ -80,11 +80,12 @@ def calc_curv(c, dr, r_vec, Rs, beta_s, particleShape):
 
 def get_elyte_disc(Nvol, L, poros, BruggExp):
     out = {}
-    # Discretization
+    # Width of each cell
     out["dxvec"] = utils.get_dxvec(L, Nvol)
+
+    # Distance between cell centers
     dxtmp = np.hstack((out["dxvec"][0], out["dxvec"], out["dxvec"][-1]))
     out["dxd1"] = utils.mean_linear(dxtmp)
-    out["dxd2"] = out["dxvec"]
 
     # The porosity vector
     out["porosvec"] = utils.get_asc_vec(poros, Nvol)
@@ -94,7 +95,6 @@ def get_elyte_disc(Nvol, L, poros, BruggExp):
     Brugg_pad = utils.pad_vec(utils.get_asc_vec(BruggExp, Nvol))
 
     # Vector of posority/tortuosity (assuming Bruggeman)
-    porostortvec_pad = porosvec_pad/porosvec_pad**(Brugg_pad)
-    out["eps_o_tau_edges"] = utils.mean_harmonic(porostortvec_pad)
+    out["eps_o_tau"] = porosvec_pad/porosvec_pad**(Brugg_pad)
 
     return out
diff --git a/mpet/io_utils.py b/mpet/io_utils.py
index d3e6716a..0721e607 100644
--- a/mpet/io_utils.py
+++ b/mpet/io_utils.py
@@ -46,7 +46,7 @@ def get_configs(paramfile="params.cfg"):
     return P_s, P_e
 
 
-def get_dicts_from_configs(P_s, P_e):
+def get_dicts_from_configs(P_s, P_e, paramfile):
     # Dictionary of dimensional system and electrode parameters
     dD_s = {}
     dD_e = {}
@@ -61,13 +61,19 @@ def get_dicts_from_configs(P_s, P_e):
     # Simulation parameters
     ndD_s["profileType"] = P_s.get('Sim Params', 'profileType')
     dD_s["Crate"] = P_s.getfloat('Sim Params', 'Crate')
+    dD_s["CrateCurr"] = P_s.getfloat('Sim Params', '1C_current_density',fallback=None)
     segs = dD_s["segments"] = ast.literal_eval(
         P_s.get('Sim Params', 'segments'))
     ndD_s["tramp"] = dD_s["tramp"] = P_s.getfloat('Sim Params', 'tramp', fallback=0)
     numsegs = dD_s["numsegments"] = len(segs)
     dD_s["Vmax"] = P_s.getfloat('Sim Params', 'Vmax')
     dD_s["Vmin"] = P_s.getfloat('Sim Params', 'Vmin')
-    ndD_s["c0_fac"] = P_s.getfloat('Electrolyte', 'c0_fac')
+    try:
+        ndD_s["c0_fac"] = P_s.getfloat('Electrolyte', 'c0_fac')
+    except configparser.NoOptionError:
+        # not available when not using solid elyte, but then
+        # also not required so ignore error
+        pass
     # Should have deprecation warnings to encourage users to
     # update their params files to mirror options in
     # configDefaults.
@@ -83,6 +89,13 @@ def get_dicts_from_configs(P_s, P_e):
     ndD_s["capFrac"] = P_s.getfloat('Sim Params', 'capFrac',fallback=1)
     dD_s["tend"] = P_s.getfloat('Sim Params', 'tend')
     ndD_s["prevDir"] = P_s.get('Sim Params', 'prevDir')
+
+    # If prevDir is a relative path, convert to absolute
+    if ndD_s["prevDir"] != "false":
+        if not os.path.isabs(ndD_s["prevDir"]):
+            dir = os.path.dirname(paramfile)
+            ndD_s["prevDir"] = os.path.normpath(os.path.join(dir,ndD_s["prevDir"]))
+
     ndD_s["tsteps"] = P_s.getint('Sim Params', 'tsteps')
     Tabs = dD_s["Tabs"] = P_s.getfloat('Sim Params', 'T')
     dD_s["Rser"] = P_s.getfloat('Sim Params', 'Rser')
@@ -133,13 +146,6 @@ def get_dicts_from_configs(P_s, P_e):
                          "s": P_s.getfloat('Geometry', 'BruggExp_s')}
 
     # Electrolyte
-    cmax = dD_s["cmax"] = P_s.getfloat('Electrolyte', 'cmax')
-    delta = ndD_s["delta"] = P_s.getfloat('Electrolyte', 'delta')
-    kr = dD_s["kr"] = P_s.getfloat('Electrolyte', 'kr')
-    kd = dD_s["kd"] = kr*cmax*delta**2/(1-delta)
-    #c0 = dD_s["c0"] = P_s.getfloat('Electrolyte', 'c0')
-    c0 = dD_s["c0"] = delta*cmax
-    a_slyte = ndD_s["a_slyte"] = P_s.getfloat('Electrolyte', 'a_slyte') #new for SE
     zp = ndD_s["zp"] = P_s.getfloat('Electrolyte', 'zp')
     zm = ndD_s["zm"] = P_s.getfloat('Electrolyte', 'zm')
     if zm > 0:
@@ -147,6 +153,17 @@ def get_dicts_from_configs(P_s, P_e):
     ndD_s["nup"] = P_s.getfloat('Electrolyte', 'nup')
     ndD_s["num"] = P_s.getfloat('Electrolyte', 'num')
     ndD_s["elyteModelType"] = P_s.get('Electrolyte', 'elyteModelType')
+    if ndD_s["elyteModelType"] == 'solid':
+        cmax = dD_s["cmax"] = P_s.getfloat('Electrolyte', 'cmax')
+        delta = ndD_s["delta"] = P_s.getfloat('Electrolyte', 'delta')
+        # a_slyte variable is never used
+        # a_slyte = ndD_s["a_slyte"] = P_s.getfloat('Electrolyte', 'a_slyte')
+        ndD_s["a_slyte"] = P_s.getfloat('Electrolyte', 'a_slyte')
+        kr = dD_s["kr"] = P_s.getfloat('Electrolyte', 'kr')
+        kd = dD_s["kd"] = kr * cmax * delta ** 2 / (1 - delta)
+        c0 = dD_s["c0"] = delta * cmax
+    else:
+        c0 = dD_s["c0"] = P_s.getfloat('Electrolyte', 'c0')
     SMset = ndD_s["SMset"] = P_s.get('Electrolyte', 'SMset')
     D_ref = dD_s["D_ref"] = dD_s["Dref"] = getattr(props_elyte,SMset)()[-1]
     ndD_s["n_refTrode"] = P_s.getfloat('Electrolyte', 'n')
@@ -212,7 +229,11 @@ def get_dicts_from_configs(P_s, P_e):
 
     # Various calculated and defined parameters
     L_ref = dD_s["L_ref"] = dD_s["Lref"] = dD_s["L"]["c"]
-    c_ref = dD_s["c_ref"] = dD_s["cref"] = ndD_s["c0_fac"]*c0  # mol/m^3 = 1 M
+    try:
+        c_ref = dD_s["c_ref"] = dD_s["cref"] = ndD_s["c0_fac"] * c0  # mol/m^3 = 1 M
+    except KeyError:
+        # use default if c0_fac not defined in config
+        c_ref = dD_s["c_ref"] = dD_s["cref"] = 1000.  # mol/m^3 = 1 M
     # Ambipolar diffusivity
     Damb = dD_s["Damb"] = ((zp-zm)*Dp*Dm)/(zp*Dp-zm*Dm)
     # Cation transference number
@@ -241,24 +262,29 @@ def get_dicts_from_configs(P_s, P_e):
         # flat plate anode with assumed infinite supply of metal
         ndD_s['z'] = 0.
     limtrode = ("c" if ndD_s["z"] < 1 else "a")
-    CrateCurr = dD_s["CrateCurr"] = dD_e[limtrode]["cap"] / 3600.  # A/m^2
-    dD_s["currset"] = CrateCurr * dD_s["Crate"]  # A/m^2
-    Rser_ref = dD_s["Rser_ref"] = (k*T_ref/e) / (curr_ref*CrateCurr)
+
+    # If 1C_current_density is not defined, use the theoretical capacity
+    theoretical_1C_current = dD_e[limtrode]["cap"] / 3600.  # A/m^2
+    if not dD_s["CrateCurr"]:
+        dD_s["CrateCurr"] = theoretical_1C_current
+
+    dD_s["currset"] = dD_s["CrateCurr"] * dD_s["Crate"]  # A/m^2
+    Rser_ref = dD_s["Rser_ref"] = (k*T_ref/e) / (curr_ref*dD_s["CrateCurr"])
 
     # Some nondimensional parameters
     ndD_s["T"] = Tabs / T_ref
     ndD_s["Rser"] = dD_s["Rser"] / Rser_ref
     ndD_s["Dp"] = Dp / D_ref
     ndD_s["Dm"] = Dm / D_ref
-    #elyte initially in equilibrium:
-    ndD_s["c0"] = delta*cmax/c_ref
-    ndD_s["cmax"] = cmax/c_ref
-    ndD_s["kr"] = kr *t_ref*c_ref#/(L_ref**3) *t_ref/N_A
-    ndD_s["kd"] = kd *t_ref
+    ndD_s["c0"] = c0 / c_ref
     ndD_s["phi_cathode"] = 0.
-    ndD_s["currset"] = dD_s["Crate"] / curr_ref
-    ndD_s["k0_foil"] = dD_s["k0_foil"] * (1./(curr_ref*CrateCurr))
+    ndD_s["currset"] = dD_s["currset"] / theoretical_1C_current / curr_ref
+    ndD_s["k0_foil"] = dD_s["k0_foil"] * (1./(curr_ref*dD_s["CrateCurr"]))
     ndD_s["Rfilm_foil"] = dD_s["Rfilm_foil"] / Rser_ref
+    if ndD_s["elyteModelType"] == 'solid':
+        ndD_s["cmax"] = cmax / c_ref
+        ndD_s["kr"] = kr * t_ref * c_ref
+        ndD_s["kd"] = kd * t_ref
 
     # parameters which depend on the electrode
     dD_s["psd_raw"] = {}
@@ -375,15 +401,17 @@ def get_dicts_from_configs(P_s, P_e):
     ndD_s["phimin"] = -((e/(k*T_ref))*dD_s["Vmax"] + ndDVref)
     ndD_s["phimax"] = -((e/(k*T_ref))*dD_s["Vmin"] + ndDVref)
 
-    #Nondimensionalize current and voltage segments
+    # Nondimensionalize current and voltage segments
     ndD_s["segments"] = []
     if ndD_s["profileType"] == "CCsegments":
         for i in range(len(dD_s["segments"])):
-            ndD_s["segments"].append((dD_s["segments"][i][0]/curr_ref, dD_s["segments"][i][1]*60/t_ref))
+            ndD_s["segments"].append((dD_s["segments"][i][0]/curr_ref,
+                                      dD_s["segments"][i][1]*60/t_ref))
     elif ndD_s["profileType"] == "CVsegments":
         for i in range(len(dD_s["segments"])):
-            ndD_s["segments"].append((-((e/(k*T_ref))*dD_s["segments"][i][0]+ndDVref), dD_s["segments"][i][1]*60/t_ref))
-        
+            ndD_s["segments"].append((-((e/(k*T_ref))*dD_s["segments"][i][0]+ndDVref),
+                                      dD_s["segments"][i][1]*60/t_ref))
+
     # Current or voltage segments profiles
     dD_s["segments_tvec"] = np.zeros(2*numsegs + 1)
     dD_s["segments_setvec"] = np.zeros(2*numsegs + 1)
@@ -550,8 +578,8 @@ def test_electrode_input(dD, ndD, dD_s, ndD_s):
     if (solidType in ["CHR", "diffn"] and solidShape not in
             ["sphere", "cylinder"]):
         raise NotImplementedError("CHR and diffn req. sphere or cylinder")
-    if ((solidType not in ndD_s["1varTypes"]) and
-            (solidType not in ndD_s["2varTypes"])):
+    if ((solidType not in ndD_s["1varTypes"])
+            and (solidType not in ndD_s["2varTypes"])):
         raise NotImplementedError("Input solidType not defined")
     if solidShape not in ["C3", "sphere", "cylinder"]:
         raise NotImplementedError("Input solidShape not defined")
diff --git a/mpet/main.py b/mpet/main.py
index e68794ba..634791f7 100644
--- a/mpet/main.py
+++ b/mpet/main.py
@@ -33,6 +33,22 @@ def run_simulation(ndD_s, ndD_e, tScale, outdir):
     # Enable reporting of all variables
     simulation.m.SetReportingOn(True)
 
+    # Turn off reporting of some variables
+    simulation.m.endCondition.ReportingOn = False
+
+    # Turn off reporting of particle ports
+    for trode in simulation.m.trodes:
+        for particle in simulation.m.particles[trode]:
+            pModel = particle[0]
+            for port in pModel.Ports:
+                for var in port.Variables:
+                    var.ReportingOn = False
+
+    # Turn off reporting of cell ports
+    for port in simulation.m.Ports:
+        for var in port.Variables:
+            var.ReportingOn = False
+
     # Set relative tolerances
     daesolver.RelativeTolerance = ndD_s["relTol"]
 
@@ -75,7 +91,7 @@ def main(paramfile, keepArchive=True):
     # Get the parameters dictionary (and the config instance) from the
     # parameter file
     P_s, P_e = IO.get_configs(paramfile)
-    dD_s, ndD_s, dD_e, ndD_e = IO.get_dicts_from_configs(P_s, P_e)
+    dD_s, ndD_s, dD_e, ndD_e = IO.get_dicts_from_configs(P_s, P_e, paramfile)
 
     # Directories we'll store output in.
     outdir_name = time.strftime("%Y%m%d_%H%M%S", time.localtime())
@@ -117,13 +133,13 @@ def main(paramfile, keepArchive=True):
     except subp.CalledProcessError:
         pass
 
-    fo=open(os.path.join(outdir, 'run_info.txt'), 'w')
-    
-    #Print mpet version
+    fo = open(os.path.join(outdir, 'run_info.txt'), 'w')
+
+    # Print mpet version
     print("mpet version:", file=fo)
     print(mpet.__version__+"\n", file=fo)
 
-    #Print git commit info if it exists
+    # Print git commit info if it exists
     if commit_hash != "":
         # Store commit info to file, as well as how to patch if
         # there's a diff
@@ -147,22 +163,23 @@ def main(paramfile, keepArchive=True):
         pyFiles = glob.glob(os.path.join(localDir, "*.py"))
         for pyFile in pyFiles:
             shutil.copy(pyFile, snapshotDir)
-    
+
     fo.close()
-    
+
     # External functions are not supported by the Compute Stack approach.
     # Activate the Evaluation Tree approach if noise, logPad, CCsegments,
     # or CVsegments are used
     cfg = dae.daeGetConfig()
-    noise=ndD_e['c']['noise']
-    logPad=ndD_e['c']['logPad']
+    noise = ndD_e['c']['noise']
+    logPad = ndD_e['c']['logPad']
     segments = ndD_s["profileType"] in ["CCsegments","CVsegments"]
-    if (noise or logPad or (segments and ndD_s["tramp"]>0)) and cfg.has_key('daetools.core.equations.evaluationMode'):
+    if (noise or logPad or (segments and ndD_s["tramp"] > 0)) \
+            and 'daetools.core.equations.evaluationMode' in cfg:
         cfg.SetString('daetools.core.equations.evaluationMode', 'evaluationTree_OpenMP')
     with open(os.path.join(outdir, "daetools_config_options.txt"), 'w') as fo:
         print(cfg, file=fo)
 
-    #Disable printStats
+    # Disable printStats
     cfg.SetString('daetools.activity.printStats','false')
 
     # Carry out the simulation
diff --git a/mpet/mod_cell.py b/mpet/mod_cell.py
index 299a360f..bd5d4f1d 100644
--- a/mpet/mod_cell.py
+++ b/mpet/mod_cell.py
@@ -8,7 +8,7 @@
  - potential drop between simulated particles
 """
 import daetools.pyDAE as dae
-from pyUnits import m, kg, s, K, Pa, mol, J, W
+from pyUnits import s
 
 import numpy as np
 
@@ -18,13 +18,14 @@
 import mpet.ports as ports
 import mpet.props_elyte as props_elyte
 import mpet.utils as utils
-from mpet.daeVariableTypes import *
+from mpet.daeVariableTypes import mole_frac_t, elec_pot_t, conc_t
 
-#Dictionary of end conditions
-endConditions={
+# Dictionary of end conditions
+endConditions = {
     1:"Vmax reached",
     2:"Vmin reached"}
 
+
 class ModCell(dae.daeModel):
     def __init__(self, Name, Parent=None, Description="", ndD_s=None,
                  ndD_e=None):
@@ -51,8 +52,8 @@ def __init__(self, Name, Parent=None, Description="", ndD_s=None,
                 "Simulated volumes in electrode {trode}".format(trode=trode))
             self.DmnPart[trode] = dae.daeDomain(
                 "Npart_{trode}".format(trode=trode), self, dae.unit(),
-                "Particles sampled in each control " +
-                "volume in electrode {trode}".format(trode=trode))
+                "Particles sampled in each control "
+                + "volume in electrode {trode}".format(trode=trode))
 
         # Variables
         self.c_lyte = {}
@@ -212,8 +213,8 @@ def DeclareEquations(self):
                     eq = self.CreateEquation(
                         "portout_pm_trode{trode}v{vInd}p{pInd}".format(
                             vInd=vInd, pInd=pInd, trode=trode))
-                    eq.Residual = (self.phi_part[trode](vInd, pInd) -
-                                   self.portsOutBulk[trode][vInd,pInd].phi_m())
+                    eq.Residual = (self.phi_part[trode](vInd, pInd)
+                                   - self.portsOutBulk[trode][vInd,pInd].phi_m())
 
             # Simulate the potential drop along the bulk electrode
             # solid phase
@@ -297,6 +298,9 @@ def DeclareEquations(self):
             # Ghost points on the left and no-gradients on the right
             ctmp = np.hstack((self.c_lyteGP_L(), cvec, cvec[-1]))
             phitmp = np.hstack((self.phi_lyteGP_L(), phivec, phivec[-1]))
+
+            Nm_edges, i_edges = get_lyte_internal_fluxes(ctmp, phitmp, disc, ndD)
+
             # If we don't have a porous anode:
             # 1) the total current flowing into the electrolyte is set
             # 2) assume we have a Li foil with BV kinetics and the specified rate constant
@@ -304,13 +308,11 @@ def DeclareEquations(self):
             eqP = self.CreateEquation("GhostPointP_L")
             if Nvol["a"] == 0:
                 # Concentration BC from mass flux
-                Nm_foil = get_lyte_internal_fluxes(
-                    ctmp[0:2], phitmp[0:2], disc["dxd1"][0], disc["eps_o_tau_edges"][0], ndD)[0]
-                eqC.Residual = Nm_foil[0]
+                eqC.Residual = Nm_edges[0]
                 # Phi BC from BV at the foil
                 # We assume BV kinetics with alpha = 0.5,
                 # exchange current density, ecd = k0_foil * c_lyte**(0.5)
-                cWall = utils.mean_harmonic(ctmp[0], ctmp[1])
+                cWall = .5*(ctmp[0] + ctmp[1])
                 ecd = ndD["k0_foil"]*cWall**0.5
                 # -current = ecd*(exp(-eta/2) - exp(eta/2))
                 # note negative current because positive current is
@@ -329,21 +331,21 @@ def DeclareEquations(self):
                 if ndD["elyteModelType"] == "dilute":
                     phiWall -= ndD["T"]*np.log(cWall)
                 # phiWall = 0.5 * (phitmp[0] + phitmp[1])
-                eqP.Residual = phiWall - utils.mean_linear(phitmp[0], phitmp[1])
+                eqP.Residual = phiWall - .5*(phitmp[0] + phitmp[1])
             # We have a porous anode -- no flux of charge or anions through current collector
             else:
                 eqC.Residual = ctmp[0] - ctmp[1]
                 eqP.Residual = phitmp[0] - phitmp[1]
 
-            Nm_edges, i_edges = get_lyte_internal_fluxes(
-                ctmp, phitmp, disc["dxd1"], disc["eps_o_tau_edges"], ndD)
-            dvgNm = np.diff(Nm_edges)/disc["dxd2"]
-            dvgi = np.diff(i_edges)/disc["dxd2"]
+            dvgNm = np.diff(Nm_edges)/disc["dxvec"]
+            dvgi = np.diff(i_edges)/disc["dxvec"]
             for vInd in range(Nlyte):
-                # Mass Conservation (done with the anion, although "c" is neutral salt conc) solid electrolyte added
+                # Mass Conservation (done with the anion, although "c" is neutral salt conc)
                 eq = self.CreateEquation("lyte_mass_cons_vol{vInd}".format(vInd=vInd))
-                eq.Residual = (disc["porosvec"][vInd]*dcdtvec[vInd] + (1./ndD["num"])*dvgNm[vInd]
-                -ndD["kd"]*(ndD["cmax"]-cvec[vInd])+ndD["kr"]*cvec[vInd]*cvec[vInd])
+                eq.Residual = disc["porosvec"][vInd]*dcdtvec[vInd] + (1./ndD["num"])*dvgNm[vInd]
+                if ndD["elyteModelType"] == "solid":
+                    eq.Residual += -ndD["kd"] * (ndD["cmax"] - cvec[vInd]) \
+                        + ndD["kr"] * cvec[vInd] ** 2
                 # Charge Conservation
                 eq = self.CreateEquation("lyte_charge_cons_vol{vInd}".format(vInd=vInd))
                 eq.Residual = -dvgi[vInd] + ndD["zp"]*Rvvec[vInd]
@@ -371,7 +373,7 @@ def DeclareEquations(self):
         if self.profileType == "CC":
             # Total Current Constraint Equation
             eq = self.CreateEquation("Total_Current_Constraint")
-            if ndD["tramp"]>0:
+            if ndD["tramp"] > 0:
                 eq.Residual = self.current() - (
                     ndD["currPrev"] + (ndD["currset"] - ndD["currPrev"])
                     * (1 - np.exp(-dae.Time()/(ndD["tend"]*ndD["tramp"]))))
@@ -380,7 +382,7 @@ def DeclareEquations(self):
         elif self.profileType == "CV":
             # Keep applied potential constant
             eq = self.CreateEquation("applied_potential")
-            if ndD["tramp"]>0:
+            if ndD["tramp"] > 0:
                 eq.Residual = self.phi_applied() - (
                     ndD["phiPrev"] + (ndD["Vset"] - ndD["phiPrev"])
                     * (1 - np.exp(-dae.Time()/(ndD["tend"]*ndD["tramp"])))
@@ -388,60 +390,60 @@ def DeclareEquations(self):
             else:
                 eq.Residual = self.phi_applied() - ndD["Vset"]
         elif self.profileType == "CCsegments":
-            if ndD["tramp"]>0:
+            if ndD["tramp"] > 0:
                 ndD["segments_setvec"][0] = ndD["currPrev"]
                 self.segSet = extern_funcs.InterpTimeScalar(
                     "segSet", self, dae.unit(), dae.Time(),
                     ndD["segments_tvec"], ndD["segments_setvec"])
                 eq = self.CreateEquation("Total_Current_Constraint")
                 eq.Residual = self.current() - self.segSet()
-            
-            #CCsegments implemented as discontinuous equations
+
+            # CCsegments implemented as discontinuous equations
             else:
-                #First segment
-                time=ndD["segments"][0][1]
-                self.IF(dae.Time()<dae.Constant(time*s), 1.e-3)
+                # First segment
+                time = ndD["segments"][0][1]
+                self.IF(dae.Time() < dae.Constant(time*s), 1.e-3)
                 eq = self.CreateEquation("Total_Current_Constraint")
                 eq.Residual = self.current() - ndD["segments"][0][0]
 
-                #Middle segments
+                # Middle segments
                 for i in range(1,len(ndD["segments"])-1):
-                    time=time+ndD["segments"][i][1]
-                    self.ELSE_IF(dae.Time()<dae.Constant(time*s), 1.e-3)
+                    time = time+ndD["segments"][i][1]
+                    self.ELSE_IF(dae.Time() < dae.Constant(time*s), 1.e-3)
                     eq = self.CreateEquation("Total_Current_Constraint")
                     eq.Residual = self.current() - ndD["segments"][i][0]
 
-                #Last segment
+                # Last segment
                 self.ELSE()
                 eq = self.CreateEquation("Total_Current_Constraint")
                 eq.Residual = self.current() - ndD["segments"][-1][0]
                 self.END_IF()
 
         elif self.profileType == "CVsegments":
-            if ndD["tramp"]>0:
+            if ndD["tramp"] > 0:
                 ndD["segments_setvec"][0] = ndD["phiPrev"]
                 self.segSet = extern_funcs.InterpTimeScalar(
                     "segSet", self, dae.unit(), dae.Time(),
                     ndD["segments_tvec"], ndD["segments_setvec"])
                 eq = self.CreateEquation("applied_potential")
                 eq.Residual = self.phi_applied() - self.segSet()
-            
-            #CVsegments implemented as discontinuous equations
+
+            # CVsegments implemented as discontinuous equations
             else:
-                #First segment
-                time=ndD["segments"][0][1]
-                self.IF(dae.Time()<dae.Constant(time*s), 1.e-3)
+                # First segment
+                time = ndD["segments"][0][1]
+                self.IF(dae.Time() < dae.Constant(time*s), 1.e-3)
                 eq = self.CreateEquation("applied_potential")
                 eq.Residual = self.phi_applied() - ndD["segments"][0][0]
 
-                #Middle segments
+                # Middle segments
                 for i in range(1,len(ndD["segments"])-1):
-                    time=time+ndD["segments"][i][1]
-                    self.ELSE_IF(dae.Time()<dae.Constant(time*s), 1.e-3)
+                    time = time+ndD["segments"][i][1]
+                    self.ELSE_IF(dae.Time() < dae.Constant(time*s), 1.e-3)
                     eq = self.CreateEquation("applied_potential")
                     eq.Residual = self.phi_applied() - ndD["segments"][i][0]
 
-                #Last segment
+                # Last segment
                 self.ELSE()
                 eq = self.CreateEquation("applied_potential")
                 eq.Residual = self.phi_applied() - ndD["segments"][-1][0]
@@ -450,27 +452,35 @@ def DeclareEquations(self):
         for eq in self.Equations:
             eq.CheckUnitsConsistency = False
 
-        #Ending conditions for the simulation
+        # Ending conditions for the simulation
         if self.profileType in ["CC", "CCsegments"]:
-            #Vmax reached
+            # Vmax reached
             self.ON_CONDITION((self.phi_applied() <= ndD["phimin"])
-                            & (self.endCondition() < 1),
+                              & (self.endCondition() < 1),
                               setVariableValues=[(self.endCondition, 1)])
-            
-            #Vmin reached
+
+            # Vmin reached
             self.ON_CONDITION((self.phi_applied() >= ndD["phimax"])
-                            & (self.endCondition() < 1),
+                              & (self.endCondition() < 1),
                               setVariableValues=[(self.endCondition, 2)])
 
 
-def get_lyte_internal_fluxes(c_lyte, phi_lyte, dxd1, eps_o_tau, ndD):
+def get_lyte_internal_fluxes(c_lyte, phi_lyte, disc, ndD):
     zp, zm, nup, num = ndD["zp"], ndD["zm"], ndD["nup"], ndD["num"]
     nu = nup + num
     T = ndD["T"]
-    c_edges_int = utils.mean_harmonic(c_lyte)
+    dxd1 = disc["dxd1"]
+    eps_o_tau = disc["eps_o_tau"]
+
+    # Get concentration at cell edges using weighted mean
+    wt = utils.pad_vec(disc["dxvec"])
+    c_edges_int = utils.weighted_linear_mean(c_lyte, wt)
+
     if ndD["elyteModelType"] == "dilute":
-        Dp = eps_o_tau * ndD["Dp"]
-        Dm = eps_o_tau * ndD["Dm"]
+        # Get porosity at cell edges using weighted harmonic mean
+        eps_o_tau_edges = utils.weighted_linear_mean(eps_o_tau, wt)
+        Dp = eps_o_tau_edges * ndD["Dp"]
+        Dm = eps_o_tau_edges * ndD["Dm"]
 #        Np_edges_int = nup*(-Dp*np.diff(c_lyte)/dxd1
 #                            - Dp*zp*c_edges_int*np.diff(phi_lyte)/dxd1)
         Nm_edges_int = num*(-Dm*np.diff(c_lyte)/dxd1
@@ -478,45 +488,38 @@ def get_lyte_internal_fluxes(c_lyte, phi_lyte, dxd1, eps_o_tau, ndD):
         i_edges_int = (-((nup*zp*Dp + num*zm*Dm)*np.diff(c_lyte)/dxd1)
                        - (nup*zp**2*Dp + num*zm**2*Dm)/T*c_edges_int*np.diff(phi_lyte)/dxd1)
 #        i_edges_int = zp*Np_edges_int + zm*Nm_edges_int
-
-    elif ndD["elyteModelType"] == "Solid":
-        D_fs, sigma_fs, thermFac, tp0 = getattr(props_elyte,ndD["SMset"])()[:-1]
-        # modify the free solution transport properties for porous media
-
-        def D(c):
-            return eps_o_tau*D_fs(c)
-
-        def sigma(c):
-            return eps_o_tau*sigma_fs(c)
-        sp, n = ndD["sp"], ndD["n_refTrode"]
-        
-        #D_fs is specified in solid_elyte_func in props_elyte.py
-        Dp = ndD["Dp"]
-        Dm = ndD["Dm"]
-        a_slyte = ndD["a_slyte"]
-        k = 1.381e-23
- 
-        i_edges_int = (-((nup*zp*D(c_edges_int)*(1/(1-c_edges_int)-a_slyte/(k*T)*2*c_edges_int) + num*zm*Dm)*np.diff(c_lyte)/dxd1)
-                       - (nup*zp**2*D(c_edges_int) + num*zm**2*Dm)/T*c_edges_int*np.diff(phi_lyte)/dxd1)
-        Nm_edges_int = num*(-D(c_edges_int)*np.diff(c_lyte)/dxd1
-                            + (1./(num*zm)*(1-tp0(c_edges_int))*i_edges_int))
-
     elif ndD["elyteModelType"] == "SM":
         D_fs, sigma_fs, thermFac, tp0 = getattr(props_elyte,ndD["SMset"])()[:-1]
-        # modify the free solution transport properties for porous media
 
-        def D(c):
-            return eps_o_tau*D_fs(c)
+        # Get diffusivity and conductivity at cell edges using weighted harmonic mean
+        D_edges = utils.weighted_harmonic_mean(eps_o_tau*D_fs(c_lyte), wt)
+        sigma_edges = utils.weighted_harmonic_mean(eps_o_tau*sigma_fs(c_lyte), wt)
 
-        def sigma(c):
-            return eps_o_tau*sigma_fs(c)
         sp, n = ndD["sp"], ndD["n_refTrode"]
-        i_edges_int = -sigma(c_edges_int)/T * (
+        i_edges_int = -sigma_edges/T * (
             np.diff(phi_lyte)/dxd1
             + nu*T*(sp/(n*nup)+tp0(c_edges_int)/(zp*nup))
             * thermFac(c_edges_int)
             * np.diff(np.log(c_lyte))/dxd1
             )
-        Nm_edges_int = num*(-D(c_edges_int)*np.diff(c_lyte)/dxd1
+        Nm_edges_int = num*(-D_edges*np.diff(c_lyte)/dxd1
                             + (1./(num*zm)*(1-tp0(c_edges_int))*i_edges_int))
+    elif ndD["elyteModelType"] == "solid":
+        D_fs, sigma_fs, thermFac, tp0 = getattr(props_elyte, ndD["SMset"])()[:-1]
+
+        # Get diffusivity at cell edges using weighted harmonic mean
+        D_edges = utils.weighted_harmonic_mean(eps_o_tau * D_fs(c_lyte), wt)
+
+        # sp, n = ndD["sp"], ndD["n_refTrode"]
+        # D_fs is specified in solid_elyte_func in props_elyte.py
+        Dm = ndD["Dm"]
+        a_slyte = ndD["a_slyte"]
+        k = 1.381e-23
+
+        i_edges_int = (-((nup*zp*D_edges*(1/(1-c_edges_int)-a_slyte/(k*T)*2*c_edges_int)
+                          + num*zm*Dm)*np.diff(c_lyte)/dxd1)
+                       - (nup * zp ** 2 * D_edges + num * zm ** 2 * Dm) / T
+                       * c_edges_int * np.diff(phi_lyte) / dxd1)
+        Nm_edges_int = num * (-D_edges * np.diff(c_lyte) / dxd1
+                              + (1. / (num * zm) * (1 - tp0(c_edges_int)) * i_edges_int))
     return Nm_edges_int, i_edges_int
diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py
index c01e3967..7c9a6df7 100644
--- a/mpet/mod_electrodes.py
+++ b/mpet/mod_electrodes.py
@@ -14,7 +14,6 @@
 import daetools.pyDAE as dae
 import numpy as np
 import scipy.sparse as sprs
-import scipy.special as spcl
 
 import mpet.extern_funcs as extern_funcs
 import mpet.geometry as geo
@@ -22,7 +21,8 @@
 import mpet.props_am as props_am
 import mpet.utils as utils
 import mpet.electrode.reactions as reactions
-from mpet.daeVariableTypes import *
+from mpet.daeVariableTypes import mole_frac_t
+
 
 class Mod2var(dae.daeModel):
     def __init__(self, Name, Parent=None, Description="", ndD=None,
@@ -61,8 +61,8 @@ def __init__(self, Name, Parent=None, Description="", ndD=None,
             self.Rxn1 = dae.daeVariable("Rxn1", dae.no_t, self, "Rate of reaction 1", [self.Dmn])
             self.Rxn2 = dae.daeVariable("Rxn2", dae.no_t, self, "Rate of reaction 2", [self.Dmn])
 
-        #Get reaction rate function from dictionary name
-        self.calc_rxn_rate=getattr(reactions,ndD["rxnType"])
+        # Get reaction rate function from dictionary name
+        self.calc_rxn_rate = getattr(reactions,ndD["rxnType"])
 
         # Ports
         self.portInLyte = ports.portFromElyte(
@@ -127,13 +127,13 @@ def DeclareEquations(self):
             eq1.Residual -= self.c1(k) * volfrac_vec[k]
             eq2.Residual -= self.c2(k) * volfrac_vec[k]
         eq = self.CreateEquation("cbar")
-        eq.Residual = self.cbar() - utils.mean_linear(self.c1bar(), self.c2bar())
+        eq.Residual = self.cbar() - .5*(self.c1bar() + self.c2bar())
 
         # Define average rate of filling of particle
         eq = self.CreateEquation("dcbardt")
         eq.Residual = self.dcbardt()
         for k in range(N):
-            eq.Residual -= utils.mean_linear(self.c1.dt(k), self.c2.dt(k)) * volfrac_vec[k]
+            eq.Residual -= .5*(self.c1.dt(k) + self.c2.dt(k)) * volfrac_vec[k]
 
         c1 = np.empty(N, dtype=object)
         c2 = np.empty(N, dtype=object)
@@ -296,8 +296,8 @@ def __init__(self, Name, Parent=None, Description="", ndD=None,
         else:
             self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn])
 
-        #Get reaction rate function from dictionary name
-        self.calc_rxn_rate=getattr(reactions,ndD["rxnType"])
+        # Get reaction rate function from dictionary name
+        self.calc_rxn_rate = getattr(reactions,ndD["rxnType"])
 
         # Ports
         self.portInLyte = ports.portFromElyte(
@@ -340,8 +340,9 @@ def DeclareEquations(self):
                 for _position_ in range(N)]
 
         # Figure out mu_O, mu of the oxidized state
+        a_slyte = self.ndD_s.get("a_slyte", None)
         mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T,
-                                   self.ndD_s["a_slyte"], self.ndD_s["elyteModelType"])
+                                   self.ndD_s["elyteModelType"], a_slyte)
 
         # Define average filling fraction in particle
         eq = self.CreateEquation("cbar")
@@ -481,8 +482,7 @@ def calc_flux_diffn(c, D, Dfunc, Flux_bc, dr, T):
     Flux_vec = np.empty(N+1, dtype=object)
     Flux_vec[0] = 0  # Symmetry at r=0
     Flux_vec[-1] = Flux_bc
-    c_edges = utils.mean_harmonic(c)
-    Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.diff(c)/dr
+    Flux_vec[1:N] = -D * np.diff(c)/dr
     return Flux_vec
 
 
@@ -491,7 +491,7 @@ def calc_flux_CHR(c, mu, D, Dfunc, Flux_bc, dr, T):
     Flux_vec = np.empty(N+1, dtype=object)
     Flux_vec[0] = 0  # Symmetry at r=0
     Flux_vec[-1] = Flux_bc
-    c_edges = utils.mean_harmonic(c)
+    c_edges = utils.mean_linear(c)
     Flux_vec[1:N] = -D/T * Dfunc(c_edges) * np.diff(mu)/dr
     return Flux_vec
 
@@ -504,23 +504,24 @@ def calc_flux_CHR2(c1, c2, mu1_R, mu2_R, D, Dfunc, Flux1_bc, Flux2_bc, dr, T):
     Flux2_vec[0] = 0.  # symmetry at r=0
     Flux1_vec[-1] = Flux1_bc
     Flux2_vec[-1] = Flux2_bc
-    c1_edges = utils.mean_harmonic(c1)
-    c2_edges = utils.mean_harmonic(c2)
+    c1_edges = utils.mean_linear(c1)
+    c2_edges = utils.mean_linear(c2)
     Flux1_vec[1:N] = -D/T * Dfunc(c1_edges) * np.diff(mu1_R)/dr
     Flux2_vec[1:N] = -D/T * Dfunc(c2_edges) * np.diff(mu2_R)/dr
     return Flux1_vec, Flux2_vec
 
 
-def calc_mu_O(c_lyte, phi_lyte, phi_sld, T, a_slyte, elyteModelType):
+def calc_mu_O(c_lyte, phi_lyte, phi_sld, T, elyteModelType, a_slyte=None):
     if elyteModelType == "SM":
         mu_lyte = phi_lyte
         act_lyte = c_lyte
-    elif elyteModelType == "Solid":
-        act_lyte = c_lyte*np.exp(a_slyte*(1-c_lyte))
-        mu_lyte = T*np.log(act_lyte) + phi_lyte
     elif elyteModelType == "dilute":
         act_lyte = c_lyte
         mu_lyte = T*np.log(act_lyte) + phi_lyte
+    elif elyteModelType == "solid":
+        assert a_slyte is not None, "a_slyte must be specified when using solid elyte model"
+        act_lyte = c_lyte * np.exp(a_slyte * (1 - c_lyte))
+        mu_lyte = T * np.log(act_lyte) + phi_lyte
     mu_O = mu_lyte - phi_sld
     return mu_O, act_lyte
 
diff --git a/mpet/plot/outmat2txt.py b/mpet/plot/outmat2txt.py
index 39db9e40..a08d550e 100644
--- a/mpet/plot/outmat2txt.py
+++ b/mpet/plot/outmat2txt.py
@@ -77,15 +77,14 @@ def main(indir, genData=True, discData=True, elyteData=True,
     ndD_s, dD_s, ndD_e, dD_e = plot_data.show_data(
         indir, plot_type="params", print_flag=False, save_flag=False,
         data_only=True)
-    limtrode = ("c" if ndD_s["z"] < 1 else "a")
     trodes = ndD_s["trodes"]
-    CrateCurr = dD_e[limtrode]["cap"] / 3600.  # A/m^2
+    CrateCurr = dD_s["CrateCurr"]  # A/m^2
     psd_len_c = dD_s["psd_len"]["c"]
     Nv_c, Np_c = psd_len_c.shape
     dlm = ","
 
-    def get_trode_str(l):
-        return ("Anode" if l == "a" else "Cathode")
+    def get_trode_str(tr):
+        return ("Anode" if tr == "a" else "Cathode")
     if "a" in trodes:
         psd_len_a = dD_s["psd_len"]["a"]
         Nv_a, Np_a = psd_len_a.shape
@@ -142,15 +141,15 @@ def get_trode_str(l):
             print(file=fo)
             print(particleIndxExpl, file=fo)
             print(particleDiscExpl, file=fo)
-            for l in trodes:
+            for tr in trodes:
                 print(file=fo)
-                Trode = get_trode_str(l)
+                Trode = get_trode_str(tr)
                 print((Trode + " particle sizes [m]"), file=fo)
-                for vind in range(ndD_s["Nvol"][l]):
-                    print(",".join(map(str, dD_s["psd_len"][l][vind,:])), file=fo)
+                for vind in range(ndD_s["Nvol"][tr]):
+                    print(",".join(map(str, dD_s["psd_len"][tr][vind,:])), file=fo)
                 print("\n" + Trode + " particle number of discr. points", file=fo)
-                for vind in range(ndD_s["Nvol"][l]):
-                    print(",".join(map(str, ndD_s["psd_num"][l][vind,:])), file=fo)
+                for vind in range(ndD_s["Nvol"][tr]):
+                    print(",".join(map(str, ndD_s["psd_num"][tr][vind,:])), file=fo)
 
     if elyteData:
         valid_current = True
@@ -187,20 +186,20 @@ def get_trode_str(l):
         dataFileName = "output_data.mat"
         dataFile = os.path.join(indir, dataFileName)
         data = sio.loadmat(dataFile)
-        for l in trodes:
-            Trode = get_trode_str(l)
+        for tr in trodes:
+            Trode = get_trode_str(tr)
             type2c = False
-            if ndD_e[l]["type"] in ndD_s["1varTypes"]:
+            if ndD_e[tr]["type"] in ndD_s["1varTypes"]:
                 str_base = partStr + "c"
-            elif ndD_e[l]["type"] in ndD_s["2varTypes"]:
+            elif ndD_e[tr]["type"] in ndD_s["2varTypes"]:
                 type2c = True
                 str1_base = partStr + "c1"
                 str2_base = partStr + "c2"
-            for i in range(ndD_s["Npart"][l]):
-                for j in range(ndD_s["Nvol"][l]):
+            for i in range(ndD_s["Npart"][tr]):
+                for j in range(ndD_s["Nvol"][tr]):
                     if type2c:
-                        sol1 = str1_base.format(l=l, i=i, j=j)
-                        sol2 = str2_base.format(l=l, i=i, j=j)
+                        sol1 = str1_base.format(l=tr, i=i, j=j)
+                        sol2 = str2_base.format(l=tr, i=i, j=j)
                         datay1 = data[sol1]
                         datay2 = data[sol2]
                         filename1 = fnameSol1Base.format(l=Trode, i=i, j=j)
@@ -210,7 +209,7 @@ def get_trode_str(l):
                         np.savetxt(os.path.join(indir, filename2), datay2,
                                    delimiter=dlm, header=sol2Hdr)
                     else:
-                        sol = str_base.format(l=l, i=i, j=j)
+                        sol = str_base.format(l=tr, i=i, j=j)
                         datay = data[sol]
                         filename = fnameSolBase.format(l=Trode, i=i, j=j)
                         np.savetxt(os.path.join(indir, filename), datay,
@@ -220,10 +219,10 @@ def get_trode_str(l):
         cbarDict = plot_data.show_data(
             indir, plot_type="cbar_full", print_flag=False,
             save_flag=False, data_only=True)
-        for l in trodes:
-            Trode = get_trode_str(l)
+        for tr in trodes:
+            Trode = get_trode_str(tr)
             fname = "cbar{l}Data.txt".format(l=Trode)
-            Nv, Np = ndD_s["Nvol"][l], ndD_s["Npart"][l]
+            Nv, Np = ndD_s["Nvol"][tr], ndD_s["Npart"][tr]
             NpartTot = Nv*Np
             cbarMat = np.zeros((ntimes, NpartTot))
             cbarHdr = cbarHdrBase
@@ -231,7 +230,7 @@ def get_trode_str(l):
             for i in range(Nv):
                 for j in range(Np):
                     cbarHdr += "{i}/{j},".format(j=j, i=i)
-                    cbarMat[:,partInd] = cbarDict[l][:,i,j]
+                    cbarMat[:,partInd] = cbarDict[tr][:,i,j]
                     partInd += 1
             np.savetxt(os.path.join(indir, fname), cbarMat,
                        delimiter=dlm, header=cbarHdr.rstrip(','))
diff --git a/mpet/plot/plot_data.py b/mpet/plot/plot_data.py
index 79f69aa6..d6fb4f6c 100644
--- a/mpet/plot/plot_data.py
+++ b/mpet/plot/plot_data.py
@@ -23,7 +23,7 @@
 # mpl.rcParams['lines.linewidth'] = 3
 # mpl.rcParams['lines.markersize'] = 10
 # mpl.rcParams['lines.markeredgewidth'] = 0.1
-#mpl.rcParams['text.usetex'] = True
+# mpl.rcParams['text.usetex'] = True
 
 
 def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOut=None, tOut=None):
@@ -156,18 +156,18 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
             elif rxnType in ["Marcus", "MHC"]:
                 print("lambda_" + trode + "/(kTref):", ndD_e[trode]["lambda"])
             if ndD_s['simBulkCond'][trode]:
-                print(trode + " bulk conductivity loss: Yes -- " +
-                      "sigma_s [S/m]: " + str(dD_s['sigma_s'][trode]))
+                print(trode + " bulk conductivity loss: Yes -- "
+                      + "sigma_s [S/m]: " + str(dD_s['sigma_s'][trode]))
             else:
                 print(trode + " bulk conductivity loss: No")
             try:
                 simSurfCond = ndD_e[trode]['simSurfCond']
                 if simSurfCond:
-                    print(trode + " surface conductivity loss: Yes -- " +
-                          "dim_scond [S]: " + str(dD_e[trode]['scond']))
+                    print(trode + " surface conductivity loss: Yes -- "
+                          + "dim_scond [S]: " + str(dD_e[trode]['scond']))
                 else:
                     print(trode + " surface conductivity loss: No")
-            except:
+            except Exception:
                 pass
 #            if ndD['simSurfCond'][l]:
 #                print (l + " surface conductivity loss: Yes -- " +
@@ -182,8 +182,8 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
 
     # Plot voltage profile
     if plot_type in ["v", "vt"]:
-        voltage = (Vstd -
-                   (k*Tref/e)*data[pfx + 'phi_applied'][0])
+        voltage = (Vstd
+                   - (k*Tref/e)*data[pfx + 'phi_applied'][0])
         ffvec = data[pfx + 'ffrac_c'][0]
         fig, ax = plt.subplots(figsize=figsize)
         if plot_type == "v":
@@ -340,7 +340,7 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
             i_edges = np.zeros((numtimes, len(facesvec)))
             for tInd in range(numtimes):
                 i_edges[tInd, :] = mod_cell.get_lyte_internal_fluxes(
-                    cmat[tInd, :], pmat[tInd, :], disc["dxd1"], disc["eps_o_tau_edges"], ndD_s)[1]
+                    cmat[tInd, :], pmat[tInd, :], disc, ndD_s)[1]
             if plot_type in ["elytei", "elyteif"]:
                 ylbl = r'Current density of electrolyte [A/m$^2$]'
                 datax = facesvec
@@ -348,7 +348,7 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
             elif plot_type in ["elytedivi", "elytedivif"]:
                 ylbl = r'Divergence of electrolyte current density [A/m$^3$]'
                 datax = cellsvec
-                datay = np.diff(i_edges, axis=1) / disc["dxd2"]
+                datay = np.diff(i_edges, axis=1) / disc["dxvec"]
                 datay *= (F*dD_s["cref"]*dD_s["Dref"]/dD_s["Lref"]**2)
         if fplot:
             datay = datay[t0ind]
diff --git a/mpet/ports.py b/mpet/ports.py
index fc56ec79..6689fd1b 100644
--- a/mpet/ports.py
+++ b/mpet/ports.py
@@ -1,7 +1,8 @@
 """This defines the ports by which mod_cell interacts with mod_electrodes."""
 import daetools.pyDAE as dae
 
-from mpet.daeVariableTypes import *
+from mpet.daeVariableTypes import mole_frac_t, elec_pot_t
+
 
 class portFromElyte(dae.daePort):
     def __init__(self, Name, PortType, Model, Description=""):
diff --git a/mpet/props_am.py b/mpet/props_am.py
index 45911ca0..42f42aad 100644
--- a/mpet/props_am.py
+++ b/mpet/props_am.py
@@ -67,13 +67,12 @@ def __init__(self, T, ndD=None, **kwargs):
         self.eokT = e/(k*Tabs)
         self.kToe = (k*Tabs)/e
 
-        #Convert "muRfunc" to a callable function
+        # Convert "muRfunc" to a callable function
         self.muRfunc = getattr(self,ndD["muRfunc"])
 
     def get_muR_from_OCV(self, OCV, muR_ref):
         return -self.eokT*OCV + muR_ref
 
-
     ######
     # Solid solution functions
     # These are all obtained directly from fitting an OCV of the
@@ -84,10 +83,10 @@ def get_muR_from_OCV(self, OCV, muR_ref):
     def LiMn2O4_ss(self, y, ybar, muR_ref, ISfuncs=None):
         """ Doyle, Newman, 1996 """
         # OCV in V vs Li/Li+
-        OCV = (4.19829 + 0.0565661*np.tanh(-14.5546*y + 8.60942) -
-               0.0275479*(1/((0.998432 - y)**(0.492465)) - 1.90111) -
-               0.157123*np.exp(-0.04738*y**8) +
-               0.810239*np.exp(-40*(y - 0.133875)))
+        OCV = (4.19829 + 0.0565661*np.tanh(-14.5546*y + 8.60942)
+               - 0.0275479*(1/((0.998432 - y)**(0.492465)) - 1.90111)
+               - 0.157123*np.exp(-0.04738*y**8)
+               + 0.810239*np.exp(-40*(y - 0.133875)))
         muR = self.get_muR_from_OCV(OCV, muR_ref)
         actR = None
         return muR, actR
@@ -95,10 +94,10 @@ def LiMn2O4_ss(self, y, ybar, muR_ref, ISfuncs=None):
     def LiMn2O4_ss2(self, y, ybar, muR_ref, ISfuncs=None):
         """ Fuller, Doyle, Newman, 1994 """
         # OCV in V vs Li/Li+
-        OCV = (4.06279 + 0.0677504*np.tanh(-21.8502*y + 12.8268) -
-               0.105734*(1/((1.00167 - y)**(0.379571)) - 1.575994) -
-               0.045*np.exp(-71.69*y**8) +
-               0.01*np.exp(-200*(y - 0.19)))
+        OCV = (4.06279 + 0.0677504*np.tanh(-21.8502*y + 12.8268)
+               - 0.105734*(1/((1.00167 - y)**(0.379571)) - 1.575994)
+               - 0.045*np.exp(-71.69*y**8)
+               + 0.01*np.exp(-200*(y - 0.19)))
         muR = self.get_muR_from_OCV(OCV, muR_ref)
         actR = None
         return muR, actR
@@ -122,11 +121,11 @@ def LiC6_coke_ss2(self, y, ybar, muR_ref, ISfuncs=None):
 
     def LiC6_ss(self, y, ybar, muR_ref, ISfuncs=None):
         """ Safari, Delacourt 2011 """
-        OCV = (0.6379 + 0.5416*np.exp(-305.5309*y) +
-               0.044*np.tanh(-(y - 0.1958)/0.1088) -
-               0.1978*np.tanh((y - 1.0571)/0.0854) -
-               0.6875*np.tanh((y + 0.0117)/0.0529) -
-               0.0175*np.tanh((y - 0.5692)/0.0875))
+        OCV = (0.6379 + 0.5416*np.exp(-305.5309*y)
+               + 0.044*np.tanh(-(y - 0.1958)/0.1088)
+               - 0.1978*np.tanh((y - 1.0571)/0.0854)
+               - 0.6875*np.tanh((y + 0.0117)/0.0529)
+               - 0.0175*np.tanh((y - 0.5692)/0.0875))
         muR = self.get_muR_from_OCV(OCV, muR_ref)
         actR = None
         return muR, actR
@@ -176,8 +175,8 @@ def NCA_ss1(self, y, ybar, muR_ref, ISfuncs=None):
         This function was obtained from Dan Cogswell's fit of Samsung
         data.
         """
-        OCV = (3.86 + 1.67*y - 9.52*y**2 + 15.04*y**3 - 7.95*y**4 -
-               0.06*np.log(y/(1-y)))
+        OCV = (3.86 + 1.67*y - 9.52*y**2 + 15.04*y**3 - 7.95*y**4
+               - 0.06*np.log(y/(1-y)))
         muR = self.get_muR_from_OCV(OCV, muR_ref)
         actR = None
         return muR, actR
@@ -198,6 +197,18 @@ def NCA_ss2(self, y, ybar, muR_ref, ISfuncs=None):
         actR = None
         return muR, actR
 
+    def LTO(self, y, ybar, muR_ref, ISfuncs=None):
+        """
+        Vasileiadis 2017
+        """
+        muRtheta = -self.eokT * 1.55
+        muRhomog = self.reg_sln(y, self.ndD["Omga"], ISfuncs)
+        muRnonHomog = self.general_non_homog(y, ybar)
+        muR = muRhomog + muRnonHomog
+        actR = np.exp(muR / self.T)
+        muR += muRtheta + muR_ref
+        return muR, actR
+
     def testIS_ss(self, y, ybar, muR_ref, ISfuncs=None):
         """Ideal solution material for testing."""
         OCV = -self.kToe*np.log(y/(1-y))
@@ -401,16 +412,6 @@ def LiC6_1param(self, y, ybar, muR_ref, ISfuncs=None):
         muR += muRtheta + muR_ref
         return muR, actR
 
-    def LTO(self, y, ybar, muR_ref, ISfuncs=None):
-        """" Vasileiadis 2017 """
-        muRtheta = -self.eokT*1.55
-        muRhomog = self.reg_sln(y, self.ndD["Omga"], ISfuncs)
-        muRnonHomog = self.general_non_homog(y, ybar)
-        muR = muRhomog + muRnonHomog
-        actR = np.exp(muR/self.T)
-        muR += muRtheta + muR_ref
-        return muR, actR
-
     def testRS(self, y, ybar, muR_ref, ISfuncs=None):
         muRtheta = 0.
         muR = self.reg_sln(y, self.ndD["Omga"], ISfuncs)
@@ -427,6 +428,41 @@ def testRS_ps(self, y, ybar, muR_ref, ISfuncs=None):
         muR += muRtheta + muR_ref
         return muR, actR
 
+    def LiCoO2_LIONSIMBA(self, y, ybar, muR_ref, ISfuncs=None):
+        """ Torchio et al, 2016. """
+        r1 = 4.656
+        r2 = 88.669
+        r3 = 401.119
+        r4 = 342.909
+        r5 = 462.471
+        r6 = 433.434
+        r7 = 1
+        r8 = 18.933
+        r9 = 79.532
+        r10 = 37.311
+        r11 = 73.083
+        r12 = 95.96
+        OCV = (-r1 + r2*y**2 - r3*y**4 + r4*y**6 - r5*y**8 + r6 * y**10) \
+            / (-r7 + r8*y**2 - r9*y**4 + r10*y**6 - r11*y**8 + r12*y**10)
+        muR = self.get_muR_from_OCV(OCV, muR_ref)
+        actR = None
+        return muR, actR
+
+    def LiC6_LIONSIMBA(self, y, ybar, muR_ref, ISfuncs=None):
+        """ Torchio et al, 2016. """
+        r1 = 0.7222
+        r2 = 0.1387
+        r3 = 0.029
+        r4 = 0.0172
+        r5 = 0.0019
+        r6 = 0.2808
+        r7 = 0.7984
+        OCV = r1 + r2*y + r3*y**0.5 - r4*y**(-1) + r5*y**(-1.5) \
+            + r6*np.exp(0.9-15*y) - r7*np.exp(0.4465*y-0.4108)
+        muR = self.get_muR_from_OCV(OCV, muR_ref)
+        actR = None
+        return muR, actR
+
 
 def step_down(x, xc, delta):
     return 0.5*(-np.tanh((x - xc)/delta) + 1)
diff --git a/mpet/props_elyte.py b/mpet/props_elyte.py
index 110c06c7..e10a6467 100644
--- a/mpet/props_elyte.py
+++ b/mpet/props_elyte.py
@@ -1,4 +1,4 @@
-"""
+r"""
 This module provides functions defining properties of the ion-conducting
 phase -- the electrolyte Manage functions for the parameters involved in
 Stefan-Maxwell based concentrated electrolyte transport theory for
@@ -20,30 +20,6 @@
 k = 1.381e-23
 e = 1.602e-19
 
-def Solid_elyte_func():
-    "Solid Electrolyte version, several sources for different params"
-    def tp0(c): 
-        return 0.9
-    def D(c): 
-        return 1.19e-11 # m^2/s 
-    def kappa(c):
-        return 1.2e-3 # S/m
-    def kappa_valoen_reimers(c):
-        (k00, k01, k02,
-                k10, k11, k12,
-                k20, k21) = (
-                -10.5, 0.0740, -6.96e-5,
-                0.668, -0.0178, 2.80e-5,
-                0.494, -8.86e-4)
-        out = c * (k00 + k01*Tref + k02*Tref**2
-                + k10*c + k11*c*Tref + k12*c*Tref**2
-                + k20*c**2 + k21*c**2*Tref)**2 # mS/cm
-        out *= 0.1 # S/m
-        return out
-    Ign1, sigma_ndim, thermFac, Ign3, Dref = valoen_reimers() 
-    D_ndim = lambda c: D(c)/Dref 
-    kappa_ndim = lambda c: kappa(c)/kappa_valoen_reimers(c)
-    return D_ndim, sigma_ndim, thermFac, tp0, Dref 
 
 def LiClO4_PC():
     """ Set of parameters from Fuller, Doyle, Newman 1994, with
@@ -143,6 +119,35 @@ def sigma_ndim(c):
     return D_ndim, sigma_ndim, therm_fac, tp0, Dref
 
 
+def solid_elyte():
+    """
+    Solid Electrolyte version, several sources for different params
+    """
+    # LCO: kappa_ndim is not actually returned, so unclear
+    # what it should be used for
+    # related functions and values commented out for now
+
+    tp0 = 0.9  # tp0 is constant but a callable must be returned
+    D = 1.19e-11  # m^2/s
+    # kappa = 1.2e-3  # S/m
+
+    # def kappa_valoen_reimers(c):
+    #     k00, k01, k02, k10, k11, k12, k20, k21 = (-10.5, 0.0740, -6.96e-5, 0.668,
+    #                                               -0.0178, 2.80e-5, 0.494, -8.86e-4)
+    #     out = c * (k00 + k01*Tref + k02*Tref**2
+    #                + k10*c + k11*c*Tref + k12*c*Tref**2
+    #                + k20*c**2 + k21*c**2*Tref)**2  # mS/cm
+    #     out *= 0.1  # S/m
+    #     return out
+
+    Ign1, sigma_ndim, thermFac, Ign3, Dref = valoen_reimers()
+    D_ndim = D / Dref
+    # kappa_ndim = lambda c: kappa / kappa_valoen_reimers(c)
+
+    # D_ndim and tp0 are constants, but a callable must be returned
+    return lambda c: D_ndim, sigma_ndim, thermFac, lambda c: tp0, Dref
+
+
 def test1():
     """Set of dilute solution parameters with zp=|zm|=nup=num=1,
     Dp = 2.2e-10 m^2/s
@@ -171,3 +176,81 @@ def sigma_ndim(c):
         return sigma(c) * (
             k*Tref/(e**2*Dref*N_A*(1000*cref)))
     return D_ndim, sigma_ndim, therm_fac, tp0, Dref
+
+
+def LIONSIMBA():
+    """ Set of parameters from LIONSIMBA validation. Torchio et al, 2016.
+    """
+    T = 298  # isothermal model
+
+    def tp0(c):
+        return 0.364
+
+    def sigma(c):
+        c_dim = c*1000  # dimensionalized c
+        r1 = -10.5
+        r2 = 0.668e-3
+        r3 = 0.494e-6
+        r4 = 0.074
+        r5 = -1.78e-5
+        r6 = -8.86e-10
+        r7 = -6.96e-5
+        r8 = 2.8e-8
+        sig_out = 1e-4 * c_dim * (r1 + r2*c_dim + r3*c_dim**2
+                                  + T*(r4 + r5*c_dim + r6*c_dim**2)
+                                  + T**2 * (r7 + r8*c_dim))**2
+        return sig_out  # m^2/s
+
+    def D(c):
+        c_dim = c*1000
+        T = 298
+        r1 = 4.43
+        r2 = 54
+        r3 = 229
+        r4 = 5e-3
+        r5 = 0.22e-3
+        D_out = 1e-4 * 10**(-r1-r2/(T-r3-r4*c_dim)-r5*c_dim)
+        return D_out
+
+    def therm_fac(c):
+        return 1.
+
+    Dref = D(cref)
+
+    def D_ndim(c):
+        return D(c) / Dref
+
+    def sigma_ndim(c):
+        return sigma(c) * (
+            k*Tref/(e**2*Dref*N_A*(1000*cref)))
+    return D_ndim, sigma_ndim, therm_fac, tp0, Dref
+
+
+def LIONSIMBA_isothermal():
+    """ Set of parameters from LIONSIMBA validation. Torchio et al, 2016.
+    """
+    # isothermal at 298 K
+
+    def tp0(c):
+        return 0.364
+
+    def sigma(c):
+        ce = c*1000  # dimensionalized c
+        return 4.1253e-2 + 5.007e-4*ce - 4.7212e-7*ce**2 \
+            + 1.5094e-10*ce**3 - 1.6018*1e-14*ce**4  # S/m
+
+    def D(c):
+        return 7.5e-10  # m^2/s
+
+    def therm_fac(c):
+        return 1.
+
+    Dref = D(cref)
+
+    def D_ndim(c):
+        return D(c) / Dref
+
+    def sigma_ndim(c):
+        return sigma(c) * (
+            k*Tref/(e**2*Dref*N_A*(1000*cref)))
+    return D_ndim, sigma_ndim, therm_fac, tp0, Dref
diff --git a/mpet/sim.py b/mpet/sim.py
index 55cc2f84..28e6097d 100644
--- a/mpet/sim.py
+++ b/mpet/sim.py
@@ -31,10 +31,10 @@ def __init__(self, ndD_s=None, ndD_e=None, tScale=None):
             ndD_s["currPrev"] = self.dataPrev["current"][0,-1]
             ndD_s["phiPrev"] = self.dataPrev["phi_applied"][0,-1]
 
-        #Set absolute tolerances for variableTypes
-        mpet.daeVariableTypes.mole_frac_t.AbsoluteTolerance=ndD_s["absTol"]
-        mpet.daeVariableTypes.conc_t.AbsoluteTolerance=ndD_s["absTol"]
-        mpet.daeVariableTypes.elec_pot_t.AbsoluteTolerance=ndD_s["absTol"]
+        # Set absolute tolerances for variableTypes
+        mpet.daeVariableTypes.mole_frac_t.AbsoluteTolerance = ndD_s["absTol"]
+        mpet.daeVariableTypes.conc_t.AbsoluteTolerance = ndD_s["absTol"]
+        mpet.daeVariableTypes.elec_pot_t.AbsoluteTolerance = ndD_s["absTol"]
 
         # Define the model we're going to simulate
         self.m = mod_cell.ModCell("mpet", ndD_s=ndD_s, ndD_e=ndD_e)
@@ -44,13 +44,13 @@ def SetUpParametersAndDomains(self):
         ndD = self.ndD_s
         if ndD["Nvol"]["s"] >= 1:
             self.m.DmnCell["s"].CreateArray(ndD["Nvol"]["s"])
-        for l in ndD["trodes"]:
-            self.m.DmnCell[l].CreateArray(ndD["Nvol"][l])
-            self.m.DmnPart[l].CreateArray(ndD["Npart"][l])
-            for i in range(ndD["Nvol"][l]):
-                for j in range(ndD["Npart"][l]):
-                    self.m.particles[l][i, j].Dmn.CreateArray(
-                        int(ndD["psd_num"][l][i,j]))
+        for tr in ndD["trodes"]:
+            self.m.DmnCell[tr].CreateArray(ndD["Nvol"][tr])
+            self.m.DmnPart[tr].CreateArray(ndD["Npart"][tr])
+            for i in range(ndD["Nvol"][tr]):
+                for j in range(ndD["Npart"][tr]):
+                    self.m.particles[tr][i, j].Dmn.CreateArray(
+                        int(ndD["psd_num"][tr][i,j]))
 
     def SetUpVariables(self):
         ndD_s = self.ndD_s
@@ -59,26 +59,26 @@ def SetUpVariables(self):
         phi_cathode = ndD_s["phi_cathode"]
         if ndD_s["prevDir"] == "false":
             # Solids
-            for l in ndD_s["trodes"]:
-                cs0 = self.ndD_s['cs0'][l]
+            for tr in ndD_s["trodes"]:
+                cs0 = self.ndD_s['cs0'][tr]
                 # Guess initial filling fractions
-                self.m.ffrac[l].SetInitialGuess(cs0)
-                for i in range(Nvol[l]):
+                self.m.ffrac[tr].SetInitialGuess(cs0)
+                for i in range(Nvol[tr]):
                     # Guess initial volumetric reaction rates
-                    self.m.R_Vp[l].SetInitialGuess(i, 0.0)
+                    self.m.R_Vp[tr].SetInitialGuess(i, 0.0)
                     # Guess initial value for the potential of the
                     # electrodes
-                    if l == "a":  # anode
-                        self.m.phi_bulk[l].SetInitialGuess(i, self.ndD_s["phiRef"]["a"])
+                    if tr == "a":  # anode
+                        self.m.phi_bulk[tr].SetInitialGuess(i, self.ndD_s["phiRef"]["a"])
                     else:  # cathode
-                        self.m.phi_bulk[l].SetInitialGuess(i, phi_cathode)
-                    for j in range(Npart[l]):
-                        Nij = ndD_s["psd_num"][l][i,j]
-                        part = self.m.particles[l][i,j]
+                        self.m.phi_bulk[tr].SetInitialGuess(i, phi_cathode)
+                    for j in range(Npart[tr]):
+                        Nij = ndD_s["psd_num"][tr][i,j]
+                        part = self.m.particles[tr][i,j]
                         # Guess initial value for the average solid
                         # concentrations and set initial value for
                         # solid concentrations
-                        solidType = self.ndD_e[l]["indvPart"][i,j]["type"]
+                        solidType = self.ndD_e[tr]["indvPart"][i,j]["type"]
                         if solidType in ndD_s["1varTypes"]:
                             part.cbar.SetInitialGuess(cs0)
                             for k in range(Nij):
@@ -95,60 +95,60 @@ def SetUpVariables(self):
                             for k in range(Nij):
                                 part.c1.SetInitialCondition(k, cs0+rnd1[k])
                                 part.c2.SetInitialCondition(k, cs0+rnd2[k])
-            
-            #Cell potential initialization
-            if ndD_s['tramp']>0:
-                phi_guess=0
-            elif ndD_s['profileType']=='CV':
+
+            # Cell potential initialization
+            if ndD_s['tramp'] > 0:
+                phi_guess = 0
+            elif ndD_s['profileType'] == 'CV':
                 phi_guess = self.ndD_s['Vset']
-            elif ndD_s['profileType']=='CVsegments':
+            elif ndD_s['profileType'] == 'CVsegments':
                 phi_guess = self.ndD_s['segments'][0][0]
             else:
                 phi_guess = 0
             self.m.phi_applied.SetInitialGuess(phi_guess)
             self.m.phi_cell.SetInitialGuess(phi_guess)
 
-            #Initialize the ghost points used for boundary conditions
+            # Initialize the ghost points used for boundary conditions
             if not self.m.SVsim:
                 self.m.c_lyteGP_L.SetInitialGuess(ndD_s["c0"])
                 self.m.phi_lyteGP_L.SetInitialGuess(0)
 
-            #Separator electrolyte initialization
+            # Separator electrolyte initialization
             for i in range(Nvol["s"]):
                 self.m.c_lyte["s"].SetInitialCondition(i, ndD_s['c0'])
                 self.m.phi_lyte["s"].SetInitialGuess(i, 0)
-            
-            #Anode and cathode electrolyte initialization
-            for l in ndD_s["trodes"]:
-                for i in range(Nvol[l]):
-                    self.m.c_lyte[l].SetInitialCondition(i, ndD_s['c0'])
-                    self.m.phi_lyte[l].SetInitialGuess(i, 0)
-
-                    #Set electrolyte concentration in each particle
-                    for j in range(Npart[l]):
-                        self.m.particles[l][i,j].c_lyte.SetInitialGuess(ndD_s["c0"])
-                        
+
+            # Anode and cathode electrolyte initialization
+            for tr in ndD_s["trodes"]:
+                for i in range(Nvol[tr]):
+                    self.m.c_lyte[tr].SetInitialCondition(i, ndD_s['c0'])
+                    self.m.phi_lyte[tr].SetInitialGuess(i, 0)
+
+                    # Set electrolyte concentration in each particle
+                    for j in range(Npart[tr]):
+                        self.m.particles[tr][i,j].c_lyte.SetInitialGuess(ndD_s["c0"])
+
         else:
             dPrev = self.dataPrev
-            for l in ndD_s["trodes"]:
-                self.m.ffrac[l].SetInitialGuess(
-                    dPrev["ffrac_" + l][0,-1])
-                for i in range(Nvol[l]):
-                    self.m.R_Vp[l].SetInitialGuess(
-                        i, dPrev["R_Vp_" + l][-1,i])
-                    self.m.phi_bulk[l].SetInitialGuess(
-                        i, dPrev["phi_bulk_" + l][-1,i])
-                    for j in range(Npart[l]):
-                        Nij = ndD_s["psd_num"][l][i,j]
-                        part = self.m.particles[l][i,j]
-                        solidType = self.ndD_e[l]["indvPart"][i, j]["type"]
+            for tr in ndD_s["trodes"]:
+                self.m.ffrac[tr].SetInitialGuess(
+                    dPrev["ffrac_" + tr][0,-1])
+                for i in range(Nvol[tr]):
+                    self.m.R_Vp[tr].SetInitialGuess(
+                        i, dPrev["R_Vp_" + tr][-1,i])
+                    self.m.phi_bulk[tr].SetInitialGuess(
+                        i, dPrev["phi_bulk_" + tr][-1,i])
+                    for j in range(Npart[tr]):
+                        Nij = ndD_s["psd_num"][tr][i,j]
+                        part = self.m.particles[tr][i,j]
+                        solidType = self.ndD_e[tr]["indvPart"][i, j]["type"]
                         partStr = "partTrode{l}vol{i}part{j}_".format(
-                            l=l, i=i, j=j)
-                        
-                        #Set the inlet port variables for each particle
-                        part.c_lyte.SetInitialGuess(dPrev[partStr+"portInLyte_c_lyte"][0,-1])
-                        part.phi_lyte.SetInitialGuess(dPrev[partStr+"portInLyte_phi_lyte"][0,-1])
-                        part.phi_m.SetInitialGuess(dPrev[partStr+"portInBulk_phi_m"][0,-1])
+                            l=tr, i=i, j=j)
+
+                        # Set the inlet port variables for each particle
+                        part.c_lyte.SetInitialGuess(dPrev["c_lyte_" + tr][-1,i])
+                        part.phi_lyte.SetInitialGuess(dPrev["phi_lyte_" + tr][-1,i])
+                        part.phi_m.SetInitialGuess(dPrev["phi_part_" + tr][-1][i,j])
 
                         if solidType in ndD_s["1varTypes"]:
                             part.cbar.SetInitialGuess(
@@ -173,23 +173,23 @@ def SetUpVariables(self):
                     i, dPrev["c_lyte_s"][-1,i])
                 self.m.phi_lyte["s"].SetInitialGuess(
                     i, dPrev["phi_lyte_s"][-1,i])
-            for l in ndD_s["trodes"]:
-                for i in range(Nvol[l]):
-                    self.m.c_lyte[l].SetInitialCondition(
-                        i, dPrev["c_lyte_" + l][-1,i])
-                    self.m.phi_lyte[l].SetInitialGuess(
-                        i, dPrev["phi_lyte_" + l][-1,i])
-            
-            #Read in the ghost point values
+            for tr in ndD_s["trodes"]:
+                for i in range(Nvol[tr]):
+                    self.m.c_lyte[tr].SetInitialCondition(
+                        i, dPrev["c_lyte_" + tr][-1,i])
+                    self.m.phi_lyte[tr].SetInitialGuess(
+                        i, dPrev["phi_lyte_" + tr][-1,i])
+
+            # Read in the ghost point values
             if not self.m.SVsim:
                 self.m.c_lyteGP_L.SetInitialGuess(dPrev["c_lyteGP_L"][0,-1])
                 self.m.phi_lyteGP_L.SetInitialGuess(dPrev["phi_lyteGP_L"][0,-1])
-            
+
             # Guess the initial cell voltage
             self.m.phi_applied.SetInitialGuess(
                 dPrev["phi_applied"][0,-1])
 
-        #The simulation runs when the endCondition is 0
+        # The simulation runs when the endCondition is 0
         self.m.endCondition.AssignValue(0)
 
     def Run(self):
@@ -197,19 +197,17 @@ def Run(self):
         Overload the simulation "Run" function so that the simulation
         terminates when the specified condition is satisfied.
         """
-        time = 0.
         tScale = self.tScale
         for nextTime in self.ReportingTimes:
             self.Log.Message(
                 "Integrating from {t0:.2f} to {t1:.2f} s ...".format(
                     t0=self.CurrentTime*tScale, t1=nextTime*tScale), 0)
-            time = self.IntegrateUntilTime(
-                nextTime, dae.eStopAtModelDiscontinuity, True)
+            self.IntegrateUntilTime(nextTime, dae.eStopAtModelDiscontinuity, True)
             self.ReportData(self.CurrentTime)
             self.Log.SetProgress(int(100. * self.CurrentTime/self.TimeHorizon))
 
-            #Break when an end condition has been met
+            # Break when an end condition has been met
             if self.m.endCondition.npyValues:
-                description=mod_cell.endConditions[int(self.m.endCondition.npyValues)]
+                description = mod_cell.endConditions[int(self.m.endCondition.npyValues)]
                 self.Log.Message("Ending condition: " + description, 0)
-                break
\ No newline at end of file
+                break
diff --git a/mpet/utils.py b/mpet/utils.py
index 84f0a4fa..33cfb6ae 100644
--- a/mpet/utils.py
+++ b/mpet/utils.py
@@ -5,20 +5,22 @@
 import daetools.pyDAE as dae
 
 
-def mean_linear(a, b=None):
-    """Calculate the linear mean along a vector or between two values."""
-    if isinstance(a, np.ndarray):
-        return 0.5*(a[1:] + a[:-1])
-    else:
-        return 0.5*(a + b)
+def mean_linear(a):
+    """Calculate the linear mean along a vector."""
+    return 0.5*(a[1:] + a[:-1])
 
 
-def mean_harmonic(a, b=None):
-    """Calculate the harmonic mean along a vector or between two values."""
-    if isinstance(a, np.ndarray):
-        return (2 * a[1:] * a[:-1]) / (a[1:] + a[:-1] + 1e-20)
-    else:
-        return (2 * a * b) / (a + b + 1e-20)
+def weighted_linear_mean(a, wt):
+    return (wt[1:]*a[1:] + wt[:-1]*a[:-1])/(wt[1:]+wt[:-1])
+
+
+def mean_harmonic(a):
+    """Calculate the harmonic mean along a vector."""
+    return (2 * a[1:] * a[:-1]) / (a[1:] + a[:-1] + 1e-20)
+
+
+def weighted_harmonic_mean(a, wt):
+    return((wt[1:]+wt[:-1])/(wt[1:]/a[1:]+wt[:-1]/a[:-1]))
 
 
 def get_cell_Ntot(Nvol):
diff --git a/mpet/version.py b/mpet/version.py
index bad32e90..66a87bb6 100644
--- a/mpet/version.py
+++ b/mpet/version.py
@@ -1,1 +1,1 @@
-__version__ = '0.1.4'
\ No newline at end of file
+__version__ = '0.1.5'
diff --git a/setup.py b/setup.py
index 332455cf..3d2889ca 100644
--- a/setup.py
+++ b/setup.py
@@ -1,20 +1,29 @@
 import setuptools
 
-exec(open('mpet/version.py').read())
+version = {}
+with open('mpet/version.py', 'r', encoding='utf-8') as fh:
+    exec(fh.read(), version)
+__version__ = version['__version__']
+
+with open('README.md', 'r', encoding='utf-8') as fh:
+    long_description = fh.read()
 
 setuptools.setup(
     name='mpet',
     version=__version__,
     description='Multiphase porous electrode theory',
+    long_description=long_description,
+    long_description_content_type="text/markdown",
     author='Dan Cogswell',
     author_email='cogswell@mit.edu',
     license='MIT',
     url='https://bitbucket.org/bazantgroup/mpet',
     packages=['mpet','mpet.plot','mpet.electrode'],
-    install_requires=["numpy","scipy","matplotlib","pyQt5"],
+    install_requires=['numpy','scipy','matplotlib','pyQt5'],
+    extras_require={'test':['pytest','coverage', 'coveralls','configparser','h5py','flake8']},
     python_requires='>=3.5,<3.8',
     scripts=['bin/mpetrun.py','bin/mpetplot.py'],
     classifiers=[
         "Programming Language :: Python :: 3",
     ],
-)
\ No newline at end of file
+)


Test patch:
diff --git a/.github/workflows/mpet-docker-regression-test.yml b/.github/workflows/mpet-docker-regression-test.yml
new file mode 100644
index 00000000..1e860a1b
--- /dev/null
+++ b/.github/workflows/mpet-docker-regression-test.yml
@@ -0,0 +1,62 @@
+name: MPET regression test
+
+on: [push, workflow_dispatch]
+
+jobs:
+  test:
+    
+    runs-on: ubuntu-latest
+    container:
+      image: vikko00/daetools:1.9.0
+
+    steps:
+    - name: update container
+      run: |
+        apt-get update
+        apt-get upgrade -y
+
+    - uses: actions/checkout@v1
+      with:
+        fetch-depth: 1
+    - uses: actions/checkout@v2
+      with:
+        fetch-depth: 1
+        path: mpet
+    - uses: actions/checkout@v2
+      with:
+        fetch-depth: 1
+        path: base-mpet
+        ref: ${{ github.base_ref }}
+
+    - name: Install additional dependencies using mpet's setup.py
+      run: |
+        cd mpet
+        pip install .[test]
+
+    - name: Set up test for modified branch
+      run: |
+        cd mpet/bin
+        rm -rf workdir
+        mkdir workdir
+        cd workdir
+
+        cp ../run_tests.py .
+        ln -s ../../mpet .
+        ln -s ../../tests .
+
+    - name: run tests for modified branch and get coverage
+      run: |
+        cd mpet/bin/workdir
+        coverage run --source=../../mpet/ run_tests.py --test_dir ./tests --output_dir ../../bin/workdir/modified > /dev/null
+
+    - name: upload Coveralls
+      env:
+        GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+      run: |
+        cd mpet/bin/workdir
+        coveralls --service=github || : #Dont fret if if fails
+
+    - name: Checks test results
+      run: |
+        cd mpet/tests
+        pytest --baseDir=ref_outputs --modDir=../bin/workdir/modified compare_tests.py
diff --git a/bin/mpettest.py b/bin/mpettest.py
index 97aa9977..8e810e71 100755
--- a/bin/mpettest.py
+++ b/bin/mpettest.py
@@ -1,18 +1,23 @@
 #!/usr/bin/env python3
 
-import os
 import os.path as osp
-import sys
-import time
+import pytest
+import run_tests
 
-import tests.test_suite as tst
+# Run the tests and save the input args
+args = run_tests.main()
 
-if len(sys.argv) > 1:
-    compareDir = osp.join(os.getcwd(), sys.argv[1])
-else:
-    compareDir = None
-timeStart = time.time()
-tst.main(compareDir)
-timeEnd = time.time()
-tTot = timeEnd - timeStart
-print("Total test time:", tTot, "s")
+# Create the list of arguments for pytest
+pytest_args = ["--baseDir="+osp.join(args.test_dir,"ref_outputs"),
+               "--modDir="+args.output_dir,
+               osp.join(osp.dirname(osp.abspath(__file__)),"../tests/compare_tests.py")
+               ]
+
+# Add addtional options if a test list was provided
+if args.tests:
+    pytest_args.append("--tests="+" ".join(args.tests))
+    pytest_args.append("--skip-analytic")
+
+
+# Run pytest
+pytest.main(pytest_args)
diff --git a/bin/run_tests.py b/bin/run_tests.py
new file mode 100644
index 00000000..0672d932
--- /dev/null
+++ b/bin/run_tests.py
@@ -0,0 +1,64 @@
+import os.path as osp
+import time
+from os import getcwd, makedirs, walk
+from tests.test_suite import run_test_sims, run_test_sims_analyt
+import tests.test_defs as defs
+from shutil import rmtree
+import argparse
+
+
+def run(test_outputs, testDir, tests=None):
+    pflag = False
+    dirDict = {}
+    # Get the default configs
+    dirDict["suite"] = testDir
+    dirDict["simOut"] = osp.join(getcwd(), "sim_output")
+    dirDict["out"] = test_outputs
+    dirDict["baseConfig"] = osp.join(dirDict["suite"], "baseConfigs")
+
+    # Dictionary containing info about the tests to run
+    # Identifier strings are associated with functions to call and
+    # whether to run that particular test.
+    if not tests:
+        # Generate a list of tests from the directories in ref_outputs
+        ref_outputs = osp.join(dirDict["suite"],"ref_outputs")
+        _, directories, _ = next(walk(osp.join(ref_outputs)))
+        tests = directories
+        tests.sort()
+
+        runInfoAnalyt = {
+            "testAnalytCylDifn": (defs.testAnalytCylDifn, defs.analytCylDifn),
+            "testAnalytSphDifn": (defs.testAnalytSphDifn, defs.analytSphDifn),
+        }
+
+    if osp.exists(dirDict["out"]):
+        rmtree(dirDict["out"])
+    makedirs(dirDict["out"])
+    dirDict["plots"] = osp.join(dirDict["out"], "plots")
+    makedirs(dirDict["plots"])
+    run_test_sims(tests, dirDict, pflag)
+    try:
+        run_test_sims_analyt(runInfoAnalyt, dirDict)
+    except Exception:
+        pass
+
+
+def main():
+    parser = argparse.ArgumentParser(description='Run test suite.')
+    parser.add_argument('--output_dir', metavar='o', type=str,
+                        default=osp.join(osp.dirname(osp.abspath(__file__)),"../tests",
+                                         "test_outputs", time.strftime("%Y%m%d_%H%M%S")),
+                        help='where to put the simulations, default is tests/tests_outputs/date')
+    parser.add_argument('--test_dir', metavar='t', type=str,
+                        default=osp.join(osp.dirname(osp.abspath(__file__)),"../tests"),
+                        help='where are the tests located?')
+    parser.add_argument('tests', nargs='*', default=[], help='which tests do I run?')
+    args = parser.parse_args()
+
+    run(args.output_dir,args.test_dir, args.tests)
+
+    return(args)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/README.md b/tests/README.md
index 77438196..d7660f35 100644
--- a/tests/README.md
+++ b/tests/README.md
@@ -1,7 +1,58 @@
-## Testing
+# Regression Testing
 
 When adding new features or making changes to the code, it's helpful to run a suite of tests to make sure various things are behaving as expected. This should not be necessary for users who are not changing the code at all, although it could still be nice to verify that the outputs users are seeing match those the developers expect for a few specific cases.
 
-To run the tests, copy `mpettest.py` to the repository root directory and run it. This will run a number of simulations and compare their outputs to those outputs from the downloaded reference outputs along with a few comparisons of simulations to their corresponding analytical results.
+## Prerequisites
 
-Note that tests may "fail" even when things are okay, resulting from small numerical differences. If tests fail, it is helpful to look at the comparison plots generated by default within `mpet/tests/test_outputs/[time-stamped-directory]/plots` to see if the differences seem significant.
\ No newline at end of file
+The testing suite has additional dependencies, which can be installed from MPET's root directory with `pip install .[test]`. These additional dependencies are: `pytest`, `coverage`, `coveralls`, `configparser`, and `h5py`.
+
+## Running
+
+To run the tests, execute `PYTHONPATH=. ./bin/mpettest.py` from the repository root. The `-h` flag shows which arguments are accepted.  This will run a number of
+simulations and test the results.
+
+To compare the output manually you can use pytest:
+```bash
+  pytest --baseDir=tests/ref_outputs/ --modDir=tests/test_outputs/20201208_154137/ tests/compare_tests.py
+```
+or you can compare specific tests:
+```bash
+  pytest --baseDir=tests/ref_outputs/ --modDir=tests/test_outputs/20201208_154137/ tests/compare_tests.py --tests test001 test002 --skip-analytic
+```
+or you can compare timings:
+```bash
+  pytest --baseDir=tests/ref_outputs/ --modDir=tests/test_outputs/20201208_154137/ tests/compare_timings.py --tests test001 test002 --skip-analytic
+```
+
+You can also compare different output folders, or against the reference solution.
+
+Note that tests may "fail" even when things are okay, resulting from small numerical differences. If tests fail, it is helpful to look at the comparison plots generated by default within `mpet/tests/test_outputs/[time-stamped-directory]/plots` to see if the differences seem significant.
+There is also the flag --tolerance which can be used in combination with pytest
+it is both tested as absolute and relative tolerance.
+
+```bash
+  pytest --tolerance=0.01 --baseDir=tests/ref_outputs/ --modDir=tests/test_outputs/20201208_154137/ tests/compare_tests.py
+```
+
+
+# List of tests
+
+ - test001: LFP ACR C3
+ - test002: LFP CHR cylinder
+ - test003: LFP CHR sphere
+ - test004: LFP CHR sphere with noise
+ - test005: LFP homog
+ - test006: LFP homog with logPad, Vmin
+ - test007: LFP homog with logPad, Vmax
+ - test008: LFP homog_sdn
+ - test009: Graphite-2param homog
+ - test010: Graphite-2param CHR cylinder
+ - test011: Graphite-2param CHR sphere
+ - test012: Solid solution, diffn sphere, homog, LiC6_coke_ss2, LiMn2O4_ss2, BV_mod01, BV_mod02, cathode + anode
+ - test013: Solid solution, diffn cylinder, homog, testIS_ss, LiMn2O4_ss2, Marcus, BV_raw, cathode + separator + anode
+ - test014: LFP homog with CCsegments, MHC, Rser
+ - test015: testRS homog with CVsegments, bulkCond, partCond
+ - test016: test CC continuation
+ - test017: test CV continuation
+ - test018: Like test014, LFP homog with CCsegments, BV, Rfilm, Rfilm_foil
+ - test019: SM electrolyte with LFP homog, sep + cathode, significant elyte polarization
diff --git a/tests/compare_plots.py b/tests/compare_plots.py
index dff111a3..54acf262 100644
--- a/tests/compare_plots.py
+++ b/tests/compare_plots.py
@@ -1,6 +1,5 @@
 import os.path as osp
 
-import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/tests/compare_tests.py b/tests/compare_tests.py
new file mode 100644
index 00000000..7e05fa73
--- /dev/null
+++ b/tests/compare_tests.py
@@ -0,0 +1,127 @@
+import os.path as osp
+import scipy.io as sio
+import numpy as np
+import mpet.io_utils as IO
+import tests.test_defs as defs
+import errno
+import h5py
+import pytest
+
+
+def test_compare(Dirs, tol):
+    refDir, testDir = Dirs
+    newDir = osp.join(testDir, "sim_output")
+    refDir = osp.join(refDir, "sim_output")
+    newDatah5 = False
+    refDatah5 = False
+    newDataFile = osp.join(newDir, "output_data.mat")
+    if not osp.exists(newDataFile):
+        newDataFile = osp.join(newDir, "output_data.hdf5")
+        newDatah5 = True
+    assert osp.exists(newDataFile), "neither output_data.{mat,hdf5} present"
+    refDataFile = osp.join(refDir, "output_data.mat")
+    if not osp.exists(refDataFile):
+        refDataFile = osp.join(refDir, "output_data.hdf5")
+        refDatah5 = True
+    assert osp.exists(refDataFile), "neither output_data.{mat,hdf5} present"
+
+    try:
+        if newDatah5:
+            newData = h5py.File(newDataFile, 'r')
+        else:
+            newData = sio.loadmat(newDataFile)
+    except IOError as exception:
+        # If it's an error _other than_ the file not being there
+        assert exception.errno == errno.ENOENT, "IO error on opening file"
+        assert False, "File %s does not exist" % (newDataFile)
+
+    if refDatah5:
+        refData = h5py.File(refDataFile, 'r')
+    else:
+        refData = sio.loadmat(refDataFile)
+    for varKey in (set(refData.keys()) & set(newData.keys())):
+        # TODO -- Consider keeping a list of the variables that fail
+
+        # Ignore certain entries not of numerical output
+        if varKey[0:2] == "__":
+            continue
+
+        # Compute the difference between the solution and the reference
+        try:
+            varDataNew = newData[varKey]
+            varDataRef = refData[varKey]
+            diffMat = np.abs(varDataNew - varDataRef)
+        except ValueError:
+            assert False, "Fail from ValueError"
+        except KeyError:
+            assert False, "Fail from KeyError"
+
+        # #Check absolute and relative error against tol
+        assert np.max(diffMat) < tol or \
+               np.max(diffMat) < tol * np.max(np.abs(varDataRef)), \
+               "Fail from tolerance\nVariable failing: %s\nMax error:\
+        %f" % (varKey, np.max(diffMat))
+
+
+@pytest.mark.analytic
+def test_cyldifn(testDir, tol):
+    _test_analytic(testDir + "/testAnalytCylDifn", tol,
+                   (defs.testAnalytCylDifn, defs.analytCylDifn))
+
+
+@pytest.mark.analytic
+def test_sphdifn(testDir, tol):
+    _test_analytic(testDir + "/testAnalytSphDifn", tol,
+                   (defs.testAnalytSphDifn, defs.analytSphDifn))
+
+
+def _test_analytic(testDir, tol, info):
+    newDir = osp.join(testDir, "sim_output")
+    newDatah5 = False
+    newDataFile = osp.join(newDir, "output_data.mat")
+    if not osp.exists(newDataFile):
+        newDataFile = osp.join(newDir, "output_data.hdf5")
+        newDatah5 = True
+    assert osp.exists(newDataFile), "neither output_data.{mat,hdf5} present"
+
+    dD_s, ndD_s = IO.read_dicts(osp.join(newDir, "input_dict_system"))
+    tmp = IO.read_dicts(osp.join(newDir, "input_dict_c"))
+    dD_e = {}
+    ndD_e = {}
+    dD_e["c"], ndD_e["c"] = tmp
+    try:
+        if newDatah5:
+            newData = h5py.File(newDataFile, 'r')
+        else:
+            newData = sio.loadmat(newDataFile)
+    except IOError as exception:
+        # If it's an error _other than_ the file not being there
+        assert exception.errno == errno.ENOENT, "IO error on opening file"
+        assert False, "File %s does not exist" % (newDataFile)
+    t_ref = dD_s["t_ref"]
+    L_part = dD_s["psd_len"]["c"][0, 0]
+    nx_part = ndD_s["psd_num"]["c"][0, 0]
+    t_refPart = L_part ** 2 / dD_e["c"]["D"]
+    # Skip first time point: analytical solution fails at t=0.
+    t0ind = 2
+    r0ind = 1
+    if newDatah5:
+        tvecA = newData["phi_applied_times"][t0ind:] * (t_ref / t_refPart)
+    else:
+        tvecA = newData["phi_applied_times"][0][t0ind:] * (t_ref / t_refPart)
+    cmat = newData["partTrodecvol0part0_c"]
+    cmin, cmax = np.min(cmat), np.max(cmat)
+    delC = cmax - cmin
+    # Skip center mesh point and first time points:
+    # analytical solutions have issues with r=0 and t=0
+    cmat = cmat[t0ind:, r0ind:]
+    nt = len(tvecA)
+    # Skip center mesh point: analytical solution as sin(r)/r
+    xvecA = np.linspace(0., 1., nx_part)[1:]
+    R, T = np.meshgrid(xvecA, tvecA)
+    theta = info[1](R, T)
+    theta = delC * theta + cmin
+    for tind in range(nt):
+        cvec = cmat[tind, :]
+        thetavec = theta[tind, :]
+        assert np.max(np.abs(thetavec - cvec)) < tol, "Fail from tolerance"
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 00000000..114e581f
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,82 @@
+import pytest
+import os.path as osp
+from os import walk
+
+
+def pytest_addoption(parser):
+    parser.addoption(
+        "--modDir",
+        action="store",
+        default=None,
+        help="The directory with the modified test results"
+    )
+    parser.addoption(
+        "--baseDir",
+        action="store",
+        default=None,
+        help="The directory with the base test results"
+    )
+    parser.addoption(
+        "--tolerance",
+        action="store",
+        default=1e-4,
+        help="tolerance for test cases"
+    )
+    parser.addoption(
+        "--slowdown_tolerance",
+        action="store",
+        default=2,
+        help="tolerance for slowdown of test cases"
+    )
+    parser.addoption(
+        "--tests",
+        action="store",
+        default=[],
+        nargs='*',
+        help="test cases to check"
+    )
+    parser.addoption(
+        "--skip-analytic", action="store_true", default=False, help="skip analytic tests"
+    )
+
+
+def pytest_configure(config):
+    config.addinivalue_line("markers", "analytic: mark test as analytic")
+
+
+def pytest_collection_modifyitems(config, items):
+    if config.getoption("--skip-analytic"):
+        skip_analytic = pytest.mark.skip(reason="--skip-analytic defined")
+        for item in items:
+            if "analytic" in item.keywords:
+                item.add_marker(skip_analytic)
+
+
+def pytest_generate_tests(metafunc):
+    if "Dirs" in metafunc.fixturenames:
+        dir_t = metafunc.config.getoption("modDir")
+        dir_b = metafunc.config.getoption("baseDir")
+        if metafunc.config.getoption("tests") == []:
+            _, directories, _ = next(walk(dir_b))
+            directories.sort()
+            try:
+                directories.remove('plots')  # Directories to ignore
+            except Exception:
+                pass
+
+            metafunc.parametrize("Dirs", [(osp.join(dir_b, dir),
+                                           osp.join(dir_t, dir)) for dir in directories])
+        else:
+            # Sometimes the tests option is a list, other times it is
+            # an array with one element (a string of tests).
+            # This handles both.
+            tests_str = " ".join(metafunc.config.getoption("tests"))
+            tests_lst = tests_str.split()
+            metafunc.parametrize("Dirs", [(dir_b + "/" + test, dir_t + "/" + test)
+                                          for test in tests_lst])
+    if "tol" in metafunc.fixturenames:
+        metafunc.parametrize("tol", [float(metafunc.config.getoption("tolerance"))])
+    if "testDir" in metafunc.fixturenames:
+        dir_t = metafunc.config.getoption("modDir")
+        metafunc.parametrize("testDir", [dir_t])
+    pytest.slowdown_tolerance = float(metafunc.config.getoption("slowdown_tolerance"))
diff --git a/tests/ref_outputs/20200430_092037 b/tests/ref_outputs/20200430_092037
deleted file mode 100644
index e69de29b..00000000
diff --git a/tests/ref_outputs/test001/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test001/sim_output/input_dict_system_dD.p
index 8a6f7397..194ce46a 100644
Binary files a/tests/ref_outputs/test001/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test001/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test001/sim_output/output_data.mat b/tests/ref_outputs/test001/sim_output/output_data.mat
index e3ae110b..b437c231 100644
Binary files a/tests/ref_outputs/test001/sim_output/output_data.mat and b/tests/ref_outputs/test001/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test001/sim_output/run_info.txt b/tests/ref_outputs/test001/sim_output/run_info.txt
index a62a0316..2f576bdd 100644
--- a/tests/ref_outputs/test001/sim_output/run_info.txt
+++ b/tests/ref_outputs/test001/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.1266016960144043 s
+Total run time: 0.14665699005126953 s
diff --git a/tests/ref_outputs/test002/sim_output/daetools_config_options.txt b/tests/ref_outputs/test002/sim_output/daetools_config_options.txt
index 89c9d5d9..5b8dc25f 100644
--- a/tests/ref_outputs/test002/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test002/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test002/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test002/sim_output/input_dict_system_dD.p
index 5b151d1c..53d0af49 100644
Binary files a/tests/ref_outputs/test002/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test002/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test002/sim_output/output_data.mat b/tests/ref_outputs/test002/sim_output/output_data.mat
index 183c6c41..9e353ac1 100644
Binary files a/tests/ref_outputs/test002/sim_output/output_data.mat and b/tests/ref_outputs/test002/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test002/sim_output/run_info.txt b/tests/ref_outputs/test002/sim_output/run_info.txt
index d755ffb3..1605fad5 100644
--- a/tests/ref_outputs/test002/sim_output/run_info.txt
+++ b/tests/ref_outputs/test002/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.13753271102905273 s
+Total run time: 0.11190462112426758 s
diff --git a/tests/ref_outputs/test003/sim_output/daetools_config_options.txt b/tests/ref_outputs/test003/sim_output/daetools_config_options.txt
index 89c9d5d9..5b8dc25f 100644
--- a/tests/ref_outputs/test003/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test003/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test003/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test003/sim_output/input_dict_system_dD.p
index 1312ae80..da75a085 100644
Binary files a/tests/ref_outputs/test003/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test003/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test003/sim_output/output_data.mat b/tests/ref_outputs/test003/sim_output/output_data.mat
index 1885609b..05bee3a1 100644
Binary files a/tests/ref_outputs/test003/sim_output/output_data.mat and b/tests/ref_outputs/test003/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test003/sim_output/run_info.txt b/tests/ref_outputs/test003/sim_output/run_info.txt
index ff01881b..42bab6e2 100644
--- a/tests/ref_outputs/test003/sim_output/run_info.txt
+++ b/tests/ref_outputs/test003/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.14322519302368164 s
+Total run time: 0.11441421508789062 s
diff --git a/tests/ref_outputs/test004/sim_output/daetools_config_options.txt b/tests/ref_outputs/test004/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test004/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test004/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test004/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test004/sim_output/input_dict_system_dD.p
index 1312ae80..da75a085 100644
Binary files a/tests/ref_outputs/test004/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test004/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test004/sim_output/output_data.mat b/tests/ref_outputs/test004/sim_output/output_data.mat
index 0d0ffe91..ad8b1ef9 100644
Binary files a/tests/ref_outputs/test004/sim_output/output_data.mat and b/tests/ref_outputs/test004/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test004/sim_output/run_info.txt b/tests/ref_outputs/test004/sim_output/run_info.txt
index 83c07555..1bca0126 100644
--- a/tests/ref_outputs/test004/sim_output/run_info.txt
+++ b/tests/ref_outputs/test004/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 2.2543089389801025 s
+Total run time: 2.438649892807007 s
diff --git a/tests/ref_outputs/test005/sim_output/daetools_config_options.txt b/tests/ref_outputs/test005/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test005/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test005/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test005/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test005/sim_output/input_dict_system_dD.p
index 4c63e08b..2d62ba7f 100644
Binary files a/tests/ref_outputs/test005/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test005/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test005/sim_output/output_data.mat b/tests/ref_outputs/test005/sim_output/output_data.mat
index 09233382..8300569d 100644
Binary files a/tests/ref_outputs/test005/sim_output/output_data.mat and b/tests/ref_outputs/test005/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test005/sim_output/run_info.txt b/tests/ref_outputs/test005/sim_output/run_info.txt
index 39ad2055..ef7ea8b7 100644
--- a/tests/ref_outputs/test005/sim_output/run_info.txt
+++ b/tests/ref_outputs/test005/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.09043455123901367 s
+Total run time: 0.06692743301391602 s
diff --git a/tests/ref_outputs/test006/sim_output/daetools_config_options.txt b/tests/ref_outputs/test006/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test006/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test006/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test006/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test006/sim_output/input_dict_system_dD.p
index a8f24381..e1517373 100644
Binary files a/tests/ref_outputs/test006/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test006/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test006/sim_output/output_data.mat b/tests/ref_outputs/test006/sim_output/output_data.mat
index 26a311ef..830fe964 100644
Binary files a/tests/ref_outputs/test006/sim_output/output_data.mat and b/tests/ref_outputs/test006/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test006/sim_output/run_info.txt b/tests/ref_outputs/test006/sim_output/run_info.txt
index 09e3890a..3ff1934f 100644
--- a/tests/ref_outputs/test006/sim_output/run_info.txt
+++ b/tests/ref_outputs/test006/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.05935525894165039 s
+Total run time: 0.04276466369628906 s
diff --git a/tests/ref_outputs/test007/sim_output/daetools_config_options.txt b/tests/ref_outputs/test007/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test007/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test007/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test007/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test007/sim_output/input_dict_system_dD.p
index 974063ac..3a9bcc90 100644
Binary files a/tests/ref_outputs/test007/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test007/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test007/sim_output/output_data.mat b/tests/ref_outputs/test007/sim_output/output_data.mat
index c5edfdcd..9595498b 100644
Binary files a/tests/ref_outputs/test007/sim_output/output_data.mat and b/tests/ref_outputs/test007/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test007/sim_output/run_info.txt b/tests/ref_outputs/test007/sim_output/run_info.txt
index 6c2e731d..5c57ad82 100644
--- a/tests/ref_outputs/test007/sim_output/run_info.txt
+++ b/tests/ref_outputs/test007/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.05792498588562012 s
+Total run time: 0.05393362045288086 s
diff --git a/tests/ref_outputs/test008/sim_output/daetools_config_options.txt b/tests/ref_outputs/test008/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test008/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test008/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test008/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test008/sim_output/input_dict_system_dD.p
index 502f5b59..df280a0c 100644
Binary files a/tests/ref_outputs/test008/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test008/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test008/sim_output/output_data.mat b/tests/ref_outputs/test008/sim_output/output_data.mat
index 3f1d6c75..228939de 100644
Binary files a/tests/ref_outputs/test008/sim_output/output_data.mat and b/tests/ref_outputs/test008/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test008/sim_output/run_info.txt b/tests/ref_outputs/test008/sim_output/run_info.txt
index 4a376fc6..512a7d91 100644
--- a/tests/ref_outputs/test008/sim_output/run_info.txt
+++ b/tests/ref_outputs/test008/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.2602696418762207 s
+Total run time: 0.26702141761779785 s
diff --git a/tests/ref_outputs/test009/sim_output/daetools_config_options.txt b/tests/ref_outputs/test009/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test009/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test009/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test009/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test009/sim_output/input_dict_system_dD.p
index 4aa63c32..8342e07f 100644
Binary files a/tests/ref_outputs/test009/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test009/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test009/sim_output/output_data.mat b/tests/ref_outputs/test009/sim_output/output_data.mat
index 53716ae9..bf33ba19 100644
Binary files a/tests/ref_outputs/test009/sim_output/output_data.mat and b/tests/ref_outputs/test009/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test009/sim_output/run_info.txt b/tests/ref_outputs/test009/sim_output/run_info.txt
index b9de4f38..698e5f72 100644
--- a/tests/ref_outputs/test009/sim_output/run_info.txt
+++ b/tests/ref_outputs/test009/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.3343966007232666 s
+Total run time: 0.34839534759521484 s
diff --git a/tests/ref_outputs/test010/params_system.cfg b/tests/ref_outputs/test010/params_system.cfg
index d1058292..93c143d5 100644
--- a/tests/ref_outputs/test010/params_system.cfg
+++ b/tests/ref_outputs/test010/params_system.cfg
@@ -14,7 +14,7 @@ prevDir = false
 tend = 1.2e3
 tsteps = 35
 relTol = 1e-6
-absTol = 1e-6
+absTol = 1e-7
 T = 298
 randomSeed = true
 seed = 10
diff --git a/tests/ref_outputs/test010/sim_output/daetools_config_options.txt b/tests/ref_outputs/test010/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test010/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test010/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test010/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test010/sim_output/input_dict_system_dD.p
index 5b151d1c..53d0af49 100644
Binary files a/tests/ref_outputs/test010/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test010/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test010/sim_output/input_dict_system_ndD.p b/tests/ref_outputs/test010/sim_output/input_dict_system_ndD.p
index f7d01d4b..a83055cc 100644
Binary files a/tests/ref_outputs/test010/sim_output/input_dict_system_ndD.p and b/tests/ref_outputs/test010/sim_output/input_dict_system_ndD.p differ
diff --git a/tests/ref_outputs/test010/sim_output/input_params_system.cfg b/tests/ref_outputs/test010/sim_output/input_params_system.cfg
index d1058292..93c143d5 100644
--- a/tests/ref_outputs/test010/sim_output/input_params_system.cfg
+++ b/tests/ref_outputs/test010/sim_output/input_params_system.cfg
@@ -14,7 +14,7 @@ prevDir = false
 tend = 1.2e3
 tsteps = 35
 relTol = 1e-6
-absTol = 1e-6
+absTol = 1e-7
 T = 298
 randomSeed = true
 seed = 10
diff --git a/tests/ref_outputs/test010/sim_output/output_data.mat b/tests/ref_outputs/test010/sim_output/output_data.mat
index 4517fbbd..491a427b 100644
Binary files a/tests/ref_outputs/test010/sim_output/output_data.mat and b/tests/ref_outputs/test010/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test010/sim_output/run_info.txt b/tests/ref_outputs/test010/sim_output/run_info.txt
index fc3e1abd..db5c70e7 100644
--- a/tests/ref_outputs/test010/sim_output/run_info.txt
+++ b/tests/ref_outputs/test010/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.3
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-0ed6bd4
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 1.634394645690918 s
+Total run time: 1.1364917755126953 s
diff --git a/tests/ref_outputs/test011/sim_output/daetools_config_options.txt b/tests/ref_outputs/test011/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test011/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test011/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test011/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test011/sim_output/input_dict_system_dD.p
index 1312ae80..da75a085 100644
Binary files a/tests/ref_outputs/test011/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test011/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test011/sim_output/output_data.mat b/tests/ref_outputs/test011/sim_output/output_data.mat
index f7c57e23..ddbb848c 100644
Binary files a/tests/ref_outputs/test011/sim_output/output_data.mat and b/tests/ref_outputs/test011/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test011/sim_output/run_info.txt b/tests/ref_outputs/test011/sim_output/run_info.txt
index 4a9c62ee..b9a6a643 100644
--- a/tests/ref_outputs/test011/sim_output/run_info.txt
+++ b/tests/ref_outputs/test011/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 1.1732561588287354 s
+Total run time: 1.518448829650879 s
diff --git a/tests/ref_outputs/test012/sim_output/daetools_config_options.txt b/tests/ref_outputs/test012/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test012/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test012/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test012/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test012/sim_output/input_dict_system_dD.p
index 2e20a9f7..3202ba53 100644
Binary files a/tests/ref_outputs/test012/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test012/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test012/sim_output/output_data.mat b/tests/ref_outputs/test012/sim_output/output_data.mat
index 4d2e9075..b183f30d 100644
Binary files a/tests/ref_outputs/test012/sim_output/output_data.mat and b/tests/ref_outputs/test012/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test012/sim_output/run_info.txt b/tests/ref_outputs/test012/sim_output/run_info.txt
index e606d084..7770d4b6 100644
--- a/tests/ref_outputs/test012/sim_output/run_info.txt
+++ b/tests/ref_outputs/test012/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.16162872314453125 s
+Total run time: 0.1365950107574463 s
diff --git a/tests/ref_outputs/test013/sim_output/daetools_config_options.txt b/tests/ref_outputs/test013/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test013/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test013/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test013/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test013/sim_output/input_dict_system_dD.p
index 84fa6392..b316ae40 100644
Binary files a/tests/ref_outputs/test013/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test013/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test013/sim_output/output_data.mat b/tests/ref_outputs/test013/sim_output/output_data.mat
index 5b031287..96399fe6 100644
Binary files a/tests/ref_outputs/test013/sim_output/output_data.mat and b/tests/ref_outputs/test013/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test013/sim_output/run_info.txt b/tests/ref_outputs/test013/sim_output/run_info.txt
index 515e3a01..fdecf0d2 100644
--- a/tests/ref_outputs/test013/sim_output/run_info.txt
+++ b/tests/ref_outputs/test013/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.18322324752807617 s
+Total run time: 0.16135025024414062 s
diff --git a/tests/ref_outputs/test014/sim_output/daetools_config_options.txt b/tests/ref_outputs/test014/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test014/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test014/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test014/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test014/sim_output/input_dict_system_dD.p
index 83a509cd..7696b197 100644
Binary files a/tests/ref_outputs/test014/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test014/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test014/sim_output/output_data.mat b/tests/ref_outputs/test014/sim_output/output_data.mat
index 5e1a9a64..92405100 100644
Binary files a/tests/ref_outputs/test014/sim_output/output_data.mat and b/tests/ref_outputs/test014/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test014/sim_output/run_info.txt b/tests/ref_outputs/test014/sim_output/run_info.txt
index e17bfc89..ca69183c 100644
--- a/tests/ref_outputs/test014/sim_output/run_info.txt
+++ b/tests/ref_outputs/test014/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.10962629318237305 s
+Total run time: 0.10097002983093262 s
diff --git a/tests/ref_outputs/test015/sim_output/daetools_config_options.txt b/tests/ref_outputs/test015/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test015/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test015/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test015/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test015/sim_output/input_dict_system_dD.p
index 9788508c..d8c4d0d2 100644
Binary files a/tests/ref_outputs/test015/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test015/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test015/sim_output/output_data.mat b/tests/ref_outputs/test015/sim_output/output_data.mat
index d2d1432a..839d4580 100644
Binary files a/tests/ref_outputs/test015/sim_output/output_data.mat and b/tests/ref_outputs/test015/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test015/sim_output/run_info.txt b/tests/ref_outputs/test015/sim_output/run_info.txt
index c81ae050..ee03e116 100644
--- a/tests/ref_outputs/test015/sim_output/run_info.txt
+++ b/tests/ref_outputs/test015/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.47542333602905273 s
+Total run time: 0.4992647171020508 s
diff --git a/tests/ref_outputs/test016/params_system.cfg b/tests/ref_outputs/test016/params_system.cfg
index cdbd98ed..39f87aa9 100644
--- a/tests/ref_outputs/test016/params_system.cfg
+++ b/tests/ref_outputs/test016/params_system.cfg
@@ -10,7 +10,7 @@ segments = [
 	(-0.5, 0.1),
 	]
 tramp = 1e-3
-prevDir = /home/cogswell/Projects/mpet/tests/test_outputs/20200430_092037/test016/../test008/sim_output
+prevDir = ../test008/sim_output
 tend = 1.2e3
 tsteps = 35
 relTol = 1e-6
diff --git a/tests/ref_outputs/test016/sim_output/daetools_config_options.txt b/tests/ref_outputs/test016/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test016/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test016/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test016/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test016/sim_output/input_dict_system_dD.p
index 07aeac28..4a594543 100644
Binary files a/tests/ref_outputs/test016/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test016/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test016/sim_output/input_dict_system_ndD.p b/tests/ref_outputs/test016/sim_output/input_dict_system_ndD.p
index 549388e4..129d74c2 100644
Binary files a/tests/ref_outputs/test016/sim_output/input_dict_system_ndD.p and b/tests/ref_outputs/test016/sim_output/input_dict_system_ndD.p differ
diff --git a/tests/ref_outputs/test016/sim_output/input_params_system.cfg b/tests/ref_outputs/test016/sim_output/input_params_system.cfg
index cdbd98ed..39f87aa9 100644
--- a/tests/ref_outputs/test016/sim_output/input_params_system.cfg
+++ b/tests/ref_outputs/test016/sim_output/input_params_system.cfg
@@ -10,7 +10,7 @@ segments = [
 	(-0.5, 0.1),
 	]
 tramp = 1e-3
-prevDir = /home/cogswell/Projects/mpet/tests/test_outputs/20200430_092037/test016/../test008/sim_output
+prevDir = ../test008/sim_output
 tend = 1.2e3
 tsteps = 35
 relTol = 1e-6
diff --git a/tests/ref_outputs/test016/sim_output/output_data.mat b/tests/ref_outputs/test016/sim_output/output_data.mat
index ce071351..f9df7de2 100644
Binary files a/tests/ref_outputs/test016/sim_output/output_data.mat and b/tests/ref_outputs/test016/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test016/sim_output/run_info.txt b/tests/ref_outputs/test016/sim_output/run_info.txt
index 52cfe2c7..d54aba50 100644
--- a/tests/ref_outputs/test016/sim_output/run_info.txt
+++ b/tests/ref_outputs/test016/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.1737227439880371 s
+Total run time: 0.15644240379333496 s
diff --git a/tests/ref_outputs/test017/params_system.cfg b/tests/ref_outputs/test017/params_system.cfg
index 57d035b2..4f864190 100644
--- a/tests/ref_outputs/test017/params_system.cfg
+++ b/tests/ref_outputs/test017/params_system.cfg
@@ -10,7 +10,7 @@ segments = [
 	(-0.5, 0.1),
 	]
 tramp = 1e-3
-prevDir = /home/cogswell/Projects/mpet/tests/test_outputs/20200430_092037/test017/../test008/sim_output
+prevDir = ../test008/sim_output
 tend = 3e3
 tsteps = 35
 relTol = 1e-6
diff --git a/tests/ref_outputs/test017/sim_output/daetools_config_options.txt b/tests/ref_outputs/test017/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test017/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test017/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test017/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test017/sim_output/input_dict_system_dD.p
index b9481559..4cfe34e8 100644
Binary files a/tests/ref_outputs/test017/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test017/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test017/sim_output/input_dict_system_ndD.p b/tests/ref_outputs/test017/sim_output/input_dict_system_ndD.p
index cdd71460..a2de2101 100644
Binary files a/tests/ref_outputs/test017/sim_output/input_dict_system_ndD.p and b/tests/ref_outputs/test017/sim_output/input_dict_system_ndD.p differ
diff --git a/tests/ref_outputs/test017/sim_output/input_params_system.cfg b/tests/ref_outputs/test017/sim_output/input_params_system.cfg
index 60efaf74..4f864190 100644
--- a/tests/ref_outputs/test017/sim_output/input_params_system.cfg
+++ b/tests/ref_outputs/test017/sim_output/input_params_system.cfg
@@ -10,7 +10,7 @@ segments = [
 	(-0.5, 0.1),
 	]
 tramp = 1e-3
-prevDir = /home/cogswell/Projects/mpet/tests/test_outputs/20200702_105612/test017/../test008/sim_output
+prevDir = ../test008/sim_output
 tend = 3e3
 tsteps = 35
 relTol = 1e-6
diff --git a/tests/ref_outputs/test017/sim_output/output_data.mat b/tests/ref_outputs/test017/sim_output/output_data.mat
index 33e5fe4a..d85d412d 100644
Binary files a/tests/ref_outputs/test017/sim_output/output_data.mat and b/tests/ref_outputs/test017/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test017/sim_output/run_info.txt b/tests/ref_outputs/test017/sim_output/run_info.txt
index f4572a28..8ed14a4a 100644
--- a/tests/ref_outputs/test017/sim_output/run_info.txt
+++ b/tests/ref_outputs/test017/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.3
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-0ed6bd4
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.15526533126831055 s
+Total run time: 0.16385197639465332 s
diff --git a/tests/ref_outputs/test018/sim_output/daetools_config_options.txt b/tests/ref_outputs/test018/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test018/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test018/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test018/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test018/sim_output/input_dict_system_dD.p
index 2c945416..1f4ea4d2 100644
Binary files a/tests/ref_outputs/test018/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test018/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test018/sim_output/output_data.mat b/tests/ref_outputs/test018/sim_output/output_data.mat
index 97a461a3..f12f3df2 100644
Binary files a/tests/ref_outputs/test018/sim_output/output_data.mat and b/tests/ref_outputs/test018/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test018/sim_output/run_info.txt b/tests/ref_outputs/test018/sim_output/run_info.txt
index 848bab70..e38fa82d 100644
--- a/tests/ref_outputs/test018/sim_output/run_info.txt
+++ b/tests/ref_outputs/test018/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.08303952217102051 s
+Total run time: 0.05601668357849121 s
diff --git a/tests/ref_outputs/test019/sim_output/daetools_config_options.txt b/tests/ref_outputs/test019/sim_output/daetools_config_options.txt
index d2bbfe95..c07d6540 100644
--- a/tests/ref_outputs/test019/sim_output/daetools_config_options.txt
+++ b/tests/ref_outputs/test019/sim_output/daetools_config_options.txt
@@ -31,7 +31,7 @@
         },
         "activity": {
             "printHeader": "false",
-            "printStats": "true",
+            "printStats": "false",
             "timeHorizon": "100.0",
             "reportingInterval": "1.0",
             "reportTimeDerivatives": "false",
diff --git a/tests/ref_outputs/test019/sim_output/input_dict_system_dD.p b/tests/ref_outputs/test019/sim_output/input_dict_system_dD.p
index 370935b0..a22c6056 100644
Binary files a/tests/ref_outputs/test019/sim_output/input_dict_system_dD.p and b/tests/ref_outputs/test019/sim_output/input_dict_system_dD.p differ
diff --git a/tests/ref_outputs/test019/sim_output/output_data.mat b/tests/ref_outputs/test019/sim_output/output_data.mat
index b722081c..a1f58260 100644
Binary files a/tests/ref_outputs/test019/sim_output/output_data.mat and b/tests/ref_outputs/test019/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test019/sim_output/run_info.txt b/tests/ref_outputs/test019/sim_output/run_info.txt
index c6bd83e2..8e99a1c2 100644
--- a/tests/ref_outputs/test019/sim_output/run_info.txt
+++ b/tests/ref_outputs/test019/sim_output/run_info.txt
@@ -1,11 +1,11 @@
 mpet version:
-0.1.2
+0.1.5
 
 branch name:
-development
+fix/finite_volumes
 
 commit hash:
-96c732a
+85d2a26
 
 to run, from the root repo directory, copy relevant files there,
 edit input_params_system.cfg to point to correct material
@@ -14,4 +14,4 @@ $ git checkout [commit hash]
 $ patch -p1 < commit.diff:
 $ python[3] mpetrun.py input_params_system.cfg
 
-Total run time: 0.2591094970703125 s
+Total run time: 0.21903109550476074 s
diff --git a/tests/test_defs.py b/tests/test_defs.py
index 41257140..8775a877 100644
--- a/tests/test_defs.py
+++ b/tests/test_defs.py
@@ -7,6 +7,7 @@
 import tests.compare_plots as cmpr
 import mpet.main as main
 import mpet.io_utils as IO
+import re
 
 
 def corePlots(testDir, dirDict):
@@ -27,517 +28,12 @@ def electrodePlots(testDir, dirDict, trode):
     cmpr.cbarLine(testDir, dirDict, "c")
 
 
-def test001(testDir, dirDict, pflag):
-    """ LFP ACR C3 """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "ACR")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test002(testDir, dirDict, pflag):
-    """ LFP CHR cylinder """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "CHR")
-    P.set("Particles", "shape", "cylinder")
-    P.set("Material", "muRfunc", "LiFePO4")
-    P.set("Material", "Dfunc", "lattice")
-    P.set("Material", "dgammadc", "5e-30")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test003(testDir, dirDict, pflag):
-    """ LFP CHR sphere """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "CHR")
-    P.set("Particles", "shape", "sphere")
-    P.set("Material", "muRfunc", "LiFePO4")
-    P.set("Material", "Dfunc", "lattice")
-    P.set("Material", "dgammadc", "-2e-30")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test004(testDir, dirDict, pflag):
-    """ LFP CHR sphere with noise  """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "CHR")
-    P.set("Particles", "discretization", "2e-9")
-    P.set("Particles", "shape", "sphere")
-    P.set("Material", "muRfunc", "LiFePO4")
-    P.set("Material", "noise", "true")
-    P.set("Material", "numnoise", "100")
-    P.set("Material", "Dfunc", "lattice")
-    P.set("Material", "dgammadc", "-2e-30")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test005(testDir, dirDict, pflag):
-    """ LFP homog """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    P_s.set("Sim Params", "Nvol_c", "2")
-    P_s.set("Sim Params", "Nvol_s", "2")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-        cmpr.bulkpf(testDir, dirDict, "c")
-
-
-def test006(testDir, dirDict, pflag):
-    """ LFP homog with logPad, Vmin """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    P_s.set("Sim Params", "Vmin", "2.8e-0")
-    P_s.set("Particles", "cs0_c", "0.5")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test007(testDir, dirDict, pflag):
-    """ LFP homog with logPad, Vmax """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "-1e-2")
-    P_s.set("Sim Params", "Vmax", "4.0e-0")
-    P_s.set("Particles", "cs0_c", "0.5")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test008(testDir, dirDict, pflag):
-    """ LFP homog_sdn """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    P_s.set("Sim Params", "Nvol_c", "3")
-    P_s.set("Sim Params", "Npart_c", "3")
-    P_s.set("Particles", "stddev_c", "25e-9")
-    P_s.set("Sim Params", "capFrac", "0.6")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog_sdn")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-        cmpr.bulkpf(testDir, dirDict, "c")
-
-
-def test009(testDir, dirDict, pflag):
-    """ Graphite-2param homog """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog2")
-    P.set("Material", "muRfunc", "LiC6")
-    P.set("Material", "noise", "true")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test010(testDir, dirDict, pflag):
-    """ Graphite-2param CHR cylinder """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "CHR2")
-    P.set("Particles", "discretization", "2e-9")
-    P.set("Particles", "shape", "cylinder")
-    P.set("Material", "muRfunc", "LiC6")
-    P.set("Material", "Dfunc", "lattice")
-    P.set("Material", "B", "0.1e9")
-    P.set("Material", "kappa", "2.5e-10")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test011(testDir, dirDict, pflag):
-    """ Graphite-2param CHR sphere """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    P_s.set("Sim Params", "relTol", "1e-7")
-    P_s.set("Sim Params", "absTol", "1e-7")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "CHR2")
-    P.set("Particles", "discretization", "2e-9")
-    P.set("Particles", "shape", "sphere")
-    P.set("Material", "muRfunc", "LiC6")
-    P.set("Material", "Dfunc", "lattice")
-    P.set("Material", "B", "0.1e9")
-    P.set("Material", "kappa", "2.5e-10")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test012(testDir, dirDict, pflag):
-    """ Solid solution, diffn sphere, homog, LiC6_coke_ss2, LiMn2O4_ss2
-    BV_mod01, BV_mod02
-    cathode + anode
-    """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_a.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrodec = osp.join(testDir, "params_c.cfg")
-    ptrodea = osp.join(testDir, "params_a.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    P_s.set("Sim Params", "capFrac", "0.67")
-    P_s.set("Sim Params", "tsteps", "25")
-    P_s.set("Sim Params", "Nvol_c", "2")
-    P_s.set("Sim Params", "Nvol_a", "2")
-    P_s.set("Sim Params", "Vmin", "2e0")
-    P_s.set("Particles", "cs0_c", "0.2")
-    P_s.set("Particles", "cs0_a", "0.495")
-    P_s.set("Electrolyte", "elyteModelType", "SM")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrodec)
-    P.set("Particles", "type", "homog")
-    P.set("Particles", "shape", "cylinder")
-    P.set("Material", "muRfunc", "LiMn2O4_ss2")
-    P.set("Reactions", "rxnType", "BV_mod01")
-    IO.write_config_file(P, ptrodec)
-    P = IO.get_config(ptrodea)
-    P.set("Particles", "discretization", "2.5e-9")
-    P.set("Particles", "shape", "sphere")
-    P.set("Material", "muRfunc", "LiC6_coke_ss2")
-    P.set("Reactions", "rxnType", "BV_mod02")
-    IO.write_config_file(P, ptrodea)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-        cmpr.bulkpf(testDir, dirDict, "c")
-        electrodePlots(testDir, dirDict, "a")
-        cmpr.bulkpf(testDir, dirDict, "a")
-
-
-def test013(testDir, dirDict, pflag):
-    """ Solid solution, diffn cylinder, homog, testIS_ss, LiMn2O4_ss2
-    Marcus, BV_raw
-    cathode + separator + anode
-    """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_a.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrodec = osp.join(testDir, "params_c.cfg")
-    ptrodea = osp.join(testDir, "params_a.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "profileType", "CV")
-    P_s.set("Sim Params", "Vset", "3.8")
-    P_s.set("Sim Params", "tsteps", "25")
-    P_s.set("Sim Params", "Nvol_c", "2")
-    P_s.set("Sim Params", "Nvol_s", "2")
-    P_s.set("Sim Params", "Nvol_a", "2")
-    P_s.set("Particles", "cs0_c", "0.2")
-    P_s.set("Particles", "cs0_a", "0.95")
-    P_s.set("Electrolyte", "elyteModelType", "dilute")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrodec)
-    P.set("Particles", "type", "homog")
-    P.set("Particles", "shape", "sphere")
-    P.set("Material", "muRfunc", "LiMn2O4_ss2")
-    P.set("Reactions", "rxnType", "Marcus")
-    IO.write_config_file(P, ptrodec)
-    P = IO.get_config(ptrodea)
-    P.set("Particles", "discretization", "2.5e-9")
-    P.set("Particles", "shape", "cylinder")
-    P.set("Material", "muRfunc", "testIS_ss")
-    P.set("Reactions", "rxnType", "BV_raw")
-    IO.write_config_file(P, ptrodea)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-        cmpr.bulkpf(testDir, dirDict, "c")
-        electrodePlots(testDir, dirDict, "a")
-        cmpr.bulkpf(testDir, dirDict, "a")
-
-
-def test014(testDir, dirDict, pflag):
-    """ LFP homog with CCsegments, MHC, Rser """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "profileType", "CCsegments")
-    P_s.set("Sim Params", "segments",
-            "[(1., 25), (-2., 10), (0., 30)]")
-    P_s.set("Sim Params", "tramp", "1e+0")
-    P_s.set("Sim Params", "Rser", "1e-3")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog")
-    P.set("Material", "muRfunc", "LiFePO4")
-    P.set("Reactions", "rxnType", "MHC")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test015(testDir, dirDict, pflag):
-    """ testRS homog with CVsegments, bulkCond, partCond """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "profileType", "CVsegments")
-    P_s.set("Sim Params", "Npart_c", "3")
-    P_s.set("Sim Params", "Nvol_c", "3")
-    P_s.set("Sim Params", "segments",
-            "[(-0.3, 25), (0., 10), (0.3, 30)]")
-    P_s.set("Sim Params", "tramp", "1e+0")
-    P_s.set("Sim Params", "absTol", "1e-5")
-    P_s.set("Conductivity", "simBulkCond_c", "true")
-    P_s.set("Conductivity", "sigma_s_c", "1e-2")
-    P_s.set("Conductivity", "simPartCond_c", "true")
-    P_s.set("Conductivity", "G_mean_c", "1e-14")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog")
-    P.set("Material", "muRfunc", "testRS")
-    P.set("Reactions", "rxnType", "BV_raw")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-        cmpr.bulkpf(testDir, dirDict, "c")
-
-
-def test016(testDir, dirDict, pflag):
-    """ test CC continuation """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Crate", "1e-2")
-    P_s.set("Sim Params", "Nvol_c", "3")
-    P_s.set("Sim Params", "Npart_c", "3")
-    P_s.set("Sim Params", "capFrac", "0.34")
-    test008dir = str(osp.join(testDir, "..", "test008", "sim_output"))
-    P_s.set("Sim Params", "prevDir", test008dir)
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog_sdn")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-        cmpr.bulkpf(testDir, dirDict, "c")
-
-
-def test017(testDir, dirDict, pflag):
-    """ test CV continuation """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "profileType", "CV")
-    P_s.set("Sim Params", "Vset", "3.45")
-    P_s.set("Sim Params", "tend", "3e3")
-    P_s.set("Sim Params", "Nvol_c", "3")
-    P_s.set("Sim Params", "Npart_c", "3")
-    test008dir = str(osp.join(testDir, "..", "test008", "sim_output"))
-    P_s.set("Sim Params", "prevDir", test008dir)
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog_sdn")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-        cmpr.bulkpf(testDir, dirDict, "c")
-
-
-def test018(testDir, dirDict, pflag):
-    """ Like test014, LFP homog with CCsegments, BV, Rfilm, Rfilm_foil """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "profileType", "CCsegments")
-    P_s.set("Sim Params", "segments",
-            "[(1., 25), (-2., 10), (0., 30)]")
-    P_s.set("Sim Params", "tramp", "1e+0")
-    P_s.set("Electrodes", "k0_foil", "3e+0")
-    P_s.set("Electrodes", "Rfilm_foil", "3e+0")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog")
-    P.set("Material", "muRfunc", "LiFePO4")
-    P.set("Reactions", "Rfilm", "1e+1")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
-
-
-def test019(testDir, dirDict, pflag):
-    """ SM electrolyte with LFP homog, sep + cathode, significant elyte polarization """
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_system.cfg"), testDir)
-    shutil.copy(osp.join(dirDict["baseConfig"], "params_c.cfg"), testDir)
-    psys = osp.join(testDir, "params_system.cfg")
-    ptrode = osp.join(testDir, "params_c.cfg")
-    P_s = IO.get_config(psys)
-    P_s.set("Sim Params", "Nvol_c", "3")
-    P_s.set("Sim Params", "Nvol_s", "3")
-    P_s.set("Geometry", "L_c", "120e-6")
-    P_s.set("Geometry", "L_s", "90e-6")
-    P_s.set("Electrolyte", "elyteModelType", "SM")
-    IO.write_config_file(P_s, psys)
-    P = IO.get_config(ptrode)
-    P.set("Particles", "type", "homog")
-    P.set("Material", "muRfunc", "LiFePO4")
-    IO.write_config_file(P, ptrode)
-    main.main(psys, keepArchive=False)
-    shutil.move(dirDict["simOut"], testDir)
-    if pflag:
-        corePlots(testDir, dirDict)
-        elytePlots(testDir, dirDict)
-        electrodePlots(testDir, dirDict, "c")
+def getNumberOfTests():
+    n = 0
+    for key in globals().keys():
+        if re.match(r'test[0-9][0-9][0-9]$', key):
+            n = n + 1
+    return n
 
 
 def testAnalytSphDifn(testDir, dirDict):
@@ -603,12 +99,12 @@ def testAnalytCylDifn(testDir, dirDict):
 def analytSphDifn(R, T):
     """ Analytical solution for diffusion in a sphere """
     n = 30
-    lmbdavec = np.pi*np.arange(1, n+1)
-    theta = 0*R
+    lmbdavec = np.pi * np.arange(1, n + 1)
+    theta = 0 * R
     for i, lmbda in enumerate(lmbdavec):
-        theta += ((2./lmbda**2) * np.sin(lmbda*R)/R
-                  * (np.sin(lmbda) - lmbda*np.cos(lmbda))
-                  * np.exp(-lmbda**2*T))
+        theta += ((2. / lmbda ** 2) * np.sin(lmbda * R) / R
+                  * (np.sin(lmbda) - lmbda * np.cos(lmbda))
+                  * np.exp(-lmbda ** 2 * T))
     return theta
 
 
@@ -616,8 +112,8 @@ def analytCylDifn(R, T):
     """ Analytical solution for diffusion in a cylinder """
     n = 30
     lmbdavec = spcl.jn_zeros(0, n)
-    theta = 0*R
+    theta = 0 * R
     for i, lmbda in enumerate(lmbdavec):
-        theta += ((2./lmbda) * spcl.j0(lmbda*R)/spcl.j1(lmbda)
-                  * np.exp(-lmbda**2*T))
+        theta += ((2. / lmbda) * spcl.j0(lmbda * R) / spcl.j1(lmbda)
+                  * np.exp(-lmbda ** 2 * T))
     return theta
diff --git a/tests/test_suite.py b/tests/test_suite.py
index b72f088b..a6e1e974 100644
--- a/tests/test_suite.py
+++ b/tests/test_suite.py
@@ -1,22 +1,36 @@
 import errno
 import os
 import os.path as osp
+from os import walk
+import shutil
 import time
 
-import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 import scipy.io as sio
 
+import mpet.main
 import mpet.io_utils as IO
 import tests.test_defs as defs
 
 
 def run_test_sims(runInfo, dirDict, pflag=True):
-    for testStr in sorted(runInfo.keys()):
+    for testStr in runInfo:
         testDir = osp.join(dirDict["out"], testStr)
         os.makedirs(testDir)
-        runInfo[testStr](testDir, dirDict, pflag)
+
+        # Copy config files from ref_outputs
+        refDir = osp.join(dirDict["suite"],"ref_outputs",testStr)
+        _, _, filenames = next(walk(osp.join(refDir)))
+        for f in filenames:
+            if ".cfg" in f:
+                shutil.copyfile(osp.join(refDir,f),osp.join(testDir,f))
+
+        # Run the simulation
+        configfile = osp.join(testDir,'params_system.cfg')
+        mpet.main.main(configfile, keepArchive=False)
+        shutil.move(dirDict["simOut"], testDir)
+
     # Remove the history directory that mpet creates.
     try:
         os.rmdir(osp.join(dirDict["suite"], "history"))
@@ -113,7 +127,7 @@ def compare_with_ref(runInfo, dirDict, tol=1e-4):
                 raise
             print("No simulation data for " + testStr)
             continue
-        
+
         refData = sio.loadmat(refDataFile)
         for varKey in (set(refData.keys()) & set(newData.keys())):
             # TODO -- Consider keeping a list of the variables that fail
@@ -122,7 +136,7 @@ def compare_with_ref(runInfo, dirDict, tol=1e-4):
             if varKey[0:2] == "__":
                 continue
 
-            #Compute the difference between the solution and the reference
+            # Compute the difference between the solution and the reference
             try:
                 varDataNew = newData[varKey]
                 varDataRef = refData[varKey]
@@ -135,14 +149,14 @@ def compare_with_ref(runInfo, dirDict, tol=1e-4):
                 print(testStr, "Fail from KeyError")
                 testFailed = True
                 continue
-            
+
             # #Check absolute and relative error against tol
             if np.max(diffMat) > tol and np.max(diffMat) > tol*np.max(np.abs(varDataRef)):
                 print(testStr, "Fail from tolerance")
                 print("variable failing:", varKey)
                 print("max error:", np.max(diffMat))
                 testFailed = True
-        
+
         if testFailed:
             failList.append(testStr)
 
@@ -197,7 +211,7 @@ def main(compareDir):
     failList = compare_with_ref(runInfo, dirDict, tol=1e-4)
     failListAnalyt = compare_with_analyt(runInfoAnalyt, dirDict, tol=1e-4)
 
-    #Print a summary of test results
+    # Print a summary of test results
     print("\n")
     print("--------------------------------------------------------------")
     if len(failList) > 0:
diff --git a/tests/timings_tests.py b/tests/timings_tests.py
new file mode 100644
index 00000000..6dfc6123
--- /dev/null
+++ b/tests/timings_tests.py
@@ -0,0 +1,17 @@
+import os.path as osp
+import pytest
+
+
+def get_sim_time(simDir):
+    with open(osp.join(simDir, "run_info.txt")) as fi:
+        simTime = float(fi.readlines()[-1].split()[-2])
+    return simTime
+
+
+def test_compare_timings(Dirs, tol):
+    refDir, testDir = Dirs
+    newDir = osp.join(testDir, "sim_output")
+    refDir = osp.join(refDir, "sim_output")
+    time_new = get_sim_time(newDir)
+    time_ref = get_sim_time(refDir)
+    assert time_new < time_ref * pytest.slowdown_tolerance, "too high slowdown with reference"