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

Issue:
Move new solid branch documentation from README to docs



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index c14538b0..ee6559fa 100644
--- a/.gitignore
+++ b/.gitignore
@@ -44,3 +44,8 @@ venv/
 # documentation
 docs/_build
 docs/apidocs
+
+# ignore python biuld files
+build/
+# ignore daetools
+daetools*
\ No newline at end of file
diff --git a/bin/create_ensemble.py b/bin/create_ensemble.py
new file mode 100755
index 00000000..0abdfd8b
--- /dev/null
+++ b/bin/create_ensemble.py
@@ -0,0 +1,52 @@
+#!/usr/bin/env python3
+
+import configparser
+import sys
+import itertools
+import os
+
+
+def ensemble_definitions():
+    # Values that need ensemble
+    ensemble = [
+        [("Geometry","L_c"), ["2","43"]],
+        [("Geometry","L_a"), ["3","5","43"]],
+    ]
+
+    # helpers
+    keys = [vals[0] for vals in ensemble]
+    val = [vals[1] for vals in ensemble]
+    return keys, val
+
+
+def create_ensemble(cff, keys=None, val=None):
+    with open('ensemble_parallel_configs.txt', "w") as ff:
+        if keys is None and val is None:
+            keys, val = ensemble_definitions()
+        cfg = configparser.ConfigParser()
+        cfg.optionxform = str
+        cfg.read(cff)
+        # Create all variations
+        combinations = list(itertools.product(*val))
+        for combination in combinations:
+            params = dict(zip(keys,combination))
+            new_cfg = cfg
+            nicename = []
+            for key, val in params.items():
+                new_cfg[key[0]][key[1]] = val
+                nicename.append(key[1] + "=" + val)
+
+            # Write config
+            cfg_dir = os.path.dirname(cff)
+            with open(cfg_dir + "/" + "-".join(nicename) + ".cfg", "w") as f:
+                new_cfg.write(f)
+                ff.write(str(cfg_dir + "/" + "-".join(nicename) + ".cfg\n"))
+    return
+
+
+if __name__ == '__main__':
+    # Read in file
+    if len(sys.argv) < 2:
+        print("need the config file [python create_enamble.py <baseconfig>]")
+        exit(1)
+    create_ensemble(sys.argv[1])
diff --git a/bin/mpet_create_runjobs_dashboard.py b/bin/mpet_create_runjobs_dashboard.py
new file mode 100644
index 00000000..692f4238
--- /dev/null
+++ b/bin/mpet_create_runjobs_dashboard.py
@@ -0,0 +1,101 @@
+from create_ensemble import create_ensemble
+from run_jobs import create_slurm_cluster, create_pbs_cluster, create_local_cluster, run_mpet
+# import mpet_plot_app
+import os
+import subprocess
+import shutil
+from dask.distributed import Client
+
+
+# --------------------------------------------------------------
+# Fill all these out
+# --------------------------------------------------------------
+''' Ensemble settings/ '''
+# The default config file that you want to adjust in the ensemble
+# (path relative to main folder)
+cff = 'configs/params_system.cfg'
+# Values that need ensemble
+ensemble = [
+    [("Sim Params","Nvol_c"), ["4", "5"]]
+]
+
+
+class ClusterSettings():
+    '''Cluster settings.'''
+    # time = 00:00:00  # Max walltime per job (hh:mm:ss format). Argument is not used with a local.
+    nproc = 1  # type=int, Number of CPU cores per job. Argument is not used with a local cluster.
+    mem = 1  # Max memory usage per job. For alocal cluster it sets the memory limit per worker.
+    queue = 1  # Queue to use. Argument is not used with a local cluster.
+    dashboard_port = 4096  # Port for dask dashboard
+
+
+class MainSettings():
+    '''Cluster script settings.'''
+    scheduler = 'local'  # choices=('slurm', 'pbs', 'local'); Scheduling system to use
+    min_jobs = 1  # int; Minimum number of jobs to launch. Argument is not used with local cluster.
+    max_jobs = 1  # int; Maximum number of jobs to launch. Argument is not used with local cluster.
+    # Text file containg the path to each MPET config file to run;
+    # 'parallel_configs.txt' is what create_ensemble saves to automatically (so don't change it)
+    mpet_configs = 'ensemble_parallel_configs.txt'
+
+
+# --------------------------------------------------------------
+'''
+End of parameters to fill out.
+To run: execute run_mpet_create_run_dashboard.sh from the terminal;
+. ./run_mpet_create_run_dashboard.sh
+To see the resulting dashboard open the http link presented in the terminal:
+" Dash is running on http://127.0.0.1:8050/ "
+'''
+# --------------------------------------------------------------
+
+
+# helpers for ensemble, do not change
+keys = [vals[0] for vals in ensemble]
+val = [vals[1] for vals in ensemble]
+
+
+def call_run_cluster(output_folder):
+    # split cluster settings from the rest
+    args = ClusterSettings()
+    main_settings = MainSettings()
+    '''
+    for arg in ["scheduler", "mpet_configs", "min_jobs", "max_jobs"]:
+        main_settings[arg] = getattr(args, arg)
+        delattr(args.__class__, arg)'''
+    cluster_settings = vars(args)
+
+    # create cluster
+    if main_settings.scheduler == 'slurm':
+        cluster = create_slurm_cluster(**cluster_settings)
+    elif main_settings.scheduler == 'pbs':
+        cluster = create_pbs_cluster(**cluster_settings)
+    elif main_settings.scheduler == 'local':
+        cluster = create_local_cluster(args.mem, args.dashboard_port)
+
+    # Scale Dask cluster automatically based on scheduler activity (only if not local cluster)
+    if main_settings.scheduler != 'local':
+        cluster.adapt(minimum_jobs=main_settings['min_jobs'],
+                      maximum_jobs=main_settings['max_jobs'])
+    client = Client(cluster)
+
+    run_mpet(client, output_folder, os.path.abspath(main_settings.mpet_configs))
+    return
+
+
+if __name__ == '__main__':
+    # Read in config file
+    create_ensemble(cff, keys, val)
+
+    # Define output folder
+    # Store output in folder this script was called from
+    output_folder = './runjobs_dashboard'
+    # remove sim output if it already exists to only keep newest output
+    if os.path.exists(os.path.join(output_folder, 'sim_output')):
+        shutil.rmtree(os.path.join(output_folder, 'sim_output'))
+    # create output folder if it does not exist yet
+    if not os.path.exists(output_folder):
+        os.mkdir(output_folder)
+    call_run_cluster(output_folder)
+    subprocess.call(["python", "./bin/mpet_plot_app.py", "-d",
+                     str(os.path.join(output_folder, 'sim_output'))])
diff --git a/bin/mpet_plot_app.py b/bin/mpet_plot_app.py
new file mode 100644
index 00000000..94c90ee3
--- /dev/null
+++ b/bin/mpet_plot_app.py
@@ -0,0 +1,586 @@
+import dash
+from dash import html, dcc
+import dash_bootstrap_components as dbc
+from dash.dependencies import Input, Output
+import plotly.express as px
+
+import numpy as np
+import plotly.graph_objects as go
+from plotly.subplots import make_subplots
+
+
+from mpet.config import constants
+
+from mpet.plot import plot_data_db
+
+dff, dff_c_sub, dff_cd_sub, dff_bulkp, dff_csld_sub, df_cbar = plot_data_db.main()
+
+###################################################################
+# Part 2
+# Define plots (define that there is a grpah and perform callback (see part 3))
+# style elements
+###################################################################
+app = dash.Dash(__name__, external_stylesheets=[dbc.themes.CERULEAN])
+
+# Define colors to be used
+colors = {
+    'background': '#FFFFFF',
+    'lithium': '#2fa4e7',
+    'text': '#000000',
+    'bg_table': '#B9B6B6',
+    'text_table': '#111111',
+}
+
+defaultmodel = (dff['Model'].unique()[0])
+# Define components of app
+app.layout = html.Div([
+    dbc.Row(dbc.Col(
+        html.H1("MPET visualisation",
+            style={'textAlign': 'center', 'background': colors['lithium'], 'color': 'white',
+                   'margin-top': '10px', 'margin-bottom': '10px'}),
+        ),),
+    dbc.Row(dbc.Col([
+        html.Div("Select models to display in all plots",
+            style={'color': colors['lithium'], 'margin-left': '20px'}),
+        dcc.Checklist(
+            options=dff['Model'].unique(),
+            value=dff['Model'].unique(),
+            id='model-selection',
+            labelStyle={'display': 'block'},
+            style={'margin-bottom': '20px', 'margin-left': '20px',},
+            )]
+        ),),
+    html.Hr(style={"color": colors['lithium'], "height":'5px'}),
+    html.H2("General properties", style={'margin-left': '20px'}),
+    html.Hr(style={"color": colors['lithium'], "height":'5px'}),
+    dbc.Row([
+        dbc.Col([
+            # v/vt
+            dbc.Card([
+                html.H4(
+                    children='Voltage',
+                    style={'textAlign': 'center', 'font-family':'Sans-serif', 'margin-top': '10px'}
+                    ),
+                dcc.Dropdown(
+                    ['Time (s)', 'Cathode Filling Fraction'],
+                    'Time (s)',
+                    id='xaxis-column',
+                    style={'font-family':'Sans-serif', 'margin-left': '10px', 'width': '80%'}),
+                dcc.Loading(children=[dcc.Graph(id='Voltage-graph-double')],
+                            color=colors['lithium'], type="dot", fullscreen=False)],
+                style={'margin-bottom': '20px', 'margin-left': '10px'}
+                ),
+            # power
+            dbc.Card([
+                html.H4(children='Power', style={
+                    'textAlign': 'center', 'font-family':'Sans-serif', 'margin-top': '10px'}),
+                dcc.Loading(children=[dcc.Graph(id="power")],
+                            color=colors['lithium'], type="dot", fullscreen=False,)
+                ], style={'margin-bottom': '20px', 'margin-left': '10px'})
+            ],
+            width=4,),
+        dbc.Col([
+            # Curr
+            dbc.Card([
+                html.H4(
+                    children='Current profile',
+                    style={'textAlign': 'center', 'font-family':'Sans-serif',
+                           'margin-bottom': '45px', 'margin-top': '10px'}),
+                dcc.Loading(children=[dcc.Graph(id="current")],
+                            color=colors['lithium'], type="dot", fullscreen=False),
+            ], style={'margin-bottom': '20px'}),
+            # elytecons
+            dbc.Card([
+                html.H4(children='Average concentration of electrolyte',
+                        style={'textAlign': 'center', 'font-family':'Sans-serif',
+                               'margin-top': '10px'}),
+                dcc.Loading(children=[dcc.Graph(id='elytecons'),],
+                            color=colors['lithium'], type="dot", fullscreen=False,),
+            ], style={'margin-bottom': '20px'})],
+            width=4,),
+        dbc.Col([
+            # soc
+            dbc.Card([
+                html.H4(children='Overall utilization',
+                        style={'textAlign': 'center', 'font-family':'Sans-serif',
+                               'margin-top': '10px'}),
+                html.H4(children='(state of charge of electrode)',
+                        style={'textAlign': 'center', 'font-family':'Sans-serif'}),
+                dcc.Loading(children=[
+                    html.H5(
+                        children='Cathode',
+                        style={'textAlign': 'center', 'font-family':'Sans-serif',
+                            'color': '#7f7f7f', 'margin-top': '10px'}),
+                    dcc.Graph(id='Cathode-filling-fraction'),
+                    html.H5(
+                        children='Anode',
+                        style={'textAlign': 'center', 'font-family':'Sans-serif',
+                               'color': '#7f7f7f', 'margin-top': '5px', 'margin-bottom': '0px'}),
+                    dcc.Graph(id='Anode-filling-fraction'),
+                    ],
+                    color=colors['lithium'], type="dot", fullscreen=False,),
+                ], style={'margin-bottom': '20px', 'margin-right': '10px'})],
+                width=4,),
+        ]
+    ),
+    html.Hr(style={"color": colors['lithium'], "height":'5px'}),
+    html.H2("Details of electrolyte", style={'margin-left': '20px'}),
+    html.Hr(style={"color": colors['lithium'], "height":'5px'}),
+    dbc.Row([
+        dbc.Col([
+            # elytec
+            dbc.Card([
+                html.H4(children='Electrolyte concentration',
+                        style={'textAlign': 'center', 'margin-top': '10px'}),
+                dcc.Loading(children=[dcc.Graph(id='electrolyte-concentration-ani')],
+                            color=colors['lithium'], type="dot", fullscreen=False,),
+            ], style={'margin-bottom': '20px', 'margin-left': '10px'}),
+            # elytep
+            dbc.Card([
+                html.H4(children='Electrolyte potential',
+                        style={'textAlign': 'center', 'margin-top': '10px'}),
+                dcc.Loading(children=[dcc.Graph(id='electrolyte-potential-ani'),],
+                            color=colors['lithium'], type="dot", fullscreen=False,),
+            ], style={'margin-bottom': '20px', 'margin-left': '10px'}),
+        ], width=6),
+        dbc.Col([
+            # elytei
+            dbc.Card([
+                html.H4(children='Electrolyte current density',
+                        style={'textAlign': 'center', 'margin-top': '10px'}),
+                dcc.Loading(children=[dcc.Graph(id='electrolyte-cd-ani'),],
+                            color=colors['lithium'], type="dot", fullscreen=False,),
+            ], style={'margin-bottom': '20px', 'margin-right': '10px'}),
+            # elytedivi
+            dbc.Card([
+                html.H4(children='Divergence of electrolyte current density',
+                        style={'textAlign': 'center', 'margin-top': '10px'}),
+                dcc.Loading(children=[dcc.Graph(id='electrolyte-decd-ani'),],
+                            color=colors['lithium'], type="dot", fullscreen=False,),
+            ], style={'margin-bottom': '20px', 'margin-right': '10px'}),
+        ], width=6),
+    ]),
+    # bulkp
+    dbc.Row(dbc.Col(dbc.Card([
+        html.H4(children='Macroscopic cathode/anode solid phase potential',
+                style={'textAlign': 'center', 'margin-top': '10px'}),
+        dcc.Loading(children=[dcc.Graph(id='bulkp'),],
+                    color=colors['lithium'], type="dot", fullscreen=False,),
+        ], style={'margin-left': '10px', 'margin-right': '10px'}), width=12)),
+    #
+    html.Hr(style={"color": colors['lithium'], "height":'5px'}),
+    html.H2("Details of concentrations for individual models", style={'margin-left': '20px'}),
+    html.Hr(style={"color": colors['lithium'], "height":'5px'}),
+    dbc.Row(dbc.Col([
+        html.Div("Select model to display",
+            style={'color': colors['lithium'], 'margin-left': '20px'}),
+        dcc.Dropdown(options=dff['Model'].unique(), value=defaultmodel,
+                     id='select_single_model',
+                     style={'width':'50%', 'margin-left': '10px', 'margin-bottom': '20px'})]
+        ),),
+    # cbar
+    dbc.Row(dbc.Col(dbc.Card([
+        html.H4(children='Average solid concentrations',
+                style={'textAlign': 'center', 'font-family':'Sans-serif', 'margin-top': '10px'}),
+        dcc.Loading(children=[dcc.Graph(id='cbar_c'),
+                              dcc.Graph(id='cbar_c2'),
+                              dcc.Graph(id='cbar_a'),
+                              dcc.Graph(id='cbar_a2'),],
+                    color=colors['lithium'], type="dot", fullscreen=False,),
+    ], style={'margin-left': '10px', 'margin-right': '10px', 'margin-bottom': '20px'}))),
+    # surf
+    dbc.Row(dbc.Col(dbc.Card([
+        html.H4(children='Solid surface concentration',
+                style={'textAlign': 'center', 'font-family':'Sans-serif', 'margin-top': '10px'}),
+        dcc.Loading(children=[dcc.Graph(id='Surface-concentration-cathode'),
+                              dcc.Graph(id='Surface-concentration-anode'),],
+                    color=colors['lithium'], type="dot", fullscreen=False,),
+    ], style={'margin-left': '10px', 'margin-right': '10px', 'margin-bottom': '20px',
+              "overflow": "scroll"}))),
+    # csld
+    dbc.Row(dbc.Col(dbc.Card([
+        html.H4(children='All solid concentrations',
+                style={'textAlign': 'center', 'font-family':'Sans-serif', 'margin-top': '10px'}),
+        html.H5(children='Time percentage',
+                style={'textAlign': 'left', 'font-family':'Sans-serif', 'margin-left': '10px'}),
+        dcc.Slider(
+            0,100,step=1,value=0,
+            id='timefraction_slider',
+            marks={
+                0: '0', 5: '5', 10: '10', 15: '15', 20: '20',
+                25: '25', 30: '30', 35: '35', 40: '40', 45: '45',
+                50: '50', 55: '55', 60: '60', 65: '65', 70: '70',
+                75: '75', 80: '80', 85: '85', 90: '90', 95: '95',
+                100: '100'
+            }, tooltip={"placement": "top", "always_visible": True}),
+        dcc.Graph(id='csld_c'),
+        dcc.Graph(id='csld_a'),
+    ], style={'margin-left': '10px', 'margin-right': '10px', 'margin-bottom': '20px',
+              "overflow": "scroll"}))),
+    # cbarline
+    dbc.Row(dbc.Col(dbc.Card([
+        html.H4(children='Average concentration in each particle of electrode',
+                style={'textAlign': 'center', 'font-family':'Sans-serif', 'margin-top': '10px'}),
+        dcc.Loading(children=[dcc.Graph(id='cbarline_c'),
+                              dcc.Graph(id='cbarline_a'),],
+                    color=colors['lithium'], type="dot", fullscreen=False,),
+    ], style={'margin-left': '10px', 'margin-right': '10px', 'margin-bottom': '20px',
+              "overflow": "scroll"}))),
+])
+
+
+###################################################################
+# Part 3
+# Define callback functions to get interactive plots
+###################################################################
+@app.callback(
+    Output('current', 'figure'),
+    Output('power', 'figure'),
+    Output('Cathode-filling-fraction', 'figure'),
+    Output('Anode-filling-fraction', 'figure'),
+    Output('elytecons', 'figure'),
+    Input('model-selection', 'value')
+    )
+def update_graphs_multimodels(model_selection):
+    m_select = dff[np.in1d(dff['Model'], model_selection)]
+    # plots
+    current = plot_multimodel(m_select, 'Current', ytitle='Current (C-rate)')
+    power = plot_multimodel(m_select, 'Power', ytitle=u'Power (W/m\u00b2)')
+    fillfrac_c = plot_multimodel(m_select, 'Cathode Filling Fraction',
+                                 ytitle='Cathode Filling Fraction')
+    fillfrac_a = plot_multimodel(m_select, 'Anode Filling Fraction',
+                                 ytitle='Anode Filling Fraction')
+    elytecons = plot_multimodel(m_select, 'cavg',
+                                ytitle='Avg. Concentration of electrolyte (nondim)')
+    return current, power, fillfrac_c, fillfrac_a, elytecons
+
+
+@app.callback(
+    Output('Voltage-graph-double', 'figure'),
+    Input('xaxis-column', 'value'),
+    Input('model-selection', 'value')
+    )
+def update_graphs_multimodels_voltage(xaxis_column_name, model_selection):
+    m_select = dff[np.in1d(dff['Model'], model_selection)]
+    voltage = plot_multimodel(m_select, 'Voltage (V)', xaxis=xaxis_column_name)
+    return voltage
+
+
+@app.callback(
+    Output('electrolyte-concentration-ani', 'figure'),
+    Output('electrolyte-potential-ani', 'figure'),
+    Output('electrolyte-cd-ani', 'figure'),
+    Output('electrolyte-decd-ani', 'figure'),
+    Output('bulkp', 'figure'),
+    Input('model-selection', 'value')
+    )
+def callback_multimodel_movies(model_selection):
+    m_select_dff_c = dff_c_sub[np.in1d(dff_c_sub['Model'], model_selection)]
+    m_select_dff_cd = dff_cd_sub[np.in1d(dff_cd_sub['Model'], model_selection)]
+    m_select_dff_bulkp = dff_bulkp[np.in1d(dff_bulkp['Model'], model_selection)]
+    electrolyte_concentration = ani_elytrolyte(m_select_dff_c, "cellsvec",
+                                               "Concentration electrolyte",
+                                               "Time fraction", 'Concentration of electrolyte (M)')
+    eletrolyte_potential = ani_elytrolyte(m_select_dff_c, "cellsvec", "Potential electrolyte",
+                                          "Time fraction", 'Potential of electrolyte (V)')
+    curr_dens = ani_elytrolyte(m_select_dff_cd, "facesvec",
+                               "Curreny density electrolyte",
+                               "Time fraction", u"Current density of electrolyte (A/m\u00b2)")
+    div_curr_dens = ani_elytrolyte(m_select_dff_c, "cellsvec",
+                                   "Divergence electrolyte curr dens",
+                                   "Time fraction",
+                                   u"Divergence electrolyte current density (A/m\u00b3)")
+    bulkp = ani_bulkp(m_select_dff_bulkp)
+    return electrolyte_concentration, eletrolyte_potential, curr_dens, div_curr_dens, bulkp
+
+
+@app.callback(
+    Output('Surface-concentration-cathode', 'figure'),
+    Output('Surface-concentration-anode', 'figure'),
+    Output('cbarline_c', 'figure'),
+    Output('cbarline_a', 'figure'),
+    Output('cbar_c', 'figure'),
+    Output('cbar_c2', 'figure'),
+    Output('cbar_a', 'figure'),
+    Output('cbar_a2', 'figure'),
+    Input('select_single_model', 'value')
+    )
+def update_graphs_single_models(select_single_model):
+    m_select = dff[dff['Model'] == select_single_model]
+    if m_select["Config trode type"].iloc[0] in constants.one_var_types:
+        trode_type = 1
+    elif m_select["Config trode type"].iloc[0] in constants.two_var_types:
+        trode_type = 2
+    # plots
+    subplt_solid_surf_con_c = subplt_solid_surf_con("Cathode", m_select, trode_type)
+    subplt_solid_surf_con_a = subplt_solid_surf_con("Anode", m_select, trode_type)
+    cbarline_c, cbarline_a = subplots_cbarlinec(m_select, trode_type)
+    cbar_c, cbar_c2, cbar_a, cbar_a2 = ani_cbar(select_single_model, trode_type)
+    return (subplt_solid_surf_con_c, subplt_solid_surf_con_a,
+            cbarline_c, cbarline_a, cbar_c, cbar_c2, cbar_a, cbar_a2
+            )
+
+
+@app.callback(
+    Output('csld_c', 'figure'),
+    Output('csld_a', 'figure'),
+    Input('select_single_model', 'value'),
+    Input('timefraction_slider', 'value')
+    )
+def update_csld(select_single_model, timefraction_slider):
+    m_select_csld = dff_csld_sub[dff_csld_sub['Model'] == select_single_model]
+    csld_c, csld_a = subplots_t_csld(m_select_csld, timefraction_slider)
+    return (csld_c, csld_a)
+
+
+# Plot functions
+def plot_multimodel(df, yaxis, ytitle=None, xaxis='Time (s)'):
+    plot = px.line(df, x=xaxis, y=yaxis, color="Model")
+    if xaxis != 'Time (s)':
+        plot.update_xaxes(title=xaxis)
+    if ytitle:
+        plot.update_yaxes(title=ytitle)
+    plot.update_layout(legend=dict(
+        yanchor="bottom",
+        y=1.02,
+        xanchor="left",
+        x=0.01
+    ))
+    return plot
+
+
+def ani_elytrolyte(df, xname, yname, ani, ytitle):
+    max_y = max(df[yname])
+    min_y = min(df[yname])
+    range_min = 0.9*min_y if min_y > 0 else 1.1*min_y
+    range_max = 1.1*max_y if max_y > 0 else 1.1*max_y
+    if not (max_y == 0 and min_y == 0):
+        fig = px.line(df,
+                      x=xname,
+                      y=yname,
+                      color="Model",
+                      animation_frame=ani,
+                      animation_group=xname)
+    else:
+        fig = px.line(title='Data not availale for selected model(s)')
+    fig.update_yaxes(title=ytitle,
+                     range=[range_min, range_max])
+    fig.update_xaxes(title=u'Battery Position (\u00B5m)')
+    fig.update_layout(legend=dict(
+        yanchor="bottom",
+        y=1.02,
+        xanchor="left",
+        x=0.01
+    ))
+    return fig
+
+
+def subplt_solid_surf_con(trodes, df, trode_type):
+    trode = trodes[0].lower()
+    try:
+        pfx = max(df["pfx"])
+        sStr = max(df["sStr"])
+        r = int(max(df["Npart"+trode]))
+        c = int(max(df["Nvol"+trode]))
+        fig = make_subplots(rows=r, cols=c, shared_xaxes=True, shared_yaxes=True,
+                            x_title='Time (s)', y_title='Solid surface concentration',
+                            row_titles=['Particle ' + str(n) for n in range(1, r+1)],
+                            column_titles=['Volume ' + str(n) for n in range(1, c+1)])
+        for rr in range(0, r):
+            for cc in range(0, c):
+                if trode_type == 2:
+                    str1_base = (pfx
+                                 + "partTrode{trode}vol{{vInd}}part{{pInd}}".format(trode=trode)
+                                 + sStr + "c1")
+                    str2_base = (pfx
+                                 + "partTrode{trode}vol{{vInd}}part{{pInd}}".format(trode=trode)
+                                 + sStr + "c2")
+                    sol1_str = str1_base.format(pInd=rr, vInd=cc)
+                    sol2_str = str2_base.format(pInd=rr, vInd=cc)
+                    datax = df['Time (s)']
+                    datay1 = df[sol1_str]
+                    datay2 = df[sol2_str]
+                    fig.add_trace(
+                        trace=go.Scatter(x=datax, y=datay1, line_color='red', name='c1'),
+                        row=rr+1, col=cc+1)
+                    fig.add_trace(
+                        trace=go.Scatter(x=datax, y=datay2, line_color='blue', name='c2'),
+                        row=rr+1, col=cc+1)
+                else:
+                    str_base = (pfx
+                                + "partTrode{trode}vol{{vInd}}part{{pInd}}".format(trode=trode)
+                                + sStr + "c")
+                    sol_str = str_base.format(pInd=rr, vInd=cc)
+                    datax = df['Time (s)']
+                    datay = df[sol_str]
+                    fig.add_trace(
+                        trace=go.Scatter(x=datax, y=datay, line_color='darkslategray'),
+                        row=rr+1, col=cc+1)
+        fig.update_yaxes(range=[0,1.01])
+        fig.update_layout(height=((r+1)*150), width=((c+1)*150), showlegend=False, title=trodes)
+    except ValueError:
+        fig = px.line(title='Selected model has no '+trodes.lower())
+        fig.update_xaxes(visible=False, showgrid=False)
+        fig.update_yaxes(visible=False, showgrid=False)
+    return fig
+
+
+def subplots_t_csld(df, tf):
+    for trodes in ["Cathode", "Anode"]:
+        try:
+            trode = trodes[0].lower()
+            partStr = "partTrode{trode}vol{vInd}part{pInd}" + max(df["sStr"])
+            r = int(max(df["Npart"+trode]))
+            c = int(max(df["Nvol"+trode]))
+            fig = make_subplots(rows=r, cols=c, shared_xaxes=True, shared_yaxes=True,
+                                x_title=u'Position (\u00B5m)',
+                                y_title='Solid Concentrations of Particles in Electrode',
+                                row_titles=['Particle ' + str(n) for n in range(1, r+1)],
+                                column_titles=['Volume ' + str(n) for n in range(1, c+1)])
+            type2c = False
+            df_select = df[df['Time fraction'] == int(tf)]
+            if df_select["Config trode type"].iloc[0] in constants.one_var_types:
+                type2c = False
+            elif df_select["Config trode type"].iloc[0] in constants.two_var_types:
+                type2c = True
+            for rr in range(0, r):
+                for cc in range(0, c):
+                    lens_str = "lens_{vInd}_{pInd}".format(vInd=cc, pInd=rr)
+                    if type2c is True:
+                        c1str_base = (max(df["pfx"])
+                                      + partStr.format(trode=trode, pInd=rr, vInd=cc) + "c1")
+                        c2str_base = (max(df["pfx"])
+                                      + partStr.format(trode=trode, pInd=rr, vInd=cc) + "c2")
+                        c3str_base = (max(df["pfx"])
+                                      + partStr.format(trode=trode, pInd=rr, vInd=cc) + "cav")
+                        datax = df_select[lens_str]
+                        datay1 = df_select[c1str_base]
+                        datay2 = df_select[c2str_base]
+                        datay3 = df_select[c3str_base]
+                        fig.add_trace(
+                            trace=go.Scatter(x=datax, y=datay1, line_color='red', name='c1'),
+                            row=rr+1, col=cc+1)
+                        fig.add_trace(
+                            trace=go.Scatter(x=datax, y=datay2, line_color='blue', name='c2'),
+                            row=rr+1, col=cc+1)
+                        fig.add_trace(
+                            trace=go.Scatter(x=datax, y=datay3, line_color='grey',
+                                             name='avg. c1 & c2'),
+                            row=rr+1, col=cc+1)
+                    else:
+                        cstr_base = (max(df["pfx"]) + partStr.format(trode=trode, vInd=cc, pInd=rr)
+                                     + "c")
+                        cstr = cstr_base
+                        datax = df_select[lens_str]
+                        datay = df_select[cstr]
+                        fig.add_trace(
+                            trace=go.Scatter(x=datax, y=datay, line_color='darkslategray'),
+                            row=rr+1, col=cc+1)
+            fig.update_yaxes(range=[0,1.01])
+            fig.update_layout(height=((r+1)*150), width=((c+1)*150), showlegend=False,
+                              title=trodes+', time = {time} s'.format(
+                              time=round(df_select["time (s)"].iloc[0])))
+        except ValueError:
+            fig = px.line(title='Selected model has no '+trodes.lower())
+            fig.update_xaxes(visible=False, showgrid=False)
+            fig.update_yaxes(visible=False, showgrid=False)
+        if trode == "c":
+            fig1 = fig
+    return fig1, fig
+
+
+def subplots_cbarlinec(df, trode_type):
+    for trodes in ["Cathode", "Anode"]:
+        trode = trodes[0].lower()
+        try:
+            partStr = ("partTrode{trode}vol{{vInd}}part{{pInd}}".format(trode=trode)
+                       + max(df["sStr"]))
+            r = int(max(df["Npart"+trode]))
+            c = int(max(df["Nvol"+trode]))
+            fig = make_subplots(rows=r, cols=c, shared_xaxes=True, shared_yaxes=True,
+                                x_title='Time (s)', y_title='Particle Average Filling Fraction',
+                                row_titles=['Particle ' + str(n) for n in range(1, r+1)],
+                                column_titles=['Volume ' + str(n) for n in range(1, c+1)])
+            for rr in range(0, r):
+                for cc in range(0, c):
+                    datax = df['Time (s)']
+                    if trode_type == 2:
+                        str1_cbar_base = max(df["pfx"]) + partStr + "c1bar"
+                        str2_cbar_base = max(df["pfx"]) + partStr + "c2bar"
+                        sol1_str = str1_cbar_base.format(pInd=rr, vInd=cc)
+                        sol2_str = str2_cbar_base.format(pInd=rr, vInd=cc)
+                        datay = df[sol1_str]
+                        datay2 = df[sol2_str]
+                        fig.add_trace(
+                            trace=go.Scatter(x=datax, y=datay, line_color='red', name='c1bar'),
+                            row=rr+1, col=cc+1)
+                        fig.add_trace(
+                            trace=go.Scatter(x=datax, y=datay2, line_color='blue', name='c2bar'),
+                            row=rr+1, col=cc+1)
+                    else:
+                        str_cbar_base = max(df["pfx"]) + partStr + "cbar"
+                        sol_str = str_cbar_base.format(pInd=rr, vInd=cc)
+                        datay = df[sol_str]
+                        fig.add_trace(
+                            trace=go.Scatter(x=datax, y=datay, line_color='darkslategray'),
+                            row=rr+1, col=cc+1)
+            fig.update_yaxes(range=[0,1.01])
+            fig.update_layout(height=((r+1)*150), width=((c+1)*150),
+                              showlegend=False, title=trodes)
+        except ValueError:
+            fig = px.line(title='Selected model has no '+trodes.lower())
+            fig.update_xaxes(visible=False, showgrid=False)
+            fig.update_yaxes(visible=False, showgrid=False)
+        if trode == "c":
+            fig1 = fig
+    return fig1, fig
+
+
+def ani_cbar(ms_cbar, trode_type):
+    def plot(cbar):
+        plot = px.scatter(df_cbar_select, x='c', y='r', animation_frame='Time',
+                          animation_group='rc',
+                          color=cbar,
+                          range_color=[0,1],
+                          color_continuous_scale='Turbo',
+                          size='Relative size',
+                          height=None,
+                          width=None
+                          )
+        return plot
+    for trodes in ["Cathode", "Anode"]:
+        trode = trodes[0].lower()
+        df_cbar_select = df_cbar[(df_cbar.Model == ms_cbar) & (df_cbar.Trode == trode)]
+        if df_cbar_select.empty:
+            cbar = px.line(title='Selected model has no '+trodes.lower())
+            cbar2 = px.line()
+        else:
+            if trode_type == 1:
+                cbar = plot('Cbar')
+                cbar2 = px.line(title='Selected model has no Cbar2')
+            elif trode_type == 2:
+                cbar = plot('Cbar1')
+                cbar2 = plot('Cbar2')
+                cbar2.update_layout(title=trodes, transition_duration=100, height=None, width=None)
+            cbar.update_layout(title=trodes, transition_duration=100)
+        cbar.update_xaxes(visible=False, showgrid=False)
+        cbar.update_yaxes(visible=False, showgrid=False)
+        cbar2.update_xaxes(visible=False, showgrid=False)
+        cbar2.update_yaxes(visible=False, showgrid=False)
+        if trode == "c":
+            cbar_c = cbar
+            cbar_c2 = cbar2
+    return cbar_c, cbar_c2, cbar, cbar2
+
+
+def ani_bulkp(df):
+    bulkp = px.line(df, x='Position in electrode', y='Potential (nondim)',
+                    animation_frame='Time fraction (%)', color='Model',
+                    line_dash='Trode', log_y=True)
+    bulkp.update_yaxes(range=[-4, 1.1*np.log10(max(df['Potential (nondim)']))])
+    bulkp.update_xaxes(title='Position in electrode (\u00B5m)')
+    return bulkp
+
+
+if __name__ == '__main__':
+    app.run_server(debug=True)
diff --git a/bin/mpetplot.py b/bin/mpetplot.py
index a010fd37..867b6659 100755
--- a/bin/mpetplot.py
+++ b/bin/mpetplot.py
@@ -13,7 +13,7 @@
 
 # Ordered dictionary of plot types
 plotTypes = OrderedDict([
-    ('v','voltage vs filling fraction'),
+    ('v','voltage vs filling fraction. Default plot type.'),
     ('vt','voltage vs time'),
     ('curr','current vs time'),
     ('power','power vs time'),
@@ -27,6 +27,7 @@
     ('surf_a','anode surface concentrations'),
     ('soc_c','cathode state of charge'),
     ('soc_a','anode state of charge'),
+    ('elytecons', 'Avg. Concentration of electrolyte'),
     ('csld_c','solid concentrations of cathode particles'),
     ('csld_a','solid concentrations of anode particles'),
     ('cbarLine_c','average concentration in each cathode particle'),
@@ -50,48 +51,52 @@
 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')
-
+parser.add_argument('-pt', '--plotType', nargs='*', metavar='plotType', help=plotTypesHelp,
+                    choices=plotTypes.keys(), default=['v'])
+parser.add_argument('-t', '--text', choices=['text'],
+                    help='Optionally just convert output to text')
+parser.add_argument('-s', '--save', choices=['save', 'saveonly'],
+                    help='Optionally save the output')
+parser.add_argument('-c', '--color_changes', choices=['discrete','smooth'], default='discrete',
+                    help='Defines discrete or smooth color changes when plotting average '
+                    + 'solid concentrations. Default: discrete')
+parser.add_argument('-st', '--smooth_type', choices=['GnYlRd_1', 'GnYlRd_2', 'GnYlRd_3'],
+                    help='Set color configuration used from colormaps_custom when using smooth '
+                    + 'color changes. Default: GnYlRd_3',
+                    default='GnYlRd_3')
 parser.add_argument('-v','--version', action='version',
                     version='%(prog)s '+__version__)
 args = parser.parse_args()
 
-
 # Get input file from script parameters
-if len(sys.argv) < 2:
+if args.directory is None:
     raise Exception("Need input data directory name")
-indir = sys.argv[1]
+indir = args.directory
 if not os.path.exists(os.path.join(os.getcwd(), indir)):
     raise Exception("Input file doesn't exist")
 # Optionally just convert output to text
-if len(sys.argv) == 3 and sys.argv[2] == "text":
+if args.text == 'text':
     outmat2txt.main(indir)
     sys.exit()
 # Get plot type from script parameters
-plots = []
-if len(sys.argv) > 2:
-    plots.append(sys.argv[2])
-else:
-    plots.append("v")
+plots = args.plotType
 # Save the plot instead of showing on screen?
 # Get from script parameters
 save_flag = False
 print_flag = True
 data_only = False
 save_only = False
-if len(sys.argv) > 3:
-    if sys.argv[3] in ["save", "saveonly"]:
-        save_flag = True
-        if sys.argv[3] == "saveonly":
-            save_only = True
-            print_flag = False
-    else:
-        for i in range(3, len(sys.argv)):
-            plots.append(sys.argv[i])
+color_changes = args.color_changes
+smooth_type = args.smooth_type
+if args.save in ["save", "saveonly"]:
+    save_flag = True
+    if args.save == "saveonly":
+        save_only = True
+        print_flag = False
 out = []
+
 for plot_type in plots:
     out.append(plot_data.show_data(
-        indir, plot_type, print_flag, save_flag, data_only))
+        indir, plot_type, print_flag, save_flag, data_only, color_changes, smooth_type))
 if not save_only:
     plt.show()
diff --git a/bin/run_jobs.py b/bin/run_jobs.py
new file mode 100644
index 00000000..e909c382
--- /dev/null
+++ b/bin/run_jobs.py
@@ -0,0 +1,125 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+import shutil
+import mpet.main as main
+
+from dask_jobqueue import SLURMCluster, PBSCluster
+from dask.distributed import Client, LocalCluster
+
+
+def create_slurm_cluster(time, nproc, mem, queue, dashboard_port):
+    """Create a SLURM cluster for use with dask"""
+    cluster = SLURMCluster(cores=nproc,
+                           processes=nproc,
+                           memory=mem,
+                           queue=queue,
+                           walltime=time,
+                           scheduler_options={'dashboard_address': dashboard_port})
+    return cluster
+
+
+def create_pbs_cluster(time, nproc, mem, queue, dashboard_port):
+    """Create a PBS cluster for use with dask"""
+    cluster = PBSCluster(cores=nproc,
+                         processes=nproc,
+                         memory=mem,
+                         resource_spec=f'nodes=1:ppn={nproc}',
+                         queue=queue,
+                         walltime=time,
+                         scheduler_options={'dashboard_address': dashboard_port})
+    return cluster
+
+
+def create_local_cluster(mem, dashboard_port):
+    """Create a local cluster for use with dask"""
+    cluster = LocalCluster(memory_limit=mem,
+                           dashboard_address=dashboard_port)
+    return cluster
+
+
+def run_mpet(client, output_folder, mpet_configs):
+    """Run MPET on each config file present in the mpet_configs folder"""
+    # In order to copy or move simulation output to current directory, remove old existing folder
+    tmpDir = os.path.join(os.getcwd(), "sim_output")
+    shutil.rmtree(tmpDir, ignore_errors=True)
+
+    with open(mpet_configs, 'r') as fp:
+        config_files = fp.readlines()
+
+    folder = os.path.dirname(mpet_configs)
+    files = [f'{os.path.join(folder, fname.strip())}' for fname in config_files]
+    print('Running mpet for these config files:', files)
+
+    # function to run MPET in directory given as input to run_mpet
+    def run_mpet_instance(*args, **kwargs):
+        print(f"Running in {output_folder} with args={args}, kwargs={kwargs}")
+        os.chdir(output_folder)
+        try:
+            main.main(*args, **kwargs)
+        except Exception as e:
+            print(f'MPET crashed with error: {e}')
+
+    futures = client.map(run_mpet_instance, files, keepFullRun=True)
+    print('Waiting for MPET to finish')
+    client.gather(futures)
+    client.close()
+    print('Done')
+    return
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description='Runs several instances of MPET on a SLURM or'
+                                     ' PBS cluster, or local cluster')
+    # cluster settings
+    parser.add_argument('--time', '-t', required=True,
+                        help='Maximum walltime per job (hh:mm:ss format)')
+    parser.add_argument('--nproc', '-n', type=int, required=True,
+                        help='Number of CPU cores per job')
+    parser.add_argument('--mem', '-m', required=True,
+                        help=('Max memory usage per job. When using a '
+                              'local cluster it sets the memory limit per worker process.'))
+    parser.add_argument('--queue', '-q', default='default',
+                        help='Queue to use (default: %(default)s)')
+    parser.add_argument('--dashboard_port', '-d', type=int, default=4096,
+                        help='Port for dask dashboard (default: %(default)s)')
+
+    # script settings
+    parser.add_argument('--scheduler', '-s', default='slurm', choices=('slurm', 'pbs', 'local'),
+                        help='Scheduling system to use (default: %(default)s)')
+    parser.add_argument('--min_jobs', type=int, default=1,
+                        help='Minimum number of jobs to launch (default: %(default)s)')
+    parser.add_argument('--max_jobs', type=int, default=1,
+                        help='Maximum number of jobs to launch (default: %(default)s)')
+    parser.add_argument('mpet_configs',
+                        help='Text file containg the path to each MPET config file to run')
+
+    args = parser.parse_args()
+    # split cluster settings from the rest
+    main_settings = {}
+    for arg in ['scheduler', 'mpet_configs', 'min_jobs', 'max_jobs']:
+        main_settings[arg] = getattr(args, arg)
+        delattr(args, arg)
+    cluster_settings = vars(args)
+
+    # create cluster
+    if main_settings['scheduler'] == 'slurm':
+        cluster = create_slurm_cluster(**cluster_settings)
+    elif main_settings['scheduler'] == 'pbs':
+        cluster = create_pbs_cluster(**cluster_settings)
+    elif main_settings['scheduler'] == 'local':
+        cluster = create_local_cluster(args.mem, args.dashboard_port)
+
+    # Scale Dask cluster automatically based on scheduler activity (only if not local cluster)
+    if main_settings['scheduler'] != 'local':
+        cluster.adapt(minimum_jobs=main_settings['min_jobs'],
+                      maximum_jobs=main_settings['max_jobs'])
+    client = Client(cluster)
+
+    # Store output in folder this script was called from
+    output_folder = os.getcwd()
+
+    run_mpet(client, output_folder, os.path.abspath(main_settings['mpet_configs']))
+
+    client.shutdown()
diff --git a/configs/params_LTO.cfg b/configs/params_LTO.cfg
new file mode 100644
index 00000000..faa02524
--- /dev/null
+++ b/configs/params_LTO.cfg
@@ -0,0 +1,29 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = CHR 
+discretization = 5e-10
+shape = sphere
+thickness = 20e-9
+
+[Material] 
+muRfunc = LTO
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+Omega_a = 1.43e-20
+kappa = 1.0e-10
+B = 0.0 
+rho_s = 1.32e28
+D = 4.0e-16
+Dfunc = lattice
+dgammadc = 0e-30
+cwet = 0.98
+
+[Reactions]
+rxnType = BV
+k0 = 3.6
+alpha = 0.5
+lambda = 6.26e-20
+Rfilm = 0e-0
\ No newline at end of file
diff --git a/configs/params_system.cfg b/configs/params_system.cfg
index ef5dca6e..b0c1a3d3 100644
--- a/configs/params_system.cfg
+++ b/configs/params_system.cfg
@@ -13,7 +13,7 @@ Crate = 1
 # 1C_current_density = 12.705
 # Voltage cutoffs, V
 Vmax = 3.6
-Vmin = 2.0
+Vmin = 2.5
 # Battery applied voltage (only used for CV), V
 Vset = 0.12
 # Battery constant power (only used for CP), W/m^2
@@ -99,6 +99,11 @@ mean_c = 100e-9
 stddev_c = 1e-9
 mean_a = 100e-9
 stddev_a = 1e-9
+# contact penalty
+fraction_of_contact = 1
+stand_dev_contact = 0
+# for ACR model 
+localized_losses = false
 # Use a specific set of particle sizes for the distribution
 # If false, a randomly generated PSD is used. Otherwise the input should
 # be a 2D list of particle radii with 'Npart' rows and 'Nvol' columns.
@@ -140,6 +145,10 @@ P_L_a = 0.69
 poros_c = 0.4
 poros_a = 0.4
 poros_s = 1.0
+# Specified volume dependent porosity
+specified_poros_c = false
+specified_poros_a = false
+specified_poros_s = false
 # Bruggeman exponent (tortuosity = porosity^bruggExp)
 BruggExp_c = -0.5
 BruggExp_a = -0.5
@@ -148,6 +157,10 @@ BruggExp_s = -0.5
 [Electrolyte]
 # Initial electrolyte conc., mol/m^3
 c0 = 1000
+# for solid electrolyte select the molarity and the maximum molarity
+# that the sublattice could in host (e.g. c0 = 40000 , cmax = 80000)
+# Maximal concentration of host sites (in solid elyte) [mol m^-3]
+cmax = 80000
 # Cation/anion charge number (e.g. 2, -1 for CaCl_2)
 zp = 1
 zm = -1
@@ -162,16 +175,21 @@ 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: use 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 = SM
+# Solid solution parameter for Solid Electrolyte
+# 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: uses the lattice gas model for solid electrolytes.
 #SMset_filename = filename.py   # optional, to load SMset from custom file instead of props_elyte.py
 SMset = valoen_bernardi
 # Reference electrode (defining the electrolyte potential) information:
@@ -184,3 +202,39 @@ sp = -1
 # e.g. for LiPF6 in EC/DMC, Dp = 2.2e-10, Dm = 2.94e-10
 Dp = 2.2e-10
 Dm = 2.94e-10
+# for example for solid
+# Dp = 5.5e-12
+# Dm = 2.94e-15
+
+
+[Interface]
+# Simulate interface region?
+# Options: true, false (Default: false)
+simInterface_c = false
+simInterface_a = false
+c0_int = 40000
+cmax_i = 80000
+# for numerical dicretization
+Nvol_i = 4
+# interface length
+L_i = 20e-9
+BruggExp_i = -0.5
+poros_i = 1.0
+# 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: use the lattice gas model for solid electrolytes.
+interfaceModelType = solid
+# Options:
+#   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: uses the lattice gas model for solid electrolytes.
+interfaceSMset = solid_elyte
+# diff coefficinets of the interface, can also be liquid
+Dp_i = 5.5e-17
+Dm_i = 2.94e-19
\ No newline at end of file
diff --git a/docs/analysis.rst b/docs/analysis.rst
index 86592c99..a51cbb8e 100644
--- a/docs/analysis.rst
+++ b/docs/analysis.rst
@@ -1,12 +1,17 @@
 Analyze the results with built-in tools
-=============================================
+=======================================
 
 
-Analyze output with mpetplot.py (pass output data directory, then plot-type as arguments)
-
- * e.g., voltage plot: ``mpetplot.py sim_output v``
- * other options (full, c, a indicate full cell, cathode, and anode):
+Analyze output with ``mpetplot.py``. Pass the output data directory, then use the optional plotting arguments. The options for ``mpetplot.py`` are:
+   #. ``-pt`` for plotting types
+   #. ``-t`` for saving output to text format
+   #. ``-s`` for options to save the plot
+   #. ``-c`` for color_map options that are used with plot type ``cbar_{full,c,a}``
+   #. ``-st`` to specify the smooth colormap used with plot type ``cbar_{full,c,a}``
 
+1.  Analyze output with plots using ``mpetplot.py``. Pass output data directory, then use ``-pt [plottype]`` with one (or more) of the plot types listed below. Default is ``v``.
+    - e.g., voltage plot: ``mpetplot.py sim_output -pt v``
+    - other options (``full``, ``c``, ``a`` indicate full cell, cathode, and anode):
    #. ``v`` or ``vt`` -- voltage vs filling fraction or vs time
    #. ``curr`` -- current vs time
    #. ``elytec{f}`` -- electrolyte concentration (movie) or final snapshot with f
@@ -16,13 +21,12 @@ Analyze output with mpetplot.py (pass output data directory, then plot-type as a
    #. ``soc_{c,a}`` -- overall utilization / state of charge of electrode
    #. ``csld_{c,a}`` -- solid concentrations of particles in electrode (movie; used with solidType_{c,a} not homog)
    #. ``cbarLine_{c,a}`` -- average concentration in each particle of electrode
-   #. ``cbar_{full,c,a}`` -- average solid concentrations as changing colors (movie)
    #. ``bulkp_{c,a}`` -- macroscopic electrode solid phase potential (movie)
+   #. ``cbar_{full,c,a}`` -- average solid concentrations as changing colors (movie)
+    - There are two options for the color map type that is used: ``smooth`` or ``discrete``. This can be set with the ``-c`` option, e.g., ``mpetplot.py sim_output -pt cbar_full -c discrete``. The default value is ``discrete``.
+    - When using the ``smooth`` color map option, the colors are selected from colormao_custom.npz, which includes three options (``GnYlRd_1``, ``GnYlRd_2``, and ``GnYlRd_3``) that can be selected with the ``st`` option, e.g., ``mpetplot.py sim_output -pt cbar_full -c discrete -st GnYlRd_1``. The default value is ``GnYlRd_3``.
 
-
-
-
-Alternatively, convert the output to plain text (csv) format using : ``mpetplot.py sim_output text`` (or replace sim_output with any subfolder in the history folder).
+2.  Alternatively, convert the output to plain text (csv) format using : ``mpetplot.py sim_output text`` (or replace sim_output with any subfolder in the history folder).
 Then analyze using whatever tools you prefer. If you want to save output to a movie (or figure), add save as an extra argument to ``mpetplot.py``: ``mpetplot.py sim_output cbar save``.
 
-Movie output requires that you have ``ffmpeg`` or ``mencoder`` (part of ``MPlayer``) installed.
\ No newline at end of file
+Movie output requires that you have ``ffmpeg`` or ``mencoder`` (part of ``MPlayer``) installed.
diff --git a/docs/cluster.rst b/docs/cluster.rst
new file mode 100644
index 00000000..4a429d8d
--- /dev/null
+++ b/docs/cluster.rst
@@ -0,0 +1,23 @@
+Running multiple simulations on a cluster
+=========================================
+
+If you have many simulations you want to run, you can use ``bin/run_jobs.py`` to run them efficiently on a cluster using `Dask <https://dask.org>`_, either locally or on a slurm or PBS cluster. Using the parallel running option requires the following package to be installed: ``dask-jobqueue``.
+
+1. Follow steps 1-4 from :doc:`run_sim` for each of the simulations you want to run. Then create a text file in your working directory containing the system parameter files for your simulations. This text file should contain the file names of each of the system parameter configuration files for which you want to run a simulation. For example, if you have all your parameter files saved in the ``configs`` directory, create: ``configs/parallel_configs.txt``, which contains the lines::
+
+    params_system.cfg
+    params_system_XX.cfg
+    params_system_YY.cfg
+    etc.
+
+2. Run multiple simulations on a cluster using ``run_jobs.py``. The simplest way to run it, is to run the script on the login node. Pass the text file containing the system parameter files (e.g. ``configs/parallel_configs.txt``) and the cluster arguments:
+
+    - ``-s``: scheduler type. Options: ``slurm``, ``pbs``, and ``local``. Default is ``slurm``.
+    - ``-t``: Maximum walltime per job (hh:mm:ss format). Argument is not used with a local cluster.
+    - ``-n``: Number of CPU cores and instances of MPET per job. Argument is not used with a local cluster.
+    - ``-m``: Max memory usage per job (e.g. 2GB). When using a local cluster it sets the memory limit per worker process.
+    - ``-q``: Queue to use. Argument is not used with a local cluster.
+    - ``-d``: Port for Dask dashboard (default 4096).
+    - ``--min_jobs``: Minimum number of jobs to launch. Default = 1. Argument is not used with a local cluster.
+    - ``--max_jobs``: Maximum number of jobs to launch. Default = 1. Argument is not used with a local cluster.
+3. The simulation output is the same as described above. For each simulation a separate output folder is created in the ``history`` folder.
diff --git a/docs/dashboard.rst b/docs/dashboard.rst
new file mode 100644
index 00000000..e0cd0a68
--- /dev/null
+++ b/docs/dashboard.rst
@@ -0,0 +1,5 @@
+Comparison of different models using Dash
+=========================================
+You can compare the result of different models using the dashboard build with `Dash <https://dash.plotly.com>`_. To create the dashboard, run the ``mpet_plot_app.py`` script (located in the bin folder) and use the ``-d`` argument to provide a directory with model outputs to include the dashbaord. Each model output should be located in a subfolder of the provided directory. For example, to compare the results of all models saved in subfolders of the folder ``history``, run the command:
+``bin/mpet_plot_app.py -d history``. It will try to open the dashbaord in your web browser, where the different models can be identified based on their subfolder name.
+Running this script requires the following packages to be installed: `dash <https://pypi.org/project/dash/>`_ and `dash_bootstrap_components <https://pypi.org/project/dash-bootstrap-components/>`_.
diff --git a/docs/ensemble.rst b/docs/ensemble.rst
new file mode 100644
index 00000000..afeb94bc
--- /dev/null
+++ b/docs/ensemble.rst
@@ -0,0 +1,6 @@
+Parameter space exploration
+===========================
+
+The script ``./bin/ensemble.py`` can be used to to automatically generate a set of config files in which one or more parameters can be varied over specified values. The script is executed by running ``./bin/ensemble.py "reference_config_file"``. Within the enseble.py you can indicate which parameter(s) you want to explore over what range. The output will be a number of config files in which all combinations of the given parameter values are used.
+
+Additionally the overarching script ``/bin/mpet_create_runjobs_dashboard.py`` combines the enseble creation, cluster run, and the dashboard. In this script you can idicate the reference config file and the parameter space you want to explore, and the settings of the cluster script you want to use. Next, execute the shell script ``bin/run_mpet_create_run_dashboard.sh``. It will create the config files in the range of inidcated parameters, execute the run cluster script for all created config files, and plot the results of all these models in the dashbaord (which will be started automatically).
diff --git a/docs/index.rst b/docs/index.rst
index 18ae359a..a308e6aa 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -18,7 +18,10 @@ Welcome to MPET's documentation!
    :caption: Tutorials:
 
    run_sim
+   cluster
    analysis
+   dashboard
+   ensemble
 
 .. toctree::
    :maxdepth: 1
diff --git a/mpet/config/configuration.py b/mpet/config/configuration.py
index 2f959267..ebe39704 100644
--- a/mpet/config/configuration.py
+++ b/mpet/config/configuration.py
@@ -140,10 +140,6 @@ def _init_from_cfg(self, paramfile):
         if self.D_s['Nvol_a'] > 0:
             trodes.append('a')
         self['trodes'] = trodes
-        # to check for separator, directly access underlying dict of system config;
-        # self['Nvol']['s'] would not work because that requires have_separator to
-        # be defined already
-        self['have_separator'] = self.D_s['Nvol_s'] > 0
 
         # load electrode parameter file(s)
         self.paramfiles = {}
@@ -417,6 +413,7 @@ def _process_config(self, prevDir=None):
         #. Scale to non-dimensional values
         #. Parse current/voltage segments
         #. Either generate particle distributions or load from previous run
+        #. Create simulation times
 
         :param bool prevDir: if True, load particle distributions from previous run,
             otherwise generate them
@@ -440,6 +437,10 @@ def _process_config(self, prevDir=None):
             # set to theoretical value
             self[param] = theoretical_1C_current
 
+        # use custom concentration when using solid electrolyte model
+        if self['elyteModelType'] == 'solid':
+            self['c0'] = self['c0']
+
         # apply scalings
         self._scale_system_parameters(theoretical_1C_current)
         self._scale_electrode_parameters()  # includes separator
@@ -463,6 +464,16 @@ def _process_config(self, prevDir=None):
             # Electrode parameters that depend on invidividual particle
             self._indvPart()
 
+        self._create_times()
+
+    def _create_times(self):
+        """
+
+        """
+        # The list of reporting times excludes the first index (zero, which is implied)
+        if not self["times"]:
+            self["times"] = list(np.linspace(0, self["tend"], self["tsteps"] + 1))[1:]
+
     def _scale_system_parameters(self, theoretical_1C_current):
         """
         Scale system parameters to non-dimensional values. This method should be called only once,
@@ -483,9 +494,19 @@ def _scale_system_parameters(self, theoretical_1C_current):
         self['k0_foil'] = self['k0_foil'] / (self['1C_current_density'] * self['curr_ref'])
         self['Rfilm_foil'] = self['Rfilm_foil'] / self['Rser_ref']
 
+        if self['elyteModelType'] == 'solid':
+            self['cmax'] = self['cmax'] / constants.c_ref
+
+        if self['simInterface_a'] or self['simInterface_c']:
+            self['Dp_i'] = self['Dp_i'] / self['D_ref']
+            if self['interfaceModelType'] != 'solid':
+                self['Dm_i'] = self['Dm_i'] / self['D_ref']
+            self['c0_int'] = self['c0_int'] / constants.c_ref
+            self['cmax_i'] = self['cmax_i'] / constants.c_ref
+
     def _scale_electrode_parameters(self):
         """
-        Scale electrode and separator parameters to non-dimensional values.
+        Scale electrode, separator, and interface region parameters to non-dimensional values.
         This method should be called only once, from :meth:`_process_config`.
         """
         kT = constants.k * constants.T_ref
@@ -506,9 +527,13 @@ def _scale_electrode_parameters(self):
                     self[trode, param] = value / kT
 
         # scalings on separator
-        if self['have_separator']:
+        if self['Nvol']['s']:
             self['L']['s'] /= self['L_ref']
 
+        # scaling related to interface region (if present)
+        if self['simInterface_a'] or self['simInterface_c']:
+            self['L_i'] /= self['L_ref']
+
     def _scale_macroscopic_parameters(self, Vref):
         """
         Scale macroscopic input parameters to non-dimensional values and add
@@ -606,6 +631,8 @@ def _distr_part(self):
         self['psd_vol'] = {}
         self['psd_vol_FracVol'] = {}
 
+        self['gamma_contact'] = {}
+
         for trode in self['trodes']:
             solidType = self[trode, 'type']
             Nvol = self['Nvol'][trode]
@@ -631,10 +658,27 @@ def _distr_part(self):
                     raise ValueError('Specified particle size distribution discretization '
                                      'of volumes inequal to the one specified in the config file')
 
+            stddev_c = self['stand_dev_contact']
+            mean_c = self['fraction_of_contact']
+
+            if 0 < mean_c < 1:
+                # Contact penalty for BV
+                mean_c = 1 - mean_c  # to make distribution start at 1 if gamma is 1
+                var_c = stddev_c**2
+                mu_c = np.log((mean_c**2) / np.sqrt(var_c + mean_c**2))
+                sigma_c = np.sqrt(np.log(var_c/(mean_c**2) + 1))
+                raw_c = 1 - np.random.lognormal(mu_c, sigma_c, size=(Nvol, Npart))
+                raw_c[raw_c <= 0.0001] = 0.0001
+                gamma_contact = raw_c
+            elif mean_c == 1:
+                gamma_contact = np.ones((Nvol, Npart))
+            else:
+                raise NotImplementedError('Contact error should be between 0 and 1')
+
             # For particles with internal profiles, convert psd to
             # integers -- number of steps
             solidDisc = self[trode, 'discretization']
-            if solidType in ['ACR']:
+            if solidType in ['ACR','ACR_Diff','ACR2']:
                 psd_num = np.ceil(raw / solidDisc).astype(int)
                 psd_len = solidDisc * psd_num
             elif solidType in ['CHR', 'diffn', 'CHR2', 'diffn2']:
@@ -674,6 +718,7 @@ def _distr_part(self):
             self['psd_area'][trode] = psd_area
             self['psd_vol'][trode] = psd_vol
             self['psd_vol_FracVol'][trode] = psd_frac_vol
+            self['gamma_contact'][trode] = gamma_contact
 
     def _G(self):
         """
@@ -721,6 +766,7 @@ def _indvPart(self):
                     plen = self['psd_len'][trode][i,j]
                     parea = self['psd_area'][trode][i,j]
                     pvol = self['psd_vol'][trode][i,j]
+                    gamma_cont = self['gamma_contact'][trode][i,j]
                     # Define a few reference scales
                     F_s_ref = plen * cs_ref_part / self['t_ref']  # part/(m^2 s)
                     i_s_ref = constants.e * F_s_ref  # A/m^2
@@ -738,6 +784,10 @@ def _indvPart(self):
                         / (constants.k * constants.N_A * constants.T_ref)
                     self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \
                         / (constants.e * F_s_ref)
+                    self[trode, 'indvPart']['gamma_con'][i, j] = gamma_cont
+                    if self['fraction_of_contact'] != 1.0 and not self['localized_losses']:
+                        self[trode, 'indvPart']['k0'][i, j] = self[trode, 'k0'] \
+                            / (constants.e * F_s_ref)*gamma_cont
                     self[trode, 'indvPart']['E_A'][i, j] = self[trode, 'E_A'] \
                         / (constants.k * constants.N_A * constants.T_ref)
                     self[trode, 'indvPart']['Rfilm'][i, j] = self[trode, 'Rfilm'] \
@@ -763,8 +813,8 @@ def _verify_config(self):
         for trode in self['trodes']:
             solidShape = self[trode, 'shape']
             solidType = self[trode, 'type']
-            if solidType in ["ACR", "homog_sdn"] and solidShape != "C3":
-                raise Exception("ACR and homog_sdn req. C3 shape")
+            if solidType in ["ACR", "ACR_Diff", "homog_sdn", "ACR2"] and solidShape != "C3":
+                raise Exception("ACR, ACR_Diff, ACR2 and homog_sdn req. C3 shape")
             if (solidType in ["CHR", "diffn"] and solidShape not in ["sphere", "cylinder"]):
                 raise NotImplementedError("CHR and diffn req. sphere or cylinder")
 
diff --git a/mpet/config/constants.py b/mpet/config/constants.py
index ed309c15..32fcc1a1 100644
--- a/mpet/config/constants.py
+++ b/mpet/config/constants.py
@@ -11,9 +11,9 @@
 #: Reference flux, C/mol
 F = e * N_A
 #: General particle classification (1 var)
-two_var_types = ["diffn2", "CHR2", "homog2", "homog2_sdn"]
+two_var_types = ["diffn2", "CHR2", "homog2", "homog2_sdn", "ACR2"]
 #: General particle classification (2 var)
-one_var_types = ["ACR", "diffn", "CHR", "homog", "homog_sdn"]
+one_var_types = ["ACR","ACR_Diff", "diffn", "CHR", "homog", "homog_sdn"]
 #: Reference concentration, mol/m^3 = 1M
 c_ref = 1000.
 #: Reaction rate epsilon for values close to zero
@@ -22,10 +22,10 @@
 #: parameter that are defined per electrode with a ``_{electrode}`` suffix
 PARAMS_PER_TRODE = ['Nvol', 'Npart', 'mean', 'stddev', 'cs0', 'simBulkCond', 'sigma_s',
                     'simPartCond', 'G_mean', 'G_stddev', 'L', 'P_L', 'poros', 'BruggExp',
-                    'specified_psd']
+                    'specified_psd','specified_poros']
 #: subset of ``PARAMS_PER_TRODE``` that is defined for the separator as well
-PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp']
+PARAMS_SEPARATOR = ['Nvol', 'L', 'poros', 'BruggExp', 'specified_poros']
 #: parameters that are defined for each particle, and their type
 PARAMS_PARTICLE = {'N': int, 'kappa': float, 'beta_s': float, 'D': float, 'k0': float,
                    'Rfilm': float, 'delta_L': float, 'Omega_a': float, 'E_D': float,
-                   'E_A': float}
+                   'E_A': float, 'gamma_con': float}
diff --git a/mpet/config/derived_values.py b/mpet/config/derived_values.py
index fb18dd11..bac31afc 100644
--- a/mpet/config/derived_values.py
+++ b/mpet/config/derived_values.py
@@ -173,8 +173,16 @@ def csmax(self, trode):
     def cap(self, trode):
         """Theoretical capacity of given electrode
         """
-        return constants.e * self.config['L'][trode] * (1 - self.config['poros'][trode]) \
+        capacity = constants.e * self.config['L'][trode] * (1 - self.config['poros'][trode]) \
             * self.config['P_L'][trode] * self.config[trode, 'rho_s']
+        if np.all(self.config['specified_poros'][trode]):
+            sum_porosity = 0
+            for i in self.config['specified_poros'][trode]:
+                sum_porosity = sum_porosity + i
+            avg_porosity = sum_porosity/np.size(self.config['specified_poros'][trode])
+            capacity = constants.e * self.config['L'][trode] * (1 - avg_porosity) \
+                * self.config['P_L'][trode] * self.config[trode, 'rho_s']
+        return capacity
 
     def numsegments(self):
         """Number of segments in voltage/current profile
@@ -191,6 +199,8 @@ def D_ref(self):
         """
         if self.config['elyteModelType'] == 'dilute':
             return self.config['Damb']
+        elif self.config['elyteModelType'] == 'solid':
+            return self.config['Damb']
         else:
             SMset = self.config["SMset"]
             elyte_function = import_function(self.config["SMset_filename"], SMset,
diff --git a/mpet/config/parameterset.py b/mpet/config/parameterset.py
index a4d1556b..65673490 100644
--- a/mpet/config/parameterset.py
+++ b/mpet/config/parameterset.py
@@ -87,7 +87,7 @@ def __getitem__(self, item):
             d = {}
             # some parameters are also defined for the separator
             trodes = self['trodes'][:]  # make a copy here to avoid adding values to the original
-            if item in PARAMS_SEPARATOR and self['have_separator']:
+            if item in PARAMS_SEPARATOR:
                 trodes.append('s')
             for trode in trodes:
                 # get the value for this electrode/separator and store it
diff --git a/mpet/config/schemas.py b/mpet/config/schemas.py
index 17e2ea76..8a1df630 100644
--- a/mpet/config/schemas.py
+++ b/mpet/config/schemas.py
@@ -51,6 +51,10 @@ def tobool(value):
     return bool(strtobool(value))
 
 
+#: Defaults for config sections that are optional
+DEFAULT_SECTIONS = {'Interface': {'simInterface_a': False, 'simInterface_c': False}}
+
+
 #: System parameters, per section
 system = {'Sim Params': {'profileType': lambda x:
                          check_allowed_values(x, ["CC", "CV", "CP", "CCsegments", "CVsegments"]),
@@ -66,6 +70,7 @@ def tobool(value):
                          Optional('prevDir', default=''): str,
                          'tend': And(Use(float), lambda x: x > 0),
                          'tsteps': And(Use(int), lambda x: x > 0),
+                         Optional('times', default=[]): Use(ast.literal_eval),
                          'relTol': And(Use(float), lambda x: x > 0),
                          'absTol': And(Use(float), lambda x: x > 0),
                          'T': Use(float),
@@ -86,6 +91,10 @@ def tobool(value):
                         'stddev_c': Use(float),
                         'mean_a': Use(float),
                         'stddev_a': Use(float),
+                        Optional('fraction_of_contact',default=1.0): Use(float),
+                        Optional('stand_dev_contact',default=0): Use(float),
+                        Optional('localized_losses', default=False):
+                            Or(Use(tobool), Use(lambda x: np.array(ast.literal_eval(x)))),
                         'cs0_c': Use(float),
                         'cs0_a': Use(float),
                         Optional('specified_psd_c', default=False):
@@ -110,6 +119,12 @@ def tobool(value):
                        'poros_c': Use(float),
                        'poros_a': Use(float),
                        'poros_s': Use(float),
+                       Optional('specified_poros_c', default=False):
+                            Or(Use(tobool), Use(lambda x: np.array(ast.literal_eval(x)))),
+                       Optional('specified_poros_a', default=False):
+                            Or(Use(tobool), Use(lambda x: np.array(ast.literal_eval(x)))),
+                       Optional('specified_poros_s', default=False):
+                            Or(Use(tobool), Use(lambda x: np.array(ast.literal_eval(x)))),
                        'BruggExp_c': Use(float),
                        'BruggExp_a': Use(float),
                        'BruggExp_s': Use(float)},
@@ -124,7 +139,21 @@ def tobool(value):
                           'n': Use(int),
                           'sp': Use(int),
                           Optional('Dp', default=None): Use(float),
-                          Optional('Dm', default=None): Use(float)}}
+                          Optional('Dm', default=None): Use(float),
+                          Optional('cmax'): Use(float),
+                          Optional('a_slyte'): Use(float)},
+          'Interface': {Optional('simInterface_a',default=False): Use(tobool),
+                        Optional('simInterface_c',default=False): Use(tobool),
+                        Optional('Nvol_i'): And(Use(int), lambda x: x > 0),
+                        Optional('L_i'): Use(float),
+                        Optional('BruggExp_i'): Use(float),
+                        Optional('poros_i'): Use(float),
+                        Optional('interfaceModelType'): str,
+                        Optional('interfaceSMset'): str,
+                        Optional('c0_int'): Use(float),
+                        Optional('cmax_i'): Use(float),
+                        Optional('Dp_i'): Use(float),
+                        Optional('Dm_i'): Use(float)}}
 
 #: Electrode parameters, per section
 electrode = {'Particles': {'type': lambda x: check_allowed_values(x,
diff --git a/mpet/electrode/materials/LTO.py b/mpet/electrode/materials/LTO.py
new file mode 100644
index 00000000..cf1a25fd
--- /dev/null
+++ b/mpet/electrode/materials/LTO.py
@@ -0,0 +1,14 @@
+import numpy as np
+
+
+def LTO(self, y, ybar, muR_ref, ISfuncs=None):
+    """
+    Vasileiadis 2017
+    """
+    muRtheta = -self.eokT * 1.55
+    muRhomog = self.reg_sln(y, self.get_trode_param("Omega_a"))
+    muRnonHomog = self.general_non_homog(y, ybar)
+    muR = muRhomog + muRnonHomog
+    actR = np.exp(muR / self.T)
+    muR += muRtheta + muR_ref
+    return muR, actR
diff --git a/mpet/electrode/reactions/BV_NMC_Park2021.py b/mpet/electrode/reactions/BV_NMC_Park2021.py
new file mode 100644
index 00000000..a3e489a3
--- /dev/null
+++ b/mpet/electrode/reactions/BV_NMC_Park2021.py
@@ -0,0 +1,10 @@
+import numpy as np
+
+
+def BV_NMC_Park2021(eta, c_sld, c_lyte, k0, E_A, T, act_R=None,
+                    act_lyte=None, lmbda=None, alpha=None):
+    """Implemented for the Park, J 2021 model
+    for the NMC electrode"""
+    ecd = k0*act_lyte/(1+np.exp(25*(c_sld-0.6)))
+    Rate = ecd * np.exp(-E_A/T + E_A/1) * (np.exp(-alpha*eta/T) - np.exp((1-alpha)*eta/T))
+    return Rate
diff --git a/mpet/electrolyte/solid_elyte.py b/mpet/electrolyte/solid_elyte.py
new file mode 100644
index 00000000..41cc5989
--- /dev/null
+++ b/mpet/electrolyte/solid_elyte.py
@@ -0,0 +1,15 @@
+from .valoen_reimers import valoen_reimers
+
+
+def solid_elyte():
+    """
+    Solid Electrolyte version
+    """
+    tp0 = 0.99999
+    D = 1.19e-11  # m^2/s
+
+    Ign1, sigma_ndim, thermFac, Ign3, Dref = valoen_reimers()
+    D_ndim = D / Dref
+
+    # D_ndim and tp0 are constants, but a callable must be returned
+    return lambda c: D_ndim, sigma_ndim, thermFac, tp0, Dref
diff --git a/mpet/geometry.py b/mpet/geometry.py
index c35dbba1..52eef752 100644
--- a/mpet/geometry.py
+++ b/mpet/geometry.py
@@ -91,3 +91,25 @@ def get_elyte_disc(Nvol, L, poros, BruggExp):
     out["eps_o_tau"] = porosvec_pad/porosvec_pad**(Brugg_pad)
 
     return out
+
+
+def get_interface_disc(Nvol, L, poros, BruggExp):
+    out = {}
+    # Width of each cell
+    out["dxvec"] = utils.get_const_vec(L/Nvol, Nvol)
+
+    # Distance between cell centers
+    dxtmp = np.hstack((out["dxvec"][0], out["dxvec"], out["dxvec"][-1]))
+    out["dxd1"] = utils.mean_linear(dxtmp)
+
+    # The porosity vector
+    out["porosvec"] = utils.get_const_vec(poros, Nvol)
+    porosvec_pad = utils.pad_vec(out["porosvec"])
+
+    # Vector of Bruggeman exponents
+    Brugg_pad = utils.pad_vec(utils.get_const_vec(BruggExp, Nvol))
+
+    # Vector of posority/tortuosity (assuming Bruggeman)
+    out["eps_o_tau"] = porosvec_pad/porosvec_pad**(Brugg_pad)
+
+    return out
diff --git a/mpet/main.py b/mpet/main.py
index 3ce8c68d..e6ed890b 100644
--- a/mpet/main.py
+++ b/mpet/main.py
@@ -10,7 +10,6 @@
 
 import daetools.pyDAE as dae
 from daetools.solvers.superlu import pySuperLU
-import numpy as np
 
 import mpet
 import mpet.data_reporting as data_reporting
@@ -37,13 +36,19 @@ def run_simulation(config, outdir):
     # Turn off reporting of some variables
     simulation.m.endCondition.ReportingOn = False
 
-    # Turn off reporting of particle ports
+    # Turn off reporting of particle and interface 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
+        if config[f"simInterface_{trode}"]:
+            for interfaces in simulation.m.interfaces[trode]:
+                for iModel in interfaces:
+                    for port in iModel.Ports:
+                        for var in port.Variables:
+                            var.ReportingOn = False
 
     # Turn off reporting of cell ports
     for port in simulation.m.Ports:
@@ -56,7 +61,7 @@ def run_simulation(config, outdir):
     # Set the time horizon and the reporting interval
     simulation.TimeHorizon = config["tend"]
     # The list of reporting times excludes the first index (zero, which is implied)
-    simulation.ReportingTimes = list(np.linspace(0, config["tend"], config["tsteps"] + 1))[1:]
+    simulation.ReportingTimes = config["times"]
 
     # Connect data reporter
     simName = simulation.m.Name + time.strftime(
@@ -69,7 +74,8 @@ def run_simulation(config, outdir):
 
     # Solve at time=0 (initialization)
     # Increase the number of Newton iterations for more robust initialization
-    dae.daeGetConfig().SetString("daetools.IDAS.MaxNumItersIC","100")
+    dae.daeGetConfig().SetString("daetools.IDAS.MaxNumItersIC","1000")
+    dae.daeGetConfig().SetString("daetools.IDAS.MaxNumSteps","100000")
     simulation.SolveInitial()
 
     # Run
@@ -86,14 +92,16 @@ def run_simulation(config, outdir):
     simulation.Finalize()
 
 
-def main(paramfile, keepArchive=True):
+def main(paramfile, keepArchive=True, keepFullRun=False):
     timeStart = time.time()
     # Get the parameters dictionary (and the config instance) from the
     # parameter file
     config = Config(paramfile)
 
     # Directories we'll store output in.
-    outdir_name = time.strftime("%Y%m%d_%H%M%S", time.localtime())
+    config_file = os.path.basename(paramfile)
+    config_base = os.path.splitext(config_file)[0]
+    outdir_name = "_".join((time.strftime("%Y%m%d_%H%M%S", time.localtime()), config_base))
     outdir_path = os.path.join(os.getcwd(), "history")
     outdir = os.path.join(outdir_path, outdir_name)
     # Make sure there's a place to store the output
@@ -194,10 +202,16 @@ def main(paramfile, keepArchive=True):
     except Exception:
         pass
 
-    # Copy or move simulation output to current directory
+    # Copy or move simulation output to current directory. If running multiple jobs,
+    # make sure to keep all sim_output
     tmpDir = os.path.join(os.getcwd(), "sim_output")
-    shutil.rmtree(tmpDir, ignore_errors=True)
+    if not keepFullRun:
+        shutil.rmtree(tmpDir, ignore_errors=True)
+        tmpsubDir = tmpDir
+    else:
+        tmpsubDir = os.path.join(tmpDir, outdir_name)
+
     if keepArchive:
-        shutil.copytree(outdir, tmpDir)
+        shutil.copytree(outdir, tmpsubDir)
     else:
-        shutil.move(outdir, tmpDir)
+        shutil.move(outdir, tmpsubDir)
diff --git a/mpet/mod_cell.py b/mpet/mod_cell.py
index 46d244d2..97533395 100644
--- a/mpet/mod_cell.py
+++ b/mpet/mod_cell.py
@@ -15,6 +15,7 @@
 import mpet.extern_funcs as extern_funcs
 import mpet.geometry as geom
 import mpet.mod_electrodes as mod_electrodes
+from mpet.mod_interface import InterfaceRegion
 import mpet.ports as ports
 import mpet.utils as utils
 from mpet.config import constants
@@ -39,7 +40,7 @@ def __init__(self, config, Name, Parent=None, Description=""):
         # Domains where variables are distributed
         self.DmnCell = {}  # domains over full cell dimensions
         self.DmnPart = {}  # domains over particles in each cell volume
-        if config['have_separator']:  # If we have a separator
+        if config['Nvol']['s']:  # If we have a separator
             self.DmnCell["s"] = dae.daeDomain(
                 "DmnCell_s", self, dae.unit(),
                 "Simulated volumes in the separator")
@@ -58,6 +59,7 @@ def __init__(self, config, Name, Parent=None, Description=""):
         self.phi_bulk = {}
         self.phi_part = {}
         self.R_Vp = {}
+        self.R_Vi = {}
         self.ffrac = {}
         for trode in trodes:
             # Concentration/potential in electrode regions of elyte
@@ -81,10 +83,15 @@ def __init__(self, config, Name, Parent=None, Description=""):
                 "R_Vp_{trode}".format(trode=trode), dae.no_t, self,
                 "Rate of reaction of positives per electrode volume",
                 [self.DmnCell[trode]])
+            if self.config[f'simInterface_{trode}']:
+                self.R_Vi[trode] = dae.daeVariable(
+                    "R_Vi_{trode}".format(trode=trode), dae.no_t, self,
+                    "Rate of reaction of positives per electrode volume with interface region",
+                    [self.DmnCell[trode]])
             self.ffrac[trode] = dae.daeVariable(
                 "ffrac_{trode}".format(trode=trode), mole_frac_t, self,
                 "Overall filling fraction of solids in electrodes")
-        if config['have_separator']:  # If we have a separator
+        if config['Nvol']['s']:  # If we have a separator
             self.c_lyte["s"] = dae.daeVariable(
                 "c_lyte_s", conc_t, self,
                 "Concentration in the electrolyte in the separator",
@@ -96,7 +103,7 @@ def __init__(self, config, Name, Parent=None, Description=""):
         # Note if we're doing a single electrode volume simulation
         # It will be in a perfect bath of electrolyte at the applied
         # potential.
-        if ('a' not in config['trodes']) and (not config['have_separator']) and Nvol["c"] == 1:
+        if ('a' not in config['trodes']) and (not config['Nvol']['s']) and Nvol["c"] == 1:
             self.SVsim = True
         else:
             self.SVsim = False
@@ -119,13 +126,17 @@ def __init__(self, config, Name, Parent=None, Description=""):
         # volumes and ports with which to talk to them.
         self.portsOutLyte = {}
         self.portsOutBulk = {}
+        self.portsInInterface = {}
         self.particles = {}
+        self.interfaces = {}
         for trode in trodes:
             Nv = Nvol[trode]
             Np = Npart[trode]
             self.portsOutLyte[trode] = np.empty(Nv, dtype=object)
             self.portsOutBulk[trode] = np.empty((Nv, Np), dtype=object)
+            self.portsInInterface[trode] = np.empty((Nv, Np), dtype=object)
             self.particles[trode] = np.empty((Nv, Np), dtype=object)
+            self.interfaces[trode] = np.empty((Nv, Np), dtype=object)
             for vInd in range(Nv):
                 self.portsOutLyte[trode][vInd] = ports.portFromElyte(
                     "portTrode{trode}vol{vInd}".format(trode=trode, vInd=vInd), dae.eOutletPort,
@@ -148,8 +159,42 @@ def __init__(self, config, Name, Parent=None, Description=""):
                         Name="partTrode{trode}vol{vInd}part{pInd}".format(
                             trode=trode, vInd=vInd, pInd=pInd),
                         Parent=self)
-                    self.ConnectPorts(self.portsOutLyte[trode][vInd],
-                                      self.particles[trode][vInd,pInd].portInLyte)
+
+                    if config[f"simInterface_{trode}"]:
+                        # instantiate interfaces between particle and electrolyte per particle
+                        self.interfaces[trode][vInd,pInd] = InterfaceRegion(
+                            Name="interfaceTrode{trode}vol{vInd}part{pInd}".format(
+                                trode=trode, vInd=vInd, pInd=pInd),
+                            Parent=self, config=config, cell=self,
+                            particle=self.particles[trode][vInd,pInd],
+                            vInd=vInd,pInd=pInd,trode=trode)
+
+                        self.portsInInterface[trode][vInd,pInd] = ports.portFromInterface(
+                            "portIface{trode}vol{vInd}part{pInd}".format(
+                                trode=trode, vInd=vInd, pInd=pInd),
+                            dae.eInletPort, self,
+                            "Interface region port to elyte")
+
+                        # connect elyte to interface, then interface to particle
+                        self.ConnectPorts(self.portsOutLyte[trode][vInd],
+                                          self.interfaces[trode][vInd,pInd].portInLyte)
+
+                        self.ConnectPorts(
+                            self.interfaces[trode][vInd,pInd].portOutInterfaceParticle,
+                            self.particles[trode][vInd,pInd].portInLyte)
+
+                        # connect interface to elyte
+                        self.ConnectPorts(self.interfaces[trode][vInd,pInd].portOutInterfaceElyte,
+                                          self.portsInInterface[trode][vInd,pInd])
+
+                        # connect particle to interface
+                        self.ConnectPorts(self.particles[trode][vInd,pInd].portOutParticle,
+                                          self.interfaces[trode][vInd,pInd].portInParticle)
+                    else:
+                        # connect elyte to particle
+                        self.ConnectPorts(self.portsOutLyte[trode][vInd],
+                                          self.particles[trode][vInd,pInd].portInLyte)
+
                     self.ConnectPorts(self.portsOutBulk[trode][vInd,pInd],
                                       self.particles[trode][vInd,pInd].portInBulk)
 
@@ -186,6 +231,11 @@ def DeclareEquations(self):
                 # Start with no reaction, then add reactions for each
                 # particle in the volume.
                 RHS = 0
+                # interface region has separate reaction rate
+                if config[f"simInterface_{trode}"]:
+                    eq_i = self.CreateEquation(
+                        "R_Vi_trode{trode}vol{vInd}".format(vInd=vInd, trode=trode))
+                    RHS_i = 0
                 # sum over particle volumes in given electrode volume
                 for pInd in range(Npart[trode]):
                     # The volume of this particular particle
@@ -193,7 +243,14 @@ def DeclareEquations(self):
                     RHS += -(config["beta"][trode] * (1-config["poros"][trode])
                              * config["P_L"][trode] * Vj
                              * self.particles[trode][vInd,pInd].dcbardt())
+                    if config[f"simInterface_{trode}"]:
+                        # Nm0 = self.portsInInterface[trode][vInd,pInd].Nm0()
+                        i0 = self.portsInInterface[trode][vInd,pInd].i0()
+                        # TODO: what is the reaction rate?
+                        RHS_i += -i0
                 eq.Residual = self.R_Vp[trode](vInd) - RHS
+                if config[f"simInterface_{trode}"]:
+                    eq_i.Residual = self.R_Vi[trode](vInd) - RHS_i
 
         # Define output port variables
         for trode in trodes:
@@ -221,6 +278,12 @@ def DeclareEquations(self):
                 phi_tmp = utils.add_gp_to_vec(utils.get_var_vec(self.phi_bulk[trode], Nvol[trode]))
                 porosvec = utils.pad_vec(utils.get_const_vec(
                     (1-self.config["poros"][trode])**(1-config["BruggExp"][trode]), Nvol[trode]))
+                if np.all(self.config['specified_poros'][trode]):
+                    specified_por = self.config['specified_poros'][trode]
+                    porosvec = (
+                        (np.ones(Nvol[trode]) - specified_por))**(
+                            (1 - self.config["BruggExp"][trode]))
+                    porosvec = utils.pad_vec(porosvec)
                 poros_walls = utils.mean_harmonic(porosvec)
                 if trode == "a":  # anode
                     # Potential at the current collector is from
@@ -236,12 +299,18 @@ def DeclareEquations(self):
                 dx = config["L"][trode]/Nvol[trode]
                 dvg_curr_dens = np.diff(-poros_walls*config["sigma_s"][trode]
                                         * np.diff(phi_tmp)/dx)/dx
+
             # Actually set up the equations for bulk solid phi
             for vInd in range(Nvol[trode]):
                 eq = self.CreateEquation(
                     "phi_ac_trode{trode}vol{vInd}".format(vInd=vInd, trode=trode))
                 if simBulkCond:
-                    eq.Residual = -dvg_curr_dens[vInd] - self.R_Vp[trode](vInd)
+                    # select reaction rate with interface region or particle
+                    if config[f"simInterface_{trode}"]:
+                        R_V = self.R_Vi
+                    else:
+                        R_V = self.R_Vp
+                    eq.Residual = -dvg_curr_dens[vInd] - R_V[trode](vInd)
                 else:
                     if trode == "a":  # anode
                         eq.Residual = self.phi_bulk[trode](vInd) - self.phi_cell()
@@ -287,11 +356,18 @@ def DeclareEquations(self):
             eq = self.CreateEquation("phi_lyte")
             eq.Residual = self.phi_lyte["c"](0) - self.phi_cell()
         else:
-            disc = geom.get_elyte_disc(Nvol, config["L"], config["poros"], config["BruggExp"])
+            if np.all(config["specified_poros"]["c"]):
+                config_poros = config["specified_poros"]
+            else:
+                config_poros = config["poros"]
+            disc = geom.get_elyte_disc(Nvol, config["L"], config_poros, config["BruggExp"])
             cvec = utils.get_asc_vec(self.c_lyte, Nvol)
             dcdtvec = utils.get_asc_vec(self.c_lyte, Nvol, dt=True)
             phivec = utils.get_asc_vec(self.phi_lyte, Nvol)
-            Rvvec = utils.get_asc_vec(self.R_Vp, Nvol)
+            if config[f"simInterface_{trode}"]:
+                Rvvec = utils.get_asc_vec(self.R_Vi, Nvol)
+            else:
+                Rvvec = utils.get_asc_vec(self.R_Vp, Nvol)
             # Apply concentration and potential boundary conditions
             # Ghost points on the left and no-gradients on the right
             ctmp = np.hstack((self.c_lyteGP_L(), cvec, cvec[-1]))
@@ -313,18 +389,16 @@ def DeclareEquations(self):
                 # exchange current density, ecd = k0_foil * c_lyte**(0.5)
                 cWall = .5*(ctmp[0] + ctmp[1])
                 ecd = config["k0_foil"]*cWall**0.5
+                # Concentration is fixed for solid
+                if config["elyteModelType"] == 'solid':
+                    ecd = config["k0_foil"]*1**0.5
                 # note negative current because positive current is
                 # oxidation here
-                BVfunc = -self.current() / ecd
-                eta_eff = 2*np.arcsinh(-BVfunc/2.)
-                eta = eta_eff + self.current()*config["Rfilm_foil"]
-                # eta = mu_R - mu_O = -mu_O (evaluated at interface)
-                # mu_O = [T*ln(c) +] phiWall - phi_cell = -eta
-                # phiWall = -eta + phi_cell [- T*ln(c)]
-                phiWall = -eta + self.phi_cell()
+                eta = self.phi_cell() - self.current()*config["Rfilm_foil"] \
+                    - .5*(phitmp[0] + phitmp[1])
                 if config["elyteModelType"] == "dilute":
-                    phiWall -= config["T"]*np.log(cWall)
-                eqP.Residual = phiWall - .5*(phitmp[0] + phitmp[1])
+                    eta -= config["T"]*np.log(cWall)
+                eqP.Residual = self.current() - ecd*2*np.sinh(eta/2)
 
             # We have a porous anode -- no flux of charge or anions through current collector
             else:
@@ -355,6 +429,7 @@ def DeclareEquations(self):
                 eq.Residual -= dx * self.R_Vp[limtrode](vInd)/rxn_scl
             else:
                 eq.Residual += dx * self.R_Vp[limtrode](vInd)/rxn_scl
+
         # Define the measured voltage, offset by the "applied" voltage
         # by any series resistance.
         # phi_cell = phi_applied - I*R
@@ -482,6 +557,7 @@ def get_lyte_internal_fluxes(c_lyte, phi_lyte, disc, config):
 
     # 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 config["elyteModelType"] == "dilute":
@@ -514,4 +590,26 @@ def get_lyte_internal_fluxes(c_lyte, phi_lyte, disc, config):
             )
         Nm_edges_int = num*(-D_edges*np.diff(c_lyte)/dxd1
                             + (1./(num*zm)*(1-tp0(c_edges_int, T))*i_edges_int))
+    elif config["elyteModelType"] == "solid":
+        SMset = config["SMset"]
+        elyte_function = utils.import_function(config["SMset_filename"], SMset,
+                                               mpet_module=f"mpet.electrolyte.{SMset}")
+        D_fs, sigma_fs, thermFac, tp0 = elyte_function()[:-1]
+        # sigma_fs and thermFac not used bc the solid system is considered linear
+        a_slyte = config["a_slyte"]
+        c_edges_int_norm = c_edges_int / config["cmax"]
+
+        # Get diffusivity at cell edges using weighted harmonic mean
+        eps_o_tau_edges = utils.weighted_linear_mean(eps_o_tau, wt)
+        Dp = eps_o_tau_edges * config["Dp"]
+        Dm = (zp * Dp - zp * Dp * tp0) / (tp0 * zm)
+        Dp0 = Dp / (1-c_edges_int_norm)  # should be c0/cmax
+        Dchemp = Dp0 * (1 - 2 * a_slyte * c_edges_int_norm + 2 * a_slyte * c_edges_int_norm**2)
+        Dchemm = Dm
+        Damb = (zp * Dp * Dchemm + zm * Dm * Dchemp) / (zp * Dp - zm * Dm)
+        i_edges_int = (-((nup*zp*Dchemp + num*zm*Dchemm)*np.diff(c_lyte)/dxd1)
+                       - (nup * zp ** 2 * Dp0 * (1 - c_edges_int_norm) + num * zm ** 2 * Dm) / T
+                       * c_edges_int * np.diff(phi_lyte) / dxd1)
+        Nm_edges_int = num * (-Damb * np.diff(c_lyte) / dxd1
+                              + (1. / (num * zm) * (1 - tp0) * i_edges_int))
     return Nm_edges_int, i_edges_int
diff --git a/mpet/mod_electrodes.py b/mpet/mod_electrodes.py
index bbb241db..b6e2d13f 100644
--- a/mpet/mod_electrodes.py
+++ b/mpet/mod_electrodes.py
@@ -11,7 +11,10 @@
 These models can be instantiated from the mod_cell module to simulate various types of active
 materials within a battery electrode.
 """
+
+
 import daetools.pyDAE as dae
+
 import numpy as np
 import scipy.sparse as sprs
 import scipy.interpolate as sintrp
@@ -75,6 +78,10 @@ def __init__(self, config, trode, vInd, pInd,
         self.phi_lyte = self.portInLyte.phi_lyte
         self.c_lyte = self.portInLyte.c_lyte
         self.phi_m = self.portInBulk.phi_m
+        if self.config[f"simInterface_{trode}"]:
+            self.portOutParticle = ports.portFromParticle(
+                "portOutParticle", dae.eOutletPort, self,
+                "Outlet port from particle")
 
     def get_trode_param(self, item):
         """
@@ -109,7 +116,7 @@ def DeclareEquations(self):
         # Figure out mu_O, mu of the oxidized state
         mu_O, act_lyte = calc_mu_O(
             self.c_lyte(), self.phi_lyte(), self.phi_m(), T,
-            self.config["elyteModelType"])
+            self.config, self.trode)
 
         # Define average filling fractions in particle
         eq1 = self.CreateEquation("c1bar")
@@ -120,28 +127,33 @@ 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() - .5*(self.c1bar() + self.c2bar())
+        eq.Residual = self.cbar() - (0.5*self.c1bar() + 0.5*self.c2bar())
 
         # Define average rate of filling of particle
         eq = self.CreateEquation("dcbardt")
         eq.Residual = self.dcbardt()
         for k in range(N):
-            eq.Residual -= .5*(self.c1.dt(k) + self.c2.dt(k)) * volfrac_vec[k]
+            eq.Residual -= (0.5*self.c1.dt(k) + 0.5*self.c2.dt(k)) * volfrac_vec[k]
 
         c1 = np.empty(N, dtype=object)
         c2 = np.empty(N, dtype=object)
         c1[:] = [self.c1(k) for k in range(N)]
         c2[:] = [self.c2(k) for k in range(N)]
-        if self.get_trode_param("type") in ["diffn2", "CHR2"]:
-            # Equations for 1D particles of 1 field varible
+        if self.get_trode_param("type") in ["diffn2", "CHR2", "ACR2"]:
+            # Equations for 1D particles of 2 field variables
             self.sld_dynamics_1D2var(c1, c2, mu_O, act_lyte, noises)
         elif self.get_trode_param("type") in ["homog2", "homog2_sdn"]:
-            # Equations for 0D particles of 1 field variables
+            # Equations for 0D particles of 2 field variables
             self.sld_dynamics_0D2var(c1, c2, mu_O, act_lyte, noises)
 
         for eq in self.Equations:
             eq.CheckUnitsConsistency = False
 
+        if self.config[f"simInterface_{self.trode}"]:
+            # Output dcbardt to interface region
+            eq = self.CreateEquation("particle_to_interface_dcbardt")
+            eq.Residual = self.portOutParticle.dcbardt() - self.dcbardt()
+
     def sld_dynamics_0D2var(self, c1, c2, muO, act_lyte, noises):
         T = self.config["T"]
         c1_surf = c1
@@ -191,9 +203,17 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
             c2_surf = c2[-1]
             mu1R_surf, act1R_surf = mu1R[-1], act1R[-1]
             mu2R_surf, act2R_surf = mu2R[-1], act2R[-1]
-        eta1 = calc_eta(mu1R_surf, muO)
-        eta2 = calc_eta(mu2R_surf, muO)
+            eta1 = calc_eta(mu1R_surf, muO)
+            eta2 = calc_eta(mu2R_surf, muO)
         if self.get_trode_param("type") in ["ACR2"]:
+            c1_surf = c1
+            c2_surf = c2
+            (mu1R, mu2R), (act1R, act2R) = calc_muR(
+                (c1, c2), (self.c1bar(), self.c2bar()), self.config, self.trode, self.ind)
+            mu1R_surf, act1R_surf = mu1R, act1R
+            mu2R_surf, act2R_surf = mu2R, act2R
+            eta1 = calc_eta(mu1R, muO)
+            eta2 = calc_eta(mu2R, muO)
             eta1_eff = np.array([eta1[i]
                                  + self.Rxn1(i)*self.get_trode_param("Rfilm") for i in range(N)])
             eta2_eff = np.array([eta2[i]
@@ -222,7 +242,10 @@ def sld_dynamics_1D2var(self, c1, c2, muO, act_lyte, noises):
             eq2.Residual = self.Rxn2() - Rxn2
 
         # Get solid particle fluxes (if any) and RHS
-        if self.get_trode_param("type") in ["diffn2", "CHR2"]:
+        if self.get_trode_param("type") in ["ACR2"]:
+            RHS1 = 0.5*np.array([self.get_trode_param("delta_L")*self.Rxn1(i) for i in range(N)])
+            RHS2 = 0.5*np.array([self.get_trode_param("delta_L")*self.Rxn2(i) for i in range(N)])
+        elif self.get_trode_param("type") in ["diffn2", "CHR2"]:
             # Positive reaction (reduction, intercalation) is negative
             # flux of Li at the surface.
             Flux1_bc = -0.5 * self.Rxn1()
@@ -276,16 +299,24 @@ def __init__(self, config, trode, vInd, pInd,
         # Domain
         self.Dmn = dae.daeDomain("discretizationDomain", self, dae.unit(),
                                  "discretization domain")
-
         # Variables
         self.c = dae.daeVariable("c", mole_frac_t, self,
                                  "Concentration in active particle",
                                  [self.Dmn])
+
+        # Creation of the ghost points to assit BC
+
+        if self.get_trode_param("type") in ["ACR_Diff"]:
+            self.c_left_GP = dae.daeVariable("c_left", mole_frac_t, self,
+                                             "Concentration on the left side of the particle")
+            self.c_right_GP = dae.daeVariable("c_right", mole_frac_t, self,
+                                              "Concentration on the right side of the particle")
+
         self.cbar = dae.daeVariable(
             "cbar", mole_frac_t, self,
             "Average concentration in active particle")
         self.dcbardt = dae.daeVariable("dcbardt", dae.no_t, self, "Rate of particle filling")
-        if config[trode, "type"] not in ["ACR"]:
+        if config[trode, "type"] not in ["ACR", "ACR_Diff"]:
             self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction")
         else:
             self.Rxn = dae.daeVariable("Rxn", dae.no_t, self, "Rate of reaction", [self.Dmn])
@@ -306,6 +337,10 @@ def __init__(self, config, trode, vInd, pInd,
         self.phi_lyte = self.portInLyte.phi_lyte
         self.c_lyte = self.portInLyte.c_lyte
         self.phi_m = self.portInBulk.phi_m
+        if config[f"simInterface_{trode}"]:
+            self.portOutParticle = ports.portFromParticle(
+                "portOutParticle", dae.eOutletPort, self,
+                "Outlet port from particle")
 
     def get_trode_param(self, item):
         """
@@ -335,7 +370,7 @@ def DeclareEquations(self):
 
         # Figure out mu_O, mu of the oxidized state
         mu_O, act_lyte = calc_mu_O(self.c_lyte(), self.phi_lyte(), self.phi_m(), T,
-                                   self.config["elyteModelType"])
+                                   self.config, self.trode)
 
         # Define average filling fraction in particle
         eq = self.CreateEquation("cbar")
@@ -350,8 +385,12 @@ def DeclareEquations(self):
             eq.Residual -= self.c.dt(k) * volfrac_vec[k]
 
         c = np.empty(N, dtype=object)
-        c[:] = [self.c(k) for k in range(N)]
-        if self.get_trode_param("type") in ["ACR", "diffn", "CHR"]:
+        if self.get_trode_param("type") in ["ACR_Diff"]:
+            ctmp = utils.get_var_vec(self.c,N)
+            c = np.hstack((self.c_left_GP(),ctmp,self.c_right_GP()))
+        else:
+            c[:] = [self.c(k) for k in range(N)]
+        if self.get_trode_param("type") in ["ACR", "ACR_Diff", "diffn", "CHR"]:
             # Equations for 1D particles of 1 field varible
             self.sld_dynamics_1D1var(c, mu_O, act_lyte, self.noise)
         elif self.get_trode_param("type") in ["homog", "homog_sdn"]:
@@ -361,6 +400,11 @@ def DeclareEquations(self):
         for eq in self.Equations:
             eq.CheckUnitsConsistency = False
 
+        if self.config[f"simInterface_{self.trode}"]:
+            # Output dcbardt to interface region
+            eq = self.CreateEquation("particle_to_interface_dcbardt")
+            eq.Residual = self.portOutParticle.dcbardt() - self.dcbardt()
+
     def sld_dynamics_0D1var(self, c, muO, act_lyte, noise):
         T = self.config["T"]
         c_surf = c
@@ -389,8 +433,20 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise):
         dr, edges = geo.get_dr_edges(self.get_trode_param('shape'), N)
 
         # Get solid particle chemical potential, overpotential, reaction rate
-        if self.get_trode_param("type") in ["ACR"]:
+        if self.get_trode_param("type") in ["ACR", "ACR_Diff"]:
             c_surf = c
+            # surface diffusion in the ACR C3 model
+            if self.get_trode_param("type") in ["ACR_Diff"]:
+                dx = 1/np.size(c)
+                beta_s = self.get_trode_param("beta_s")
+                eqL = self.CreateEquation("leftBC")
+                eqL.Residual = (c_surf[0] - c_surf[1] +
+                                - dx*beta_s*(c_surf[1]+0.008)*(1-c_surf[1]-0.008))
+                eqR = self.CreateEquation("rightBC")
+                eqR.Residual = (c_surf[-1] - c_surf[-2] +
+                                - dx*beta_s*(c_surf[-2]+0.008)*(1-c_surf[-2]-0.008))
+
+        if self.get_trode_param("type") in ["ACR", "ACR_Diff"]:
             muR_surf, actR_surf = calc_muR(
                 c_surf, self.cbar(), self.config, self.trode, self.ind)
         elif self.get_trode_param("type") in ["diffn", "CHR"]:
@@ -401,17 +457,27 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise):
                 actR_surf = None
             else:
                 actR_surf = actR[-1]
-        eta = calc_eta(muR_surf, muO)
-        if self.get_trode_param("type") in ["ACR"]:
+        # surface diffusion in the ACR C3 model
+        if self.get_trode_param("type") in ["ACR_Diff"]:
+            eta = calc_eta(muR_surf[1:-1], muO)
+        else:
+            eta = calc_eta(muR_surf, muO)
+        if self.get_trode_param("type") in ["ACR", "ACR_Diff"]:
             eta_eff = np.array([eta[i] + self.Rxn(i)*self.get_trode_param("Rfilm")
                                 for i in range(N)])
         else:
             eta_eff = eta + self.Rxn()*self.get_trode_param("Rfilm")
-        Rxn = self.calc_rxn_rate(
-            eta_eff, c_surf, self.c_lyte(), self.get_trode_param("k0"),
-            self.get_trode_param("E_A"), T, actR_surf, act_lyte,
-            self.get_trode_param("lambda"), self.get_trode_param("alpha"))
-        if self.get_trode_param("type") in ["ACR"]:
+        if self.get_trode_param("type") in ["ACR_Diff"]:
+            Rxn = self.calc_rxn_rate(
+                eta_eff, c_surf[1:-1], self.c_lyte(), self.get_trode_param("k0"),
+                self.get_trode_param("E_A"), T, actR_surf[1:-1], act_lyte,
+                self.get_trode_param("lambda"), self.get_trode_param("alpha"))
+        else:
+            Rxn = self.calc_rxn_rate(
+                eta_eff, c_surf, self.c_lyte(), self.get_trode_param("k0"),
+                self.get_trode_param("E_A"), T, actR_surf, act_lyte,
+                self.get_trode_param("lambda"), self.get_trode_param("alpha"))
+        if self.get_trode_param("type") in ["ACR", "ACR_Diff"]:
             for i in range(N):
                 eq = self.CreateEquation("Rxn_{i}".format(i=i))
                 eq.Residual = self.Rxn(i) - Rxn[i]
@@ -420,8 +486,34 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise):
             eq.Residual = self.Rxn() - Rxn
 
         # Get solid particle fluxes (if any) and RHS
-        if self.get_trode_param("type") in ["ACR"]:
-            RHS = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(N)])
+        if self.get_trode_param("type") in ["ACR", "ACR_Diff"]:
+            # For ACR model localized contact loss.
+            if self.config['localized_losses']:
+                gamma_con = self.get_trode_param('gamma_con')
+                if gamma_con == 1:
+                    RHS = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(N)])
+                else:
+                    N_cont = int((gamma_con)*N)-12
+                    RHS = np.zeros(N, dtype=object)
+                    position = int(np.random.uniform()*N)
+                    # random position of contact + 6 points on the sides to facilitate wetting
+                    if 6+position+N_cont < N-6:
+                        for i in range(0,6):
+                            RHS[i] = self.get_trode_param("delta_L")*self.Rxn(i)
+                        for i in range(N-6,N):
+                            RHS[i] = self.get_trode_param("delta_L")*self.Rxn(i)
+                        for i in range(position+6,position+N_cont,1):
+                            RHS[i] = self.get_trode_param("delta_L")*self.Rxn(i)
+                    else:
+                        for i in range(0,6):
+                            RHS[i] = self.get_trode_param("delta_L")*self.Rxn(i)
+                        for i in range(position, N):
+                            RHS[i] = self.get_trode_param("delta_L")*self.Rxn(i)
+                        for i in range(6,N_cont-(N-position)):
+                            RHS[i] = self.get_trode_param("delta_L")*self.Rxn(i)
+            else:
+                RHS = np.array([self.get_trode_param("delta_L")*self.Rxn(i) for i in range(N)])
+
         elif self.get_trode_param("type") in ["diffn", "CHR"]:
             # Positive reaction (reduction, intercalation) is negative
             # flux of Li at the surface.
@@ -445,9 +537,25 @@ def sld_dynamics_1D1var(self, c, muO, act_lyte, noise):
         dcdt_vec = np.empty(N, dtype=object)
         dcdt_vec[0:N] = [self.c.dt(k) for k in range(N)]
         LHS_vec = MX(Mmat, dcdt_vec)
+        if self.get_trode_param("type") in ["ACR_Diff"]:
+            # surface diffusion in the ACR C3 model
+            surf_diff_vec = calc_surf_diff(c_surf, muR_surf, self.get_trode_param("D"))
         for k in range(N):
             eq = self.CreateEquation("dcsdt_discr{k}".format(k=k))
-            eq.Residual = LHS_vec[k] - RHS[k]
+            if self.get_trode_param("type") in ["ACR_Diff"]:
+                eq.Residual = LHS_vec[k] - RHS[k] - surf_diff_vec[k]
+            else:
+                eq.Residual = LHS_vec[k] - RHS[k]
+
+
+# surface diffusion in the ACR C3 model
+def calc_surf_diff(c_surf, muR_surf, D):
+    N_2 = np.size(c_surf)
+    dxs = 1./N_2
+    c_surf_long = c_surf
+    c_surf_short = (c_surf_long[0:-1] + c_surf_long[1:])/2
+    surf_diff = D*(np.diff(c_surf_short*(1-c_surf_short)*np.diff(muR_surf)))/(dxs**2)
+    return surf_diff
 
 
 def calc_eta(muR, muO):
@@ -524,13 +632,23 @@ def calc_flux_CHR2(c1, c2, mu1_R, mu2_R, D, Dfunc, E_D, Flux1_bc, Flux2_bc, dr,
     return Flux1_vec, Flux2_vec
 
 
-def calc_mu_O(c_lyte, phi_lyte, phi_sld, T, elyteModelType):
+def calc_mu_O(c_lyte, phi_lyte, phi_sld, T, config, trode):
+    elyteModelType = config["elyteModelType"]
+
+    if config[f"simInterface_{trode}"]:
+        elyteModelType = config["interfaceModelType"]
+
     if elyteModelType == "SM":
         mu_lyte = phi_lyte
         act_lyte = c_lyte
     elif elyteModelType == "dilute":
         act_lyte = c_lyte
         mu_lyte = T*np.log(act_lyte) + phi_lyte
+    elif elyteModelType == "solid":
+        a_slyte = config['a_slyte']
+        cmax = config['cmax']
+        act_lyte = (c_lyte / cmax) / (1 - c_lyte / cmax)*np.exp(a_slyte*(1 - 2*c_lyte))
+        mu_lyte = phi_lyte
     mu_O = mu_lyte - phi_sld
     return mu_O, act_lyte
 
diff --git a/mpet/mod_interface.py b/mpet/mod_interface.py
new file mode 100644
index 00000000..a637e127
--- /dev/null
+++ b/mpet/mod_interface.py
@@ -0,0 +1,228 @@
+import numpy as np
+
+import daetools.pyDAE as dae
+from mpet import ports, utils
+import mpet.geometry as geom
+from mpet.daeVariableTypes import mole_frac_t, elec_pot_t
+
+
+"""
+Model for the interface between the electrolyte and active particles
+"""
+
+
+class InterfaceRegion(dae.daeModel):
+    def __init__(self, config, Name, Parent=None, Description="",
+                 cell=None, particle=None, vInd=None, pInd=None,
+                 trode=None):
+        super().__init__(Name, Parent, Description)
+        self.config = config
+
+        # Domain
+        self.Dmn = dae.daeDomain("discretizationDomain", self, dae.unit(),
+                                 "discretization domain")
+
+        # Variables
+        self.c = dae.daeVariable("c", mole_frac_t, self,
+                                 "Concentration in interface",
+                                 [self.Dmn])
+
+        self.phi = dae.daeVariable("phi", elec_pot_t, self,
+                                   "Electrical potential in interface",
+                                   [self.Dmn])
+
+        # Ports
+        self.portInLyte = ports.portFromElyte(
+            "portInLyte", dae.eInletPort, self,
+            "Inlet port from electrolyte")
+
+        self.portInParticle = ports.portFromParticle(
+            "portInParticle", dae.eInletPort, self,
+            "Inlet port from particle")
+
+        # Note: using portFromElyte here because the port from the
+        # interface to the particles transfers the same variables as
+        # the port from the elyte to the interface, hence the same
+        # class can be used
+        self.portOutInterfaceParticle = ports.portFromElyte(
+            "portOutInterfaceParticle", dae.eOutletPort, self,
+            "Port from interface to particles")
+
+        self.portOutInterfaceElyte = ports.portFromInterface(
+            "portOutInterfaceElyte", dae.eOutletPort, self,
+            "Port from interface to elyte")
+
+        # Particle
+        self.particle = particle
+
+        # Cell
+        self.cell = cell
+
+        # Volume and particle indices
+        self.vInd = vInd
+        self.pInd = pInd
+        self.trode = trode
+
+    def DeclareEquations(self):
+        super().DeclareEquations()
+        config = self.config
+        Nvol = config["Nvol_i"]
+
+        disc = geom.get_interface_disc(Nvol, config["L_i"],
+                                       config["poros_i"], config["BruggExp_i"])
+        cvec = utils.get_var_vec(self.c, Nvol)
+        dcdtvec = utils.get_var_vec(self.c, Nvol, dt=True)
+        phivec = utils.get_var_vec(self.phi, Nvol)
+
+        # Apply concentration and potential boundary conditions
+        # Elyte value on the left and no-gradients on the right
+        T = self.config["T"]
+
+        # Concentration: if elyte solid not continuous over
+        # interface if liquid continous over interface.
+        if config["interfaceModelType"] == "solid" or config["elyteModelType"] == "solid":
+            ctmp = np.hstack((cvec[0], cvec, cvec[-1]))
+        else:
+            ctmp = np.hstack((self.portInLyte.c_lyte(), cvec, cvec[-1]))
+
+        # Electrical potential
+        if config["elyteModelType"] == "dilute":
+            phitmp = np.hstack((self.portInLyte.phi_lyte()
+                                + T * np.log(self.portInLyte.c_lyte()), phivec, phivec[-1]))
+        elif config["interfaceModelType"] == "dilute":
+            phitmp = np.hstack((self.portInLyte.phi_lyte()
+                                - T * np.log(cvec[0]), phivec, phivec[-1]))
+        else:
+            phitmp = np.hstack((self.portInLyte.phi_lyte(), phivec, phivec[-1]))
+
+        Nm_edges, i_edges = get_interface_internal_fluxes(ctmp, phitmp, disc, config)
+
+        # The reaction rate per volume (Rvp) is normalized to the total length of the electrode.
+        dlc = config["L"][self.trode]/config["Nvol"][self.trode]
+        disc["dxvec"][:] = 1
+
+        dvgNm = np.diff(Nm_edges) / disc["dxvec"]
+        dvgi = np.diff(i_edges) / disc["dxvec"]
+
+        for vInd in range(Nvol):
+            # Mass Conservation (done with the anion, although "c" is neutral salt conc)
+            eq = self.CreateEquation("interface_mass_cons_vol{vInd}".format(vInd=vInd))
+            eq.Residual = disc["porosvec"][vInd]*dcdtvec[vInd] + (1./config["num"])*dvgNm[vInd]
+
+            # Charge Conservation
+            eq = self.CreateEquation("interface_charge_cons_vol{vInd}".format(vInd=vInd))
+            eq.Residual = -dvgi[vInd]
+            # Reaction out interface from last volume
+            if vInd == Nvol - 1:
+                # The volume of this particular particle
+                Vj = config["psd_vol_FracVol"][self.trode][self.vInd,self.pInd]
+                eq.Residual += dlc * config["zp"] * -(config["beta"][self.trode]
+                                                      * (1-config["poros"][self.trode])
+                                                      * config["P_L"][self.trode] * Vj
+                                                      * self.portInParticle.dcbardt())
+
+            # Reaction entering the interface
+            if vInd == 0:
+                # The volume of this particular particle
+                Vj = config["psd_vol_FracVol"][self.trode][self.vInd,self.pInd]
+                eq.Residual -= dlc * config["zp"] * -(config["beta"][self.trode]
+                                                      * (1-config["poros"][self.trode])
+                                                      * config["P_L"][self.trode] * Vj
+                                                      * self.portInParticle.dcbardt())
+
+        # last grid point of interface is output to particle
+        eq = self.CreateEquation("c_interface_to_particle")
+        eq.Residual = self.portOutInterfaceParticle.c_lyte() - ctmp[-1]
+
+        eq = self.CreateEquation("phi_interface_to_particle")
+        eq.Residual = self.portOutInterfaceParticle.phi_lyte() - phitmp[-1]
+
+        eq = self.CreateEquation("Nm0_interface_to_elyte")
+        eq.Residual = self.portOutInterfaceElyte.Nm0() - Nm_edges[0]
+
+        eq = self.CreateEquation("i0_interface_to_elyte")
+        # eq.Residual = self.portOutInterfaceElyte.i0() - i_edges[1]
+        eq.Residual = self.portOutInterfaceElyte.i0() - i_edges[1] / dlc
+
+        for eq in self.Equations:
+            eq.CheckUnitsConsistency = False
+
+
+def get_interface_internal_fluxes(c, phi, disc, config):
+    zp, zm, nup, num = config["zp"], config["zm"], config["nup"], config["num"]
+    nu = nup + num
+    T = config["T"]
+    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, wt)
+
+    if config["interfaceModelType"] == "dilute":
+        # 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 * config["Dp_i"]
+        Dm = eps_o_tau_edges * config["Dm_i"]
+#        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)/dxd1
+                            - Dm/T*zm*c_edges_int*np.diff(phi)/dxd1)
+        i_edges_int = (-((nup*zp*Dp + num*zm*Dm)*np.diff(c)/dxd1)
+                       - (nup*zp**2*Dp + num*zm**2*Dm)/T*c_edges_int*np.diff(phi)/dxd1)
+#        i_edges_int = zp*Np_edges_int + zm*Nm_edges_int
+    elif config["interfaceModelType"] == "SM":
+        SMset = config["SMset"]
+        elyte_function = utils.import_function(config["SMset_filename"], SMset,
+                                               mpet_module=f"mpet.electrolyte.{SMset}")
+        D_fs, sigma_fs, thermFac, tp0 = elyte_function()[:-1]
+
+        # Get diffusivity and conductivity at cell edges using weighted harmonic mean
+        D_edges = utils.weighted_harmonic_mean(eps_o_tau*D_fs(c, T), wt)
+        sigma_edges = utils.weighted_harmonic_mean(eps_o_tau*sigma_fs(c, T), wt)
+
+        sp, n = config["sp"], config["n"]
+        # there is an error in the MPET paper, temperature dependence should be
+        # in sigma and not outside of sigma
+        i_edges_int = -sigma_edges * (
+            np.diff(phi)/dxd1
+            + nu*T*(sp/(n*nup)+tp0(c_edges_int, T)/(zp*nup))
+            * thermFac(c_edges_int, T)
+            * np.diff(np.log(c))/dxd1
+            )
+        Nm_edges_int = num*(-D_edges*np.diff(c)/dxd1
+                            + (1./(num*zm)*(1-tp0(c_edges_int, T))*i_edges_int))
+
+    elif config["interfaceModelType"] == "solid":
+        SMset = config["SMset"]
+        elyte_function = utils.import_function(config["SMset_filename"], SMset,
+                                               mpet_module=f"mpet.electrolyte.{SMset}")
+        D_fs, sigma_fs, thermFac, tp0 = elyte_function()[:-1]
+
+        a_slyte = config["a_slyte"]
+        tp0 = 0.99999
+
+        c_edges_int_norm = c_edges_int / config["cmax_i"]
+
+        # Get diffusivity at cell edges using weighted harmonic mean
+        # D_edges = utils.weighted_harmonic_mean(eps_o_tau * D_fs(c_lyte), wt)
+        eps_o_tau_edges = utils.weighted_linear_mean(eps_o_tau, wt)
+        # sp, n = ndD["sp"], ndD["n_refTrode"]
+        # D_fs is specified in solid_elyte_func in props_elyte.py
+        Dp = eps_o_tau_edges * config["Dp_i"]
+        Dm = (zp * Dp - zp * Dp * tp0) / (tp0 * zm)
+
+        Dp0 = Dp / (1-c_edges_int_norm)  # should be c0/cmax
+
+        Dchemp = Dp0 * (1 - 2 * a_slyte * c_edges_int_norm + 2 * a_slyte * c_edges_int_norm**2)
+        Dchemm = Dm
+
+        Damb = (zp * Dp * Dchemm + zm * Dm * Dchemp) / (zp * Dp - zm * Dm)
+
+        i_edges_int = (-((nup*zp*Dchemp + num*zm*Dchemm)*np.diff(c)/dxd1)
+                       - (nup * zp ** 2 * Dp0 * (1 - c_edges_int_norm) + num * zm ** 2 * Dm) / T
+                       * c_edges_int * np.diff(phi) / dxd1)
+
+        Nm_edges_int = num * (-Damb * np.diff(c) / dxd1
+                              + (1. / (num * zm) * (1 - tp0) * i_edges_int))
+    return Nm_edges_int, i_edges_int
diff --git a/mpet/plot/outmat2txt.py b/mpet/plot/outmat2txt.py
index d802466b..e97d2807 100644
--- a/mpet/plot/outmat2txt.py
+++ b/mpet/plot/outmat2txt.py
@@ -79,7 +79,7 @@ def main(indir, genData=True, discData=True, elyteData=True,
          csldData=True, cbarData=True, bulkpData=True):
     config = plot_data.show_data(
         indir, plot_type="params", print_flag=False, save_flag=False,
-        data_only=True)
+        data_only=True, color_changes=None, smooth_type=None)
     trodes = config["trodes"]
     CrateCurr = config["1C_current_density"]  # A/m^2
     psd_len_c = config["psd_len"]["c"]
@@ -93,25 +93,25 @@ def get_trode_str(tr):
         Nv_a, Np_a = psd_len_a.shape
     tVec, vVec = plot_data.show_data(
         indir, plot_type="vt", print_flag=False, save_flag=False,
-        data_only=True)
+        data_only=True, color_changes=None, smooth_type=None)
     ntimes = len(tVec)
 
     if genData:
         ffVec_c = plot_data.show_data(
             indir, plot_type="soc_c", print_flag=False,
-            save_flag=False, data_only=True)[1]
+            save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
         if "a" in trodes:
             ffVec_a = plot_data.show_data(
                 indir, plot_type="soc_a", print_flag=False,
-                save_flag=False, data_only=True)[1]
+                save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
         else:
             ffVec_a = np.ones(len(tVec))
         currVec = plot_data.show_data(
             indir, plot_type="curr", print_flag=False,
-            save_flag=False, data_only=True)[1]
+            save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
         powerVec = plot_data.show_data(
             indir, plot_type="power", print_flag=False,
-            save_flag=False, data_only=True)[1]
+            save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
         genMat = np.zeros((ntimes, 7))
         genMat[:,0] = tVec
         genMat[:,1] = ffVec_a
@@ -126,7 +126,7 @@ def get_trode_str(tr):
     if discData:
         cellCentersVec, facesVec = plot_data.show_data(
             indir, plot_type="discData", print_flag=False,
-            save_flag=False, data_only=True)
+            save_flag=False, data_only=True, color_changes=None, smooth_type=None)
         with open(os.path.join(indir, "discData.txt"), "w") as fo:
             print(discCCbattery, file=fo)
             print(",".join(map(str, cellCentersVec)), file=fo)
@@ -164,15 +164,16 @@ def get_trode_str(tr):
         # so we'll get a KeyError in attempting to "plot" the electrolyte current density.
         try:
             plot_data.show_data(
-                indir, plot_type="elytei", print_flag=False, save_flag=False, data_only=True)
+                indir, plot_type="elytei", print_flag=False, save_flag=False, data_only=True,
+                color_changes=None, smooth_type=None)
         except KeyError:
             valid_current = False
         elytecMat = plot_data.show_data(
             indir, plot_type="elytec", print_flag=False,
-            save_flag=False, data_only=True)[1]
+            save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
         elytepMat = plot_data.show_data(
             indir, plot_type="elytep", print_flag=False,
-            save_flag=False, data_only=True)[1]
+            save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
         np.savetxt(os.path.join(indir, "elyteConcData.txt"),
                    elytecMat, delimiter=dlm, header=elytecHdr)
         np.savetxt(os.path.join(indir, "elytePotData.txt"),
@@ -180,10 +181,10 @@ def get_trode_str(tr):
         if valid_current:
             elyteiMat = plot_data.show_data(
                 indir, plot_type="elytei", print_flag=False,
-                save_flag=False, data_only=True)[1]
+                save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
             elytediviMat = plot_data.show_data(
                 indir, plot_type="elytedivi", print_flag=False,
-                save_flag=False, data_only=True)[1]
+                save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
             np.savetxt(os.path.join(indir, "elyteCurrDensData.txt"),
                        elyteiMat, delimiter=dlm, header=elyteiHdr)
             np.savetxt(os.path.join(indir, "elyteDivCurrDensData.txt"),
@@ -229,7 +230,7 @@ def get_trode_str(tr):
     if cbarData:
         cbarDict = plot_data.show_data(
             indir, plot_type="cbar_full", print_flag=False,
-            save_flag=False, data_only=True)
+            save_flag=False, data_only=True, color_changes='discrete', smooth_type=None)
         for tr in trodes:
             Trode = get_trode_str(tr)
             fname = "cbar{l}Data.txt".format(l=Trode)
@@ -250,13 +251,13 @@ def get_trode_str(tr):
         if "a" in trodes:
             bulkp_aData = plot_data.show_data(
                 indir, plot_type="bulkp_a", print_flag=False,
-                save_flag=False, data_only=True)[1]
+                save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
             fname = fnameBulkpBase.format(l="Anode")
             np.savetxt(os.path.join(indir, fname), bulkp_aData,
                        delimiter=dlm, header=bulkpHdr)
         bulkp_cData = plot_data.show_data(
             indir, plot_type="bulkp_c", print_flag=False,
-            save_flag=False, data_only=True)[1]
+            save_flag=False, data_only=True, color_changes=None, smooth_type=None)[1]
         fname = fnameBulkpBase.format(l="Cathode")
         np.savetxt(os.path.join(indir, fname), bulkp_cData,
                    delimiter=dlm, header=bulkpHdr)
diff --git a/mpet/plot/plot_data.py b/mpet/plot/plot_data.py
index c9550bc4..c9c0e8a8 100644
--- a/mpet/plot/plot_data.py
+++ b/mpet/plot/plot_data.py
@@ -25,7 +25,8 @@
 # mpl.rcParams['text.usetex'] = True
 
 
-def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOut=None, tOut=None):
+def show_data(indir, plot_type, print_flag, save_flag, data_only, color_changes, smooth_type,
+              vOut=None, pOut=None, tOut=None):
     pfx = 'mpet.'
     sStr = "_"
     ttl_fmt = "% = {perc:2.1f}"
@@ -69,7 +70,7 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
     dxvec = np.array(Nvol["c"] * [dxc])
     porosvec = np.array(Nvol["c"] * [config["poros"]["c"]])
     cellsvec = dxc*np.arange(Nvol["c"]) + dxc/2.
-    if config["have_separator"]:
+    if config["Nvol"]["s"]:
         dxs = config["L"]["s"]/Nvol["s"]
         dxvec_s = np.array(Nvol["s"] * [dxs])
         dxvec = np.hstack((dxvec_s, dxvec))
@@ -113,7 +114,7 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
 #            print "Actual psd_stddev [nm]:", np.std(psd_len[l])
         print("Cell structure:")
         print(("porous anode | " if "a" in config["trodes"] else "flat anode | ")
-              + ("sep | " if config["have_separator"] else "") + "porous cathode")
+              + ("sep | " if config["Nvol"]["s"] else "") + "porous cathode")
         if "a" in config["trodes"]:
             print("capacity ratio cathode:anode, 'z':", config["z"])
         for trode in trodes:
@@ -136,7 +137,7 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
         print("Specified psd_stddev, c [{unit}]:".format(unit=Lunit),
               np.array(config['stddev']["c"])*Lfac)
         print("ndim B_c:", config["c", "B"])
-        if config["have_separator"]:
+        if config["Nvol"]["s"]:
             print("Nvol_s:", Nvol["s"])
         print("Nvol_c:", Nvol["c"])
         if 'a' in config["trodes"]:
@@ -236,6 +237,8 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
                 ax[pInd,vInd].xaxis.set_major_locator(plt.NullLocator())
                 datay = utils.get_dict_key(data, sol_str, squeeze=False)[:,-1]
                 line, = ax[pInd,vInd].plot(times, datay)
+        if save_flag:
+            fig.savefig("mpet_surf.pdf", bbox_inches="tight")
         return fig, ax
 
     # Plot SoC profile
@@ -280,6 +283,8 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
             return times*td, cavg
         np.set_printoptions(precision=8)
         ax.plot(times*td, cavg)
+        if save_flag:
+            fig.savefig("mpet_elytecons.pdf", bbox_inches="tight")
         return fig, ax
 
     # Plot current profile
@@ -327,7 +332,7 @@ def show_data(indir, plot_type, print_flag, save_flag, data_only, vOut=None, pOu
         datay_p = utils.get_dict_key(data, p_cath, squeeze=False)
         L_c = config['L']["c"] * config['L_ref'] * Lfac
         Ltot = L_c
-        if config["have_separator"]:
+        if config["Nvol"]["s"]:
             datay_s_c = utils.get_dict_key(data, c_sep, squeeze=False)
             datay_s_p = utils.get_dict_key(data, p_sep, squeeze=False)
             datay_c = np.hstack((datay_s_c, datay_c))
@@ -577,7 +582,10 @@ def animate(tind):
                 ax[pInd,vInd].set_ylim(ylim)
                 ax[pInd,vInd].set_xlim((0, lens[pInd,vInd] * Lfac))
                 if plt_axlabels:
-                    ax[pInd, vInd].set_xlabel(r"$r$ [{Lunit}]".format(Lunit=Lunit))
+                    if config[trode, "type"] in ["ACR", "ACr_diff", "ACR2"]:
+                        ax[pInd, vInd].set_xlabel(r"$x$ [{Lunit}]".format(Lunit=Lunit))
+                    else:
+                        ax[pInd, vInd].set_xlabel(r"$r$ [{Lunit}]".format(Lunit=Lunit))
                     if plot_type[0] == "c":
                         ax[pInd, vInd].set_ylabel(r"$\widetilde{{c}}$")
                     elif plot_type[:2] == "mu":
@@ -665,12 +673,11 @@ def animate(tind):
         if data_only:
             return dataCbar
         # Set up colors.
-        # Define if you want smooth or discrete color changes
-        # Option: "smooth" or "discrete"
-        color_changes = "discrete"
-#        color_changes = "smooth"
+        # Uses either discrete or smooth colors
+        # Define if you want smooth or discrete color changes in plot settings (-c)
+        # Option: "discrete" or "smooth"
         # Discrete color changes:
-        if color_changes == "discrete":
+        if color_changes == 'discrete':
             # Make a discrete colormap that goes from green to yellow
             # to red instantaneously
             cdict = {
@@ -687,10 +694,11 @@ def animate(tind):
             cmap = mpl.colors.LinearSegmentedColormap(
                 "discrete", cdict)
         # Smooth colormap changes:
-        if color_changes == "smooth":
+        if color_changes == 'smooth':
             # generated with colormap.org
-            cmaps = np.load("colormaps_custom.npz")
-            cmap_data = cmaps["GnYlRd_3"]
+            cmap_location = os.path.dirname(os.path.abspath(__file__)) + r'\colormaps_custom.npz'
+            cmaps = np.load(cmap_location)
+            cmap_data = cmaps[smooth_type]
             cmap = mpl.colors.ListedColormap(cmap_data/255.)
 
         size_frac_min = 0.10
diff --git a/mpet/plot/plot_data_db.py b/mpet/plot/plot_data_db.py
new file mode 100644
index 00000000..14939d62
--- /dev/null
+++ b/mpet/plot/plot_data_db.py
@@ -0,0 +1,435 @@
+import pandas as pd
+import numpy as np
+import os
+import argparse
+from argparse import RawTextHelpFormatter
+
+import mpet.geometry as geom
+from mpet import mod_cell
+from mpet import utils
+from mpet.config import Config, constants
+from mpet.exceptions import UnknownParameterError
+
+
+###################################################################
+# Part 1
+# Import data from all subfolders in provided dataDir
+# Prepare all data (get everything in correct dataframe) for plotting
+# This is based on plot_data.py
+###################################################################
+def main():
+    desc = """ Dashboard that shows all plots and compares the resutls of different models."""
+    parser = argparse.ArgumentParser(description=desc, formatter_class=RawTextHelpFormatter)
+    parser.add_argument('-d', '--dataDir',
+                        help='Directory that contains subfolders with simulation output')
+    args = parser.parse_args()
+
+    dataFiles = [os.path.join(args.dataDir, f) for f in os.listdir(args.dataDir)
+                 if os.path.isdir(os.path.join(args.dataDir, f))]
+    dff = pd.DataFrame()
+    dff_c_sub = pd.DataFrame()
+    dff_cd_sub = pd.DataFrame()
+    dff_csld_sub = pd.DataFrame()
+    dff_bulkp = pd.DataFrame()
+    df_cbar = pd.DataFrame()
+
+    for indir in dataFiles:
+        print('Loading data from:', indir)
+        pfx = 'mpet.'
+        sStr = "_"
+        # Read in the simulation results and calcuations data
+        model = os.path.basename(indir)
+        dataFileName = "output_data"
+        dataFile = os.path.join(indir, dataFileName)
+        data = utils.open_data_file(dataFile)
+        try:
+            utils.get_dict_key(data, pfx + 'current')
+        except KeyError:
+            pfx = ''
+        try:
+            utils.get_dict_key(data, pfx + "partTrodecvol0part0" + sStr + "cbar")
+        except KeyError:
+            sStr = "."
+        # Read in the parameters used to define the simulation
+        config = Config.from_dicts(indir)
+        # simulated (porous) electrodes
+        trodes = config["trodes"]
+        # Pick out some useful calculated values
+        limtrode = config["limtrode"]
+        k = constants.k                      # Boltzmann constant, J/(K Li)
+        Tref = constants.T_ref               # Temp, K
+        e = constants.e                      # Charge of proton, C
+        F = constants.F                      # C/mol
+        c_ref = constants.c_ref
+        td = config["t_ref"]
+        Etheta = {"a": 0.}
+        cap = config[limtrode, "cap"]
+        for trode in trodes:
+            Etheta[trode] = -(k*Tref/e) * config[trode, "phiRef"]
+        Vstd = Etheta["c"] - Etheta["a"]
+        Nvol = config["Nvol"]
+        Npart = config["Npart"]
+        psd_len = config["psd_len"]
+        # Discretization (and associated porosity)
+        Lfac = 1e6
+        dxc = config["L"]["c"]/Nvol["c"]
+        dxvec = np.array(Nvol["c"] * [dxc])
+        porosvec = np.array(Nvol["c"] * [config["poros"]["c"]])
+        cellsvec = dxc*np.arange(Nvol["c"]) + dxc/2.
+        if config["have_separator"]:
+            dxs = config["L"]["s"]/Nvol["s"]
+            dxvec_s = np.array(Nvol["s"] * [dxs])
+            dxvec = np.hstack((dxvec_s, dxvec))
+            poros_s = np.array(Nvol["s"] * [config["poros"]["s"]])
+            porosvec = np.hstack((poros_s, porosvec))
+            cellsvec += config["L"]["s"] / config["L"]["c"]
+            cellsvec_s = dxs*np.arange(Nvol["s"]) + dxs/2.
+            cellsvec = np.hstack((cellsvec_s, cellsvec))
+        if "a" in trodes:
+            dxa = config["L"]["a"]/Nvol["a"]
+            dxvec_a = np.array(Nvol["a"] * [dxa])
+            dxvec = np.hstack((dxvec_a, dxvec))
+            poros_a = np.array(Nvol["a"] * [config["poros"]["a"]])
+            porosvec = np.hstack((poros_a, porosvec))
+            cellsvec += config["L"]["a"] / config["L"]["c"]
+            cellsvec_a = dxa*np.arange(Nvol["a"]) + dxa/2.
+            cellsvec = np.hstack((cellsvec_a, cellsvec))
+        cellsvec *= config["L_ref"] * Lfac
+        facesvec = np.insert(np.cumsum(dxvec), 0, 0.) * config["L_ref"] * Lfac
+        # Extract the reported simulation times
+        times = utils.get_dict_key(data, pfx + 'phi_applied_times')
+        numtimes = len(times)
+        tmin = np.min(times)
+        tmax = np.max(times)
+        # Voltage profile
+        timestd = times*td
+        voltage = (Vstd - (k*Tref/e)*utils.get_dict_key(data, pfx + 'phi_applied'))
+        # surface concentration
+        # soc profile
+        ffvec_c = utils.get_dict_key(data, pfx + 'ffrac_c')
+        if "a" in trodes:
+            ffvec_a = utils.get_dict_key(data, pfx + 'ffrac_a')
+            nparta = Npart["a"]
+            nvola = Nvol["a"]
+        else:
+            ffvec_a = 0
+            nparta = 0
+            nvola = 0
+        # Elytecons
+        # current
+        theoretical_1C_current = config[config['limtrode'], "cap"] / 3600.  # A/m^2
+        current = (utils.get_dict_key(data, pfx + 'current')
+                   * theoretical_1C_current / config['1C_current_density'] * config['curr_ref'])
+        # Power
+        current_p = utils.get_dict_key(data, pfx + 'current') * (3600/td) * (cap/3600)  # in A/m^2
+        voltage_p = (Vstd - (k*Tref/e)*utils.get_dict_key(data, pfx + 'phi_applied'))  # in V
+        power = np.multiply(current_p, voltage_p)
+
+        # Electrolyte concetration / potential
+        datax = cellsvec
+        c_sep, p_sep = pfx + 'c_lyte_s', pfx + 'phi_lyte_s'
+        c_anode, p_anode = pfx + 'c_lyte_a', pfx + 'phi_lyte_a'
+        c_cath, p_cath = pfx + 'c_lyte_c', pfx + 'phi_lyte_c'
+        datay_c = utils.get_dict_key(data, c_cath, squeeze=False)
+        datay_p = utils.get_dict_key(data, p_cath, squeeze=False)
+        L_c = config['L']["c"] * config['L_ref'] * Lfac
+        Ltot = L_c
+        if config["have_separator"]:
+            datay_s_c = utils.get_dict_key(data, c_sep, squeeze=False)
+            datay_s_p = utils.get_dict_key(data, p_sep, squeeze=False)
+            datay_c = np.hstack((datay_s_c, datay_c))
+            datay_p = np.hstack((datay_s_p, datay_p))
+            L_s = config['L']["s"] * config['L_ref'] * Lfac
+            Ltot += L_s
+        else:
+            L_s = 0
+        if "a" in trodes:
+            datay_a_c = utils.get_dict_key(data, c_anode, squeeze=False)
+            datay_a_p = utils.get_dict_key(data, p_anode, squeeze=False)
+            datay_c = np.hstack((datay_a_c, datay_c))
+            datay_p = np.hstack((datay_a_p, datay_p))
+            L_a = config['L']["a"] * config['L_ref'] * Lfac
+            Ltot += L_a
+        else:
+            L_a = 0
+        # elytec
+        datay_ce = datay_c * c_ref / 1000.
+        # elytep
+        datay_pe = datay_p*(k*Tref/e) - Vstd
+        i_edges = np.zeros((numtimes, len(facesvec)))
+        # elytei & elytedivi
+        try:
+            cGP_L = utils.get_dict_key(data, "c_lyteGP_L")
+            pGP_L = utils.get_dict_key(data, "phi_lyteGP_L")
+            cmat = np.hstack((cGP_L.reshape((-1,1)), datay_c, datay_c[:,-1].reshape((-1,1))))
+            pmat = np.hstack((pGP_L.reshape((-1,1)), datay_p, datay_p[:,-1].reshape((-1,1))))
+            disc = geom.get_elyte_disc(Nvol, config["L"], config["poros"], config["BruggExp"])
+            for tInd in range(numtimes):
+                i_edges[tInd, :] = mod_cell.get_lyte_internal_fluxes(
+                    cmat[tInd, :], pmat[tInd, :], disc, config)[1]
+            datay_cd = i_edges * (F*constants.c_ref*config["D_ref"]/config["L_ref"])
+            datay_d = np.diff(i_edges, axis=1) / disc["dxvec"]
+            datay_d *= (F*constants.c_ref*config["D_ref"]/config["L_ref"]**2)
+        except (UnknownParameterError, KeyError):
+            datay_cd = i_edges
+            datay_d = np.zeros((numtimes, len(cellsvec)))
+            datay_d *= (F*constants.c_ref*config["D_ref"]/config["L_ref"]**2)
+        # fraction
+        t_current = times
+        tfrac = (t_current - tmin)/(tmax - tmin) * 100
+        # elytecons
+        sep = pfx + 'c_lyte_s'
+        anode = pfx + 'c_lyte_a'
+        cath = pfx + 'c_lyte_c'
+        cvec = utils.get_dict_key(data, cath)
+        if config["have_separator"]:
+            cvec_s = utils.get_dict_key(data, sep)
+            cvec = np.hstack((cvec_s, cvec))
+        if "a" in trodes:
+            cvec_a = utils.get_dict_key(data, anode)
+            cvec = np.hstack((cvec_a, cvec))
+        try:
+            cavg = np.sum(porosvec*dxvec*cvec, axis=1)/np.sum(porosvec*dxvec)
+        except np.AxisError:
+            cavg = np.sum(porosvec*dxvec*cvec, axis=0)/np.sum(porosvec*dxvec)
+        # Get all data in dataframes
+        df = pd.DataFrame({
+            "Model": model,
+            "sStr": sStr,
+            "pfx": pfx,
+            "Config trode type": config[trode, "type"],
+            "Voltage (V)": voltage,
+            "Cathode Filling Fraction": ffvec_c,
+            "Anode Filling Fraction": ffvec_a,
+            "Time (s)": timestd,
+            "Npartc": Npart["c"],
+            "Nvolc": Nvol["c"],
+            "Nparta": nparta,
+            "Nvola": nvola,
+            "Current": current,
+            "Power": power,
+            "cavg": cavg
+        })
+        for trode in trodes:
+            partStr = "partTrode{trode}vol{{vInd}}part{{pInd}}".format(trode=trode) + sStr
+            lens = psd_len[trode]
+            size_fracs = 0.4*np.ones((Nvol[trode], Npart[trode]))
+            if np.max(lens) != np.min(lens):
+                size_fracs = (lens - np.min(lens))/(np.max(lens) - np.min(lens))
+            size_frac_min = 0.2
+            sizes = (size_fracs*(1-size_frac_min) + size_frac_min) / Nvol[trode]
+            for pInd in range(Npart[trode]):
+                for vInd in range(Nvol[trode]):
+                    # for c data and cbar data
+                    if config[trode, "type"] in constants.one_var_types:
+                        str_base = (pfx + partStr + "c")
+                        sol_str = str_base.format(pInd=pInd, vInd=vInd)
+                        sol_str_data = utils.get_dict_key(data, sol_str, squeeze=False)[:,-1]
+                        str_cbar_base = pfx + partStr + "cbar"
+                        sol_cbar_str = str_cbar_base.format(pInd=pInd, vInd=vInd)
+                        sol_cbar_str_data = utils.get_dict_key(data, sol_cbar_str)
+                        df = pd.concat((df, pd.DataFrame({sol_str: sol_str_data,
+                                                          sol_cbar_str: sol_cbar_str_data})),
+                                       axis=1)
+                        # for cbar movie
+                        df_c = pd.DataFrame({
+                            'Model': model,
+                            "Config trode type": 1,
+                            'Time': np.round(timestd),
+                            'Cbar': sol_cbar_str_data,
+                            'r': pInd,
+                            'c': vInd,
+                            'rc': str(pInd)+str(vInd),
+                            'Relative size': sizes[vInd,pInd],
+                            'Trode': trode
+                            })
+                    elif config[trode, "type"] in constants.two_var_types:
+                        str1_base = (pfx + partStr + "c1")
+                        str2_base = (pfx + partStr + "c2")
+                        sol1_str = str1_base.format(pInd=pInd, vInd=vInd)
+                        sol2_str = str2_base.format(pInd=pInd, vInd=vInd)
+                        sol1_str_data = utils.get_dict_key(data, sol1_str, squeeze=False)[:,-1]
+                        sol2_str_data = utils.get_dict_key(data, sol2_str, squeeze=False)[:,-1]
+                        str1_cbar_base = pfx + partStr + "c1bar"
+                        str2_cbar_base = pfx + partStr + "c2bar"
+                        sol1_cbar_str = str1_cbar_base.format(pInd=pInd, vInd=vInd)
+                        sol2_cbar_str = str2_cbar_base.format(pInd=pInd, vInd=vInd)
+                        sol1_cbar_str_data = utils.get_dict_key(data, sol1_cbar_str)
+                        sol2_cbar_str_data = utils.get_dict_key(data, sol2_cbar_str)
+                        df = pd.concat((df, pd.DataFrame({sol1_str: sol1_str_data,
+                                                          sol2_str: sol2_str_data,
+                                                          sol1_cbar_str: sol1_cbar_str_data,
+                                                          sol2_cbar_str: sol2_cbar_str_data})),
+                                       axis=1)
+                        df_c = pd.DataFrame({
+                            'Model': model,
+                            "Config trode type": 2,
+                            'Time': np.round(timestd),
+                            'Cbar1': sol1_cbar_str_data,
+                            'Cbar2': sol2_cbar_str_data,
+                            'r': pInd,
+                            'c': vInd,
+                            'rc': str(pInd)+str(vInd),
+                            'Relative size': sizes[vInd,pInd],
+                            'Trode': trode
+                            })
+                    df_cbar = pd.concat([df_cbar, df_c])
+        dff = pd.concat([dff, df], ignore_index=True)
+        # build dataframe for plots electrolyte concentration or potential
+        # and for csld subplot animation (time, pind, vind, y)
+        dff_c = pd.DataFrame({"Model": model,
+                              "Time fraction": np.round(np.repeat(tfrac, np.shape(datay_ce)[1])),
+                              "fraction orig": np.repeat(tfrac, np.shape(datay_ce)[1]),
+                              "cellsvec": np.tile(cellsvec, np.shape(datay_ce)[0]),
+                              "Concentration electrolyte": datay_ce.flatten(),
+                              "Potential electrolyte": datay_pe.flatten(),
+                              "Divergence electrolyte curr dens": datay_d.flatten()
+                              })
+        dff_cd = pd.DataFrame({"Model": model,
+                               "Time fraction": np.round(np.repeat(tfrac, np.shape(datay_cd)[1])),
+                               "fraction orig": np.repeat(tfrac, np.shape(datay_cd)[1]),
+                               "facesvec": np.tile(facesvec, np.shape(datay_cd)[0]),
+                               "Curreny density electrolyte": datay_cd.flatten()
+                               })
+        # Build dataframes for bulkp and csld
+        dff_bulkp_c = pd.DataFrame()
+        dff_csld = pd.DataFrame()
+        # cstr can have varying length, determine maximum length
+        if config[trode, "type"] in constants.one_var_types:
+            partStr = "partTrode{trode}vol{vInd}part{pInd}" + sStr
+            cstr_base = pfx + partStr + "c"
+            maxlength = max([np.shape(utils.get_dict_key(data, cstr_base.format(
+                            trode=t, pInd=p, vInd=v)))[1]
+                for t in trodes for p in range(Npart[t]) for v in range(Nvol[t])])
+        else:
+            partStr = "partTrode{trode}vol{vInd}part{pInd}" + sStr
+            maxlength = max([np.shape(utils.get_dict_key(data, (pfx + partStr + "c1").format(
+                            trode=t, pInd=p, vInd=v)))[1]
+                for t in trodes for p in range(Npart[t]) for v in range(Nvol[t])])
+        for trode in trodes:
+            bulkp = pfx + 'phi_bulk_{trode}'.format(trode=trode)
+            dataybulkp = utils.get_dict_key(data, bulkp).flatten()
+            if trode == "a":
+                dataxbulkp = cellsvec[:Nvol["a"]]
+            elif trode == "c":
+                dataxbulkp = cellsvec[-Nvol["c"]:]
+            datat = np.repeat(timestd, len(dataxbulkp))
+            datatfrac = np.repeat(tfrac, len(dataxbulkp))
+            dataxbulkp = np.repeat([dataxbulkp], len(timestd), axis=0).flatten()
+            df_b = pd.DataFrame({
+                "Model": model,
+                "Time (s)": datat,
+                "Time fraction (%)": np.round(datatfrac),
+                "fraction orig": datatfrac,
+                "Trode": trode,
+                "Potential (nondim)": dataybulkp,
+                "Position in electrode": dataxbulkp})
+            dff_bulkp_c = pd.concat([dff_bulkp_c, df_b])
+            partStr = "partTrode{trode}vol{{vInd}}part{{pInd}}".format(trode=trode) + sStr
+            for pInd in range(Npart[trode]):
+                for vInd in range(Nvol[trode]):
+                    lens_str = "lens_{vInd}_{pInd}".format(vInd=vInd, pInd=pInd)
+                    if config[trode, "type"] in constants.one_var_types:
+                        cstr_base = pfx + partStr + "c"
+                        cstr = cstr_base.format(trode=trode, pInd=pInd, vInd=vInd)
+                        datay = np.empty([len(timestd), maxlength]) + np.nan
+                        yy = utils.get_dict_key(data, cstr)
+                        datay[0:len(timestd), 0:np.shape(yy)[1]] = yy
+                        datax = np.empty(maxlength) + np.nan
+                        datax[0:np.shape(yy)[1]] = np.linspace(0, psd_len[trode][vInd,pInd] * Lfac,
+                                                               np.shape(yy)[1])
+                        if trode == trodes[0] and pInd == 0 and vInd == 0:
+                            df_csld = pd.DataFrame({"Model": model,
+                                                    "sStr": sStr,
+                                                    "pfx": pfx,
+                                                    "Config trode type": config[trode, "type"],
+                                                    "Npartc": Npart["c"],
+                                                    "Nvolc": Nvol["c"],
+                                                    "Nparta": nparta,
+                                                    "Nvola": nvola,
+                                                    "time (s)": np.repeat(timestd,
+                                                                          np.shape(datay)[1]),
+                                                    "Time fraction": np.repeat(np.round(tfrac),
+                                                                               np.shape(datay)[1]),
+                                                    "fraction orig": np.repeat(tfrac,
+                                                                               np.shape(datay)[1]),
+                                                    lens_str: np.repeat([datax], len(datay),
+                                                                        axis=0).flatten(),
+                                                    cstr: datay.flatten()
+                                                    })
+                        else:
+                            if lens_str not in df_csld.columns.to_numpy():
+                                lens_str_data = np.repeat([datax], len(datay), axis=0).flatten()
+                                df_csld = pd.concat((df_csld, pd.DataFrame({lens_str:
+                                                                           lens_str_data})),
+                                                    axis=1)
+                            df_csld = pd.concat((df_csld, pd.DataFrame({cstr: datay.flatten()})),
+                                                axis=1)
+                    elif config[trode, "type"] in constants.two_var_types:
+                        c1str_base = pfx + partStr + "c1"
+                        c2str_base = pfx + partStr + "c2"
+                        c3str_base = pfx + partStr + "cav"
+                        c1str = c1str_base.format(trode=trode, pInd=pInd, vInd=vInd)
+                        c2str = c2str_base.format(trode=trode, pInd=pInd, vInd=vInd)
+                        c3str = c3str_base.format(trode=trode, pInd=pInd, vInd=vInd)
+                        datay1 = datay2 = datay3 = np.empty([len(timestd), maxlength]) + np.nan
+                        yy1 = utils.get_dict_key(data, c1str)
+                        datay1[0:len(timestd), 0:np.shape(yy1)[1]] = yy1
+                        yy2 = utils.get_dict_key(data, c2str)
+                        datay2[0:len(timestd), 0:np.shape(yy2)[1]] = yy2
+                        datay3 = 0.5*(datay1 + datay2)
+                        datax = np.empty(maxlength) + np.nan
+                        numy = np.shape(yy1)[1] if isinstance(yy1, np.ndarray) else 1
+                        datax[0:np.shape(yy1)[1]] = np.linspace(0,
+                                                                psd_len[trode][vInd,pInd] * Lfac,
+                                                                numy)
+                        if trode == trodes[0] and pInd == 0 and vInd == 0:
+                            df_csld = pd.DataFrame({"Model": model,
+                                                    "sStr": sStr,
+                                                    "pfx": pfx,
+                                                    "Config trode type": config[trode, "type"],
+                                                    "Npartc": Npart["c"],
+                                                    "Nvolc": Nvol["c"],
+                                                    "Nparta": nparta,
+                                                    "Nvola": nvola,
+                                                    "time (s)": np.repeat(timestd,
+                                                                          np.shape(datay1)[1]),
+                                                    "Time fraction":
+                                                        np.repeat(np.round(tfrac),
+                                                                  np.shape(datay1)[1]),
+                                                    "fraction orig":
+                                                        np.repeat(tfrac,
+                                                                  np.shape(datay1)[1]),
+                                                    lens_str: np.repeat([datax], len(datay1),
+                                                                        axis=0).flatten(),
+                                                    c1str: datay1.flatten(),
+                                                    c2str: datay2.flatten(),
+                                                    c3str: datay3.flatten()
+                                                    })
+                        else:
+                            if lens_str not in df_csld.columns.to_numpy():
+                                df_csld[lens_str] = np.repeat([datax], len(datay1),
+                                                              axis=0).flatten(),
+                            df_csld[c1str] = datay1.flatten()
+                            df_csld[c2str] = datay2.flatten()
+                            df_csld[c3str] = datay3.flatten()
+            dff_csld = pd.concat([dff_csld, df_csld], ignore_index=True)
+        # make subselection dataframe with one fraction per rounded fraction
+        for i in np.unique(dff_c["Time fraction"]):
+            df_sub = dff_c[dff_c["Time fraction"] == i]
+            md = 1.0
+            mdx = 0.0
+            for j in np.unique(df_sub["fraction orig"]):
+                dx = abs(i-j)
+                if dx < md:
+                    md = dx
+                    mdx = j
+            select = dff_c[dff_c["fraction orig"] == mdx]
+            dff_c_sub = pd.concat([dff_c_sub, select], ignore_index=True)
+            select = dff_cd[dff_cd["fraction orig"] == mdx]
+            dff_cd_sub = pd.concat([dff_cd_sub, select], ignore_index=True)
+            select = dff_csld[dff_csld["fraction orig"] == mdx]
+            dff_csld_sub = pd.concat([dff_csld_sub, select], ignore_index=True)
+            select = dff_bulkp_c[dff_bulkp_c["fraction orig"] == mdx]
+            dff_bulkp = pd.concat([dff_bulkp, select], ignore_index=True)
+    return dff, dff_c_sub, dff_cd_sub, dff_bulkp, dff_csld_sub, df_cbar
diff --git a/mpet/ports.py b/mpet/ports.py
index 6689fd1b..ce9e0d2d 100644
--- a/mpet/ports.py
+++ b/mpet/ports.py
@@ -21,3 +21,23 @@ def __init__(self, Name, PortType, Model, Description=""):
         self.phi_m = dae.daeVariable(
             "phi_m", elec_pot_t, self,
             "Electric potential in the e- conducting phase")
+
+
+class portFromParticle(dae.daePort):
+    def __init__(self, Name, PortType, Model, Description=""):
+        dae.daePort.__init__(self, Name, PortType, Model, Description)
+        self.dcbardt = dae.daeVariable(
+            "dcbardt", dae.no_t, self,
+            "Rate of particle filling")
+
+
+class portFromInterface(dae.daePort):
+    def __init__(self, Name, PortType, Model, Description=""):
+        dae.daePort.__init__(self, Name, PortType, Model, Description)
+        self.i0 = dae.daeVariable(
+            "i0", dae.no_t, self,
+            "Current density in first interface region volume")
+
+        self.Nm0 = dae.daeVariable(
+            "Nm0", dae.no_t, self,
+            "Flux in first interface region volume")
diff --git a/mpet/props_am.py b/mpet/props_am.py
index 41dd6398..9f377bf3 100644
--- a/mpet/props_am.py
+++ b/mpet/props_am.py
@@ -154,6 +154,19 @@ def non_homog_rect_fixed_csurf(self, y, ybar, B, kappa, ywet):
         muR_nh = -kappa*curv + B*(y - ybar)
         return muR_nh
 
+    def non_homog_rect_variational(self, y, ybar, B, kappa):
+        """ Helper function """
+        # the taylor expansion at the edges is used
+        N_2 = len(y)
+        ytmp = np.empty(N_2+2, dtype=object)
+        dxs = 1./N_2
+        ytmp[1:-1] = y
+        ytmp[0] = y[0] + np.diff(y)[0]*dxs + 0.5*np.diff(y,2)[0]*dxs**2
+        ytmp[-1] = y[-1] + np.diff(y)[-1]*dxs + 0.5*np.diff(y,2)[-1]*dxs**2
+        curv = np.diff(ytmp, 2)/(dxs**2)
+        muR_nh = -kappa*curv + B*(y - ybar)
+        return muR_nh
+
     def non_homog_round_wetting(self, y, ybar, B, kappa, beta_s, shape, r_vec):
         """ Helper function """
         dr = r_vec[1] - r_vec[0]
@@ -177,16 +190,26 @@ def general_non_homog(self, y, ybar):
             raise Exception("Unknown input type")
         if ("homog" not in ptype) and (N > 1):
             shape = self.get_trode_param("shape")
-            kappa = self.get_trode_param("kappa")
-            B = self.get_trode_param("B")
             if shape == "C3":
                 if mod1var:
+                    kappa = self.get_trode_param("kappa")
+                    B = self.get_trode_param("B")
+                    if self.get_trode_param("type") in ["ACR"]:
+                        cwet = self.get_trode_param("cwet")
+                        muR_nh = self.non_homog_rect_fixed_csurf(
+                            y, ybar, B, kappa, cwet)
+                    elif self.get_trode_param("type") in ["ACR_Diff"]:
+                        muR_nh = self.non_homog_rect_variational(
+                            y, ybar, B, kappa)
+                elif mod2var:
+                    kappa = self.get_trode_param("kappa")
+                    B = self.get_trode_param("B")
                     cwet = self.get_trode_param("cwet")
                     muR_nh = self.non_homog_rect_fixed_csurf(
                         y, ybar, B, kappa, cwet)
-                elif mod2var:
-                    raise NotImplementedError("no 2param C3 model known")
             elif shape in ["cylinder", "sphere"]:
+                kappa = self.get_trode_param("kappa")
+                B = self.get_trode_param("B")
                 beta_s = self.get_trode_param("beta_s")
                 r_vec = geo.get_unit_solid_discr(shape, N)[0]
                 if mod1var:
diff --git a/mpet/sim.py b/mpet/sim.py
index c9130bc4..15908633 100644
--- a/mpet/sim.py
+++ b/mpet/sim.py
@@ -46,7 +46,7 @@ def __init__(self, config, tScale=None):
     def SetUpParametersAndDomains(self):
         # Domains
         config = self.config
-        if config["have_separator"]:
+        if config["Nvol"]["s"]:
             self.m.DmnCell["s"].CreateArray(config["Nvol"]["s"])
         for tr in config["trodes"]:
             self.m.DmnCell[tr].CreateArray(config["Nvol"][tr])
@@ -55,6 +55,9 @@ def SetUpParametersAndDomains(self):
                 for j in range(config["Npart"][tr]):
                     self.m.particles[tr][i, j].Dmn.CreateArray(
                         int(config["psd_num"][tr][i,j]))
+                    if config[f"simInterface_{tr}"]:
+                        self.m.interfaces[tr][i, j].Dmn.CreateArray(
+                            int(config["Nvol_i"]))
 
     def SetUpVariables(self):
         config = self.config
@@ -87,6 +90,9 @@ def SetUpVariables(self):
                             part.cbar.SetInitialGuess(cs0)
                             for k in range(Nij):
                                 part.c.SetInitialCondition(k, cs0)
+                                if self.config["c","type"] in ["ACR_Diff"]:
+                                    part.c_left_GP.SetInitialGuess(cs0)
+                                    part.c_right_GP.SetInitialGuess(cs0)
                         elif solidType in constants.two_var_types:
                             part.c1bar.SetInitialGuess(cs0)
                             part.c2bar.SetInitialGuess(cs0)
@@ -118,7 +124,7 @@ def SetUpVariables(self):
                 self.m.phi_lyteGP_L.SetInitialGuess(0)
 
             # Separator electrolyte initialization
-            if config["have_separator"]:
+            if config["Nvol"]["s"]:
                 for i in range(Nvol["s"]):
                     self.m.c_lyte["s"].SetInitialCondition(i, config['c0'])
                     self.m.phi_lyte["s"].SetInitialGuess(i, 0)
@@ -132,6 +138,12 @@ def SetUpVariables(self):
                     # Set electrolyte concentration in each particle
                     for j in range(Npart[tr]):
                         self.m.particles[tr][i,j].c_lyte.SetInitialGuess(config["c0"])
+                        # Set concentration and potential in interface region
+                        if config[f"simInterface_{tr}"]:
+                            for k in range(config["Nvol_i"]):
+                                self.m.interfaces[tr][i,j].c.SetInitialCondition(k,
+                                                                                 config["c0_int"])
+                                self.m.interfaces[tr][i,j].phi.SetInitialGuess(k, 0)
 
         else:
             dPrev = self.dataPrev
@@ -174,7 +186,7 @@ def SetUpVariables(self):
                                     k, data[partStr + "c1"][-1,k])
                                 part.c2.SetInitialCondition(
                                     k, data[partStr + "c2"][-1,k])
-            if config["have_separator"]:
+            if config["Nvol"]["s"]:
                 for i in range(Nvol["s"]):
                     self.m.c_lyte["s"].SetInitialCondition(
                         i, data["c_lyte_s"][-1,i])
diff --git a/mpet/utils.py b/mpet/utils.py
index d74a51b4..0ed23ce7 100644
--- a/mpet/utils.py
+++ b/mpet/utils.py
@@ -70,6 +70,8 @@ def get_asc_vec(var, Nvol, dt=False):
             if isinstance(var[sectn], dae.pyCore.daeVariable):
                 varout[sectn] = get_var_vec(var[sectn], Nvol[sectn], dt)
             # Otherwise, it's a parameter that varies with electrode section
+            elif isinstance(var[sectn], np.ndarray):
+                varout[sectn] = var[sectn]
             else:
                 varout[sectn] = get_const_vec(var[sectn], Nvol[sectn])
         # Otherwise, fill with zeros
@@ -88,7 +90,7 @@ def get_dxvec(L, Nvol):
         dxa = Nvol["a"] * [L["a"]/Nvol["a"]]
     else:
         dxa = []
-    if "s" in Nvol:
+    if Nvol["s"]:
         dxs = Nvol["s"] * [L["s"]/Nvol["s"]]
     else:
         dxs = []
diff --git a/setup.py b/setup.py
index 7e9b1038..9f21f434 100644
--- a/setup.py
+++ b/setup.py
@@ -25,11 +25,13 @@
         'mpet.electrode.reactions',
         'mpet.electrolyte','mpet.config'
     ],
-    install_requires=['numpy','scipy','matplotlib','pyQt5', 'h5py', 'configparser', 'schema'],
+    install_requires=['numpy','scipy','matplotlib','pyQt5', 'h5py', 'configparser', 'schema',
+                      'dask-jobqueue', 'dash', 'dash_bootstrap_components', 'bokeh'],
     extras_require={'test':['pytest','coverage', 'coveralls', 'flake8'],
                     'doc':['sphinx','sphinx_rtd_theme']},
     python_requires='>=3.6',
-    scripts=['bin/mpetrun.py','bin/mpetplot.py'],
+    scripts=['bin/mpetrun.py','bin/mpetplot.py','bin/run_jobs.py', 'bin/create_ensemble.py',
+             'bin/mpet_create_runjobs_dashboard.py', 'bin/mpet_plot_app.py'],
     classifiers=[
         "Programming Language :: Python :: 3",
     ],


Test patch:
diff --git a/tests/compare_plots.py b/tests/compare_plots.py
index 54acf262..48b1eb66 100644
--- a/tests/compare_plots.py
+++ b/tests/compare_plots.py
@@ -62,7 +62,8 @@ def xyCmp(testDir, dirDict, ptype, xlbl, ylbl, ttl, fname):
     newSimOutDir = osp.join(testDir, "sim_output")
     oldSimOutDir = osp.join(dirDict["refs"], testName, "sim_output")
     plotsDir = dirDict["plots"]
-    kwargs = {"print_flag": False, "save_flag": False, "data_only": True}
+    kwargs = {"print_flag": False, "save_flag": False, "data_only": True,
+              "color_changes": None, "smooth_type": None}
     new_xy = pd.show_data(newSimOutDir, plot_type=ptype, **kwargs)
     old_xy = pd.show_data(oldSimOutDir, plot_type=ptype, **kwargs)
     scl = 1.3
@@ -87,7 +88,8 @@ def xyPartsCmp(testDir, dirDict, ptype, xlbl, ylbl, ttl, fname):
     newSimOutDir = osp.join(testDir, "sim_output")
     oldSimOutDir = osp.join(dirDict["refs"], testName, "sim_output")
     plotsDir = dirDict["plots"]
-    kwargs = {"print_flag": False, "save_flag": False, "data_only": True}
+    kwargs = {"print_flag": False, "save_flag": False, "data_only": True,
+              "color_changes": 'discrete', "smooth_type": None}
     new_dataDict = pd.show_data(newSimOutDir, plot_type=ptype, **kwargs)
     old_dataDict = pd.show_data(oldSimOutDir, plot_type=ptype, **kwargs)
     new_t = pd.show_data(newSimOutDir, plot_type="vt", **kwargs)[0]
diff --git a/tests/ref_outputs/test025/params_c.cfg b/tests/ref_outputs/test025/params_c.cfg
new file mode 100644
index 00000000..a65405ea
--- /dev/null
+++ b/tests/ref_outputs/test025/params_c.cfg
@@ -0,0 +1,31 @@
+# Default parameters for simulating LFP in 1D using the ACR model.
+# See params_electrodes.cfg for parameter explanations.
+
+[Particles]
+type = ACR
+discretization = 1e-9
+shape = C3
+thickness = 20e-9
+
+[Material]
+muRfunc = LiFePO4
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+Omega_a = 1.8560e-20
+kappa = 5.0148e-10
+B = 0.1916e9
+rho_s = 1.3793e28
+D = 1e-18
+Dfunc = lattice
+dgammadc = 1e-29
+cwet = 0.98
+
+[Reactions]
+rxnType = BV
+k0 = 0.6e-1
+E_A = 13000
+alpha = 0.5
+# Fraggedakis et al. 2020, lambda = 8.3kBT
+lambda = 3.4113e-20
+Rfilm = 0e-0
diff --git a/tests/ref_outputs/test025/params_system.cfg b/tests/ref_outputs/test025/params_system.cfg
new file mode 100644
index 00000000..e1b62aa7
--- /dev/null
+++ b/tests/ref_outputs/test025/params_system.cfg
@@ -0,0 +1,200 @@
+# 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 = 2
+# Voltage cutoffs, V
+Vmax = 3.6
+Vmin = 2.9
+# 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 = 1
+# Series resistance, [Ohm m^2]
+Rser = 0.
+# Cathode, anode, and separator numer disc. in x direction (volumes in electrodes)
+# - Nvol_c must be >= 1
+# - If Nvol_c = 1 & Nvol_a = 0 & Nvol_s = 0, simulate a single volume with no
+#   separator and infinitely fast counter electrode
+# - If Nvol_a = 0, simulate a Li foil electrode
+Nvol_c = 10
+Nvol_s = 20
+Nvol_a = 0
+# Number of particles per volume for cathode and anode
+Npart_c = 10
+Npart_a = 0
+
+[Electrodes]
+# The name of the parameter file describing the cathode particles
+cathode = params_c.cfg
+# The name of the parameter file describing the anode particles
+# Used only if Nvol_a > 0
+anode = params_a.cfg
+# Rate constant of the Li foil electrode, A/m^2
+# Used only if Nvol_a = 0
+k0_foil = 1e0
+# Film resistance on the Li foil, Ohm m^2
+Rfilm_foil = 0e-0
+
+[Particles]
+# electrode particle size distribution info, m
+#   C3 -- size along [100] direction
+#   sphere or cylinder -- radius
+# If using stddev = 0, set Npart = 1 for that electrode to avoid
+# wasting computational effort.
+mean_c = 150e-9
+stddev_c = 50e-9
+mean_a = 100e-9
+stddev_a = 1e-9
+# contact penalty
+fraction_of_contact = 0.5
+stand_dev_contact = 0.2
+# 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 = 0e-6
+L_s = 100e-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 = 0.8
+# Bruggeman exponent (tortuosity = porosity^bruggExp)
+BruggExp_c = -1.5
+BruggExp_a = -0.5
+BruggExp_s = -0.5
+
+[Electrolyte]
+# Initial electrolyte conc., mol/m^3  = cmax*delta
+# c0 and c0fa below are not needed
+c0 = 40000
+#Maximal concentration of host sites (in solid elyte) [mol m^-3]
+cmax = 80000
+
+# 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
+# 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: 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 = 5.5e-12
+Dm = 2.94e-15
+
+[Interface]
+# Simulate interface region ? 
+# Options: true, false (Default: false)
+simInterface_c = true
+simInterface_a = false
+Nvol_i = 4
+L_i = 20e-9
+BruggExp_i = -0.5
+poros_i = 1.0
+interfaceModelType = solid
+interfaceSMset = solid_elyte
+c0_int = 40000
+cmax_i = 80000
+Dp_i = 5.5e-17
+# Dm_i = 2.94e-19
diff --git a/tests/ref_outputs/test025/sim_output/commit.diff b/tests/ref_outputs/test025/sim_output/commit.diff
new file mode 100644
index 00000000..e04a8c49
--- /dev/null
+++ b/tests/ref_outputs/test025/sim_output/commit.diff
@@ -0,0 +1,60 @@
+diff --git a/configs/params_LFP.cfg b/configs/params_LFP.cfg
+index 1dc1a98..f763c88 100644
+--- a/configs/params_LFP.cfg
++++ b/configs/params_LFP.cfg
+@@ -16,14 +16,14 @@ Omega_a = 1.8560e-20
+ kappa = 5.0148e-10
+ B = 0.1916e9
+ rho_s = 1.3793e28
+-D = 5.3e-19
++D = 1e-18
+ Dfunc = lattice
+-dgammadc = 0e-30
++dgammadc = 1e-29
+ cwet = 0.98
+ 
+ [Reactions]
+ rxnType = BV
+-k0 = 1.6e-1
++k0 = 0.6e-1
+ E_A = 13000
+ alpha = 0.5
+ # Fraggedakis et al. 2020, lambda = 8.3kBT
+diff --git a/configs/params_system.cfg b/configs/params_system.cfg
+index f763a6e..e307f18 100644
+--- a/configs/params_system.cfg
++++ b/configs/params_system.cfg
+@@ -62,7 +62,7 @@ seed = 0
+ # printing internal variable concentrations) files. hdf5 files
+ # are better for cycling, as they store less information and there is less
+ # opening/rewriting of files. Default is mat
+-dataReporter = hdf5
++dataReporter = mat
+ # Series resistance, [Ohm m^2]
+ Rser = 0.
+ # Cathode, anode, and separator numer disc. in x direction (volumes in electrodes)
+@@ -209,13 +209,13 @@ Dm = 2.94e-10
+ # Simulate interface region?
+ # Options: true, false (Default: false)
+ simInterface_c = false
+-simInterface_a = false
++simInterface_a = true
+ c0_int = 40000
+ cmax_i = 80000
+ # for numerical dicretization
+ Nvol_i = 4
+ # interface length
+-L_i = 20e-9
++L_i = 20e-8
+ BruggExp_i = -0.5
+ poros_i = 1.0
+ # Options: dilute, SM
+@@ -234,5 +234,5 @@ interfaceModelType = solid
+ #   solid_elyte: uses the lattice gas model for solid electrolytes.
+ interfaceSMset = solid_elyte
+ # diff coefficinets of the interface, can also be liquid
+-Dp_i = 5.5e-17
++Dp_i = 5.5e-18
+ Dm_i = 2.94e-19
+\ No newline at end of file
+
diff --git a/tests/ref_outputs/test025/sim_output/daetools_config_options.txt b/tests/ref_outputs/test025/sim_output/daetools_config_options.txt
new file mode 100644
index 00000000..bd8d80ba
--- /dev/null
+++ b/tests/ref_outputs/test025/sim_output/daetools_config_options.txt
@@ -0,0 +1,137 @@
+{
+    "daetools": {
+        "core": {
+            "checkForInfiniteNumbers": "false",
+            "eventTolerance": "1E-7",
+            "logIndent": "    ",
+            "pythonIndent": "    ",
+            "checkUnitsConsistency": "true",
+            "resetLAMatrixAfterDiscontinuity": "true",
+            "printInfo": "false",
+            "nodes": {
+                "useNodeMemoryPools": "false",
+                "deleteNodesThreshold": "1000000"
+            },
+            "equations": {
+                "info": [
+                    "If simplifyExpressions is true equation expressions will be simplified.",
+                    "evaluationMode specifies the mode of evaluation: evaluationTree_OpenMP or computeStack_OpenMP.",
+                    "computeStack_External is set by specifying the adComputeStackEvaluator_t object to the simulation.",
+                    "If numThreads is 0 the default number of threads will be used (typically the number of cores in the system)."
+                ],
+                "simplifyExpressions": "false",
+                "evaluationMode": "computeStack_OpenMP",
+                "evaluationTree_OpenMP": {
+                    "numThreads": "0"
+                },
+                "computeStack_OpenMP": {
+                    "numThreads": "0"
+                }
+            }
+        },
+        "activity": {
+            "printHeader": "false",
+            "printStats": "false",
+            "timeHorizon": "100.0",
+            "reportingInterval": "1.0",
+            "reportTimeDerivatives": "false",
+            "reportSensitivities": "false",
+            "stopAtModelDiscontinuity": "true",
+            "reportDataAroundDiscontinuities": "true",
+            "objFunctionAbsoluteTolerance": "1E-8",
+            "constraintsAbsoluteTolerance": "1E-8",
+            "measuredVariableAbsoluteTolerance": "1E-8"
+        },
+        "datareporting": {
+            "tcpipDataReceiverAddress": "127.0.0.1",
+            "tcpipDataReceiverPort": "50000",
+            "tcpipNumberOfRetries": "10",
+            "tcpipRetryAfterMilliSecs": "1000"
+        },
+        "logging": {
+            "tcpipLogAddress": "127.0.0.1",
+            "tcpipLogPort": "51000"
+        },
+        "minlpsolver": {
+            "printInfo": "false"
+        },
+        "IDAS": {
+            "relativeTolerance": "1E-5",
+            "integrationMode": "Normal",
+            "reportDataInOneStepMode": "false",
+            "nextTimeAfterReinitialization": "1E-7",
+            "printInfo": "false",
+            "numberOfSTNRebuildsDuringInitialization": "1000",
+            "SensitivitySolutionMethod": "Staggered",
+            "SensErrCon": "false",
+            "sensRelativeTolerance": "1E-5",
+            "sensAbsoluteTolerance": "1E-5",
+            "MaxOrd": "5",
+            "MaxNumSteps": "1000",
+            "InitStep": "0.0",
+            "MaxStep": "0.0",
+            "MaxErrTestFails": "10",
+            "MaxNonlinIters": "4",
+            "MaxConvFails": "10",
+            "NonlinConvCoef": "0.33",
+            "SuppressAlg": "false",
+            "NoInactiveRootWarn": "false",
+            "NonlinConvCoefIC": "0.0033",
+            "MaxNumStepsIC": "5",
+            "MaxNumJacsIC": "4",
+            "MaxNumItersIC": "10",
+            "LineSearchOffIC": "false",
+            "gmres": {
+                "kspace": "30",
+                "EpsLin": "0.05",
+                "JacTimesVecFn": "DifferenceQuotient",
+                "DQIncrementFactor": "1.0",
+                "MaxRestarts": "5",
+                "GSType": "MODIFIED_GS"
+            }
+        },
+        "superlu": {
+            "factorizationMethod": "SamePattern_SameRowPerm",
+            "useUserSuppliedWorkSpace": "false",
+            "workspaceSizeMultiplier": "3.0",
+            "workspaceMemoryIncrement": "1.5"
+        },
+        "superlu_mt": {
+            "numThreads": "0"
+        },
+        "intel_pardiso": {
+            "numThreads": "0"
+        },
+        "BONMIN": {
+            "IPOPT": {
+                "print_level": "0",
+                "tol": "1E-5",
+                "linear_solver": "mumps",
+                "hessianApproximation": "limited-memory",
+                "mu_strategy": "adaptive"
+            }
+        },
+        "NLOPT": {
+            "printInfo": "false",
+            "xtol_rel": "1E-6",
+            "xtol_abs": "1E-6",
+            "ftol_rel": "1E-6",
+            "ftol_abs": "1E-6",
+            "constr_tol": "1E-6"
+        },
+        "deal_II": {
+            "printInfo": "false",
+            "assembly": {
+                "info": [
+                    "parallelAssembly can be: Sequential or OpenMP.",
+                    "If numThreads is 0 the default number of threads will be used (typically the number of cores in the system).",
+                    "queueSize specifies the size of the internal queue; when this size is reached the local data are copied to the global matrices."
+                ],
+                "parallelAssembly": "OpenMP",
+                "numThreads": "0",
+                "queueSize": "32"
+            }
+        }
+    }
+}
+
diff --git a/tests/ref_outputs/test025/sim_output/input_dict_cathode.p b/tests/ref_outputs/test025/sim_output/input_dict_cathode.p
new file mode 100644
index 00000000..756ceb66
Binary files /dev/null and b/tests/ref_outputs/test025/sim_output/input_dict_cathode.p differ
diff --git a/tests/ref_outputs/test025/sim_output/input_dict_derived_values.p b/tests/ref_outputs/test025/sim_output/input_dict_derived_values.p
new file mode 100644
index 00000000..668e40b9
Binary files /dev/null and b/tests/ref_outputs/test025/sim_output/input_dict_derived_values.p differ
diff --git a/tests/ref_outputs/test025/sim_output/input_dict_system.p b/tests/ref_outputs/test025/sim_output/input_dict_system.p
new file mode 100644
index 00000000..3c42a6dc
Binary files /dev/null and b/tests/ref_outputs/test025/sim_output/input_dict_system.p differ
diff --git a/tests/ref_outputs/test025/sim_output/output_data b/tests/ref_outputs/test025/sim_output/output_data
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/ref_outputs/test025/sim_output/output_data.mat b/tests/ref_outputs/test025/sim_output/output_data.mat
new file mode 100644
index 00000000..08a9781b
Binary files /dev/null and b/tests/ref_outputs/test025/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test025/sim_output/run_info.txt b/tests/ref_outputs/test025/sim_output/run_info.txt
new file mode 100644
index 00000000..d7260a6d
--- /dev/null
+++ b/tests/ref_outputs/test025/sim_output/run_info.txt
@@ -0,0 +1,17 @@
+mpet version:
+0.1.8
+
+branch name:
+feature/solid
+
+commit hash:
+05c7393
+
+to run, from the root repo directory, copy relevant files there,
+edit input_params_system.cfg to point to correct material
+params files, and:
+$ git checkout [commit hash]
+$ patch -p1 < commit.diff:
+$ python[3] mpetrun.py input_params_system.cfg
+
+Total run time: 14.246185541152954 s
diff --git a/tests/ref_outputs/test026/params_c.cfg b/tests/ref_outputs/test026/params_c.cfg
new file mode 100644
index 00000000..a29c7791
--- /dev/null
+++ b/tests/ref_outputs/test026/params_c.cfg
@@ -0,0 +1,25 @@
+#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
+noise = false
+noise_prefac = 1e-6
+numnoise = 200
+rho_s = 3.1036e28
+D = 1e-14
+Dfunc = constant
+E_D = 5000
+
+[Reactions]
+rxnType = BV_mod01
+k0 = 3.671
+E_A = 5000
+alpha = 0.5
+Rfilm = 0e-0
diff --git a/tests/ref_outputs/test026/params_system.cfg b/tests/ref_outputs/test026/params_system.cfg
new file mode 100644
index 00000000..2d5f6b2a
--- /dev/null
+++ b/tests/ref_outputs/test026/params_system.cfg
@@ -0,0 +1,200 @@
+# 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 = -3
+# Voltage cutoffs, V
+Vmax = 4.6
+Vmin = 2.8
+# 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 = 1
+# Series resistance, [Ohm m^2]
+Rser = 0.
+# Cathode, anode, and separator numer disc. in x direction (volumes in electrodes)
+# - Nvol_c must be >= 1
+# - If Nvol_c = 1 & Nvol_a = 0 & Nvol_s = 0, simulate a single volume with no
+#   separator and infinitely fast counter electrode
+# - If Nvol_a = 0, simulate a Li foil electrode
+Nvol_c = 10
+Nvol_s = 20
+Nvol_a = 0
+# Number of particles per volume for cathode and anode
+Npart_c = 10
+Npart_a = 0
+
+[Electrodes]
+# The name of the parameter file describing the cathode particles
+cathode = params_c.cfg
+# The name of the parameter file describing the anode particles
+# Used only if Nvol_a > 0
+anode = params_a.cfg
+# Rate constant of the Li foil electrode, A/m^2
+# Used only if Nvol_a = 0
+k0_foil = 1e0
+# Film resistance on the Li foil, Ohm m^2
+Rfilm_foil = 0e-0
+
+[Particles]
+# electrode particle size distribution info, m
+#   C3 -- size along [100] direction
+#   sphere or cylinder -- radius
+# If using stddev = 0, set Npart = 1 for that electrode to avoid
+# wasting computational effort.
+mean_c = 500e-9
+stddev_c = 100e-9
+mean_a = 100e-9
+stddev_a = 1e-9
+# contact penalty
+fraction_of_contact = 1
+stand_dev_contact = 0
+# Initial electrode filling fractions
+# (for disch, anode starts full, cathode starts empty)
+cs0_c = 0.99
+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 = 0e-6
+L_s = 100e-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 = 0.8
+# Bruggeman exponent (tortuosity = porosity^bruggExp)
+BruggExp_c = -1.5
+BruggExp_a = -0.5
+BruggExp_s = -0.5
+
+[Electrolyte]
+# Initial electrolyte conc., mol/m^3  = cmax*delta
+# c0 and c0fa below are not needed
+c0 = 40000
+#Maximal concentration of host sites (in solid elyte) [mol m^-3]
+cmax = 80000
+
+# 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
+# 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: 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 = 5.5e-12
+Dm = 2.94e-15
+
+[Interface]
+# Simulate interface region ? 
+# Options: true, false (Default: false)
+simInterface_c = true
+simInterface_a = false
+Nvol_i = 4
+L_i = 20e-9
+BruggExp_i = -0.5
+poros_i = 1.0
+interfaceModelType = solid
+interfaceSMset = solid_elyte
+c0_int = 40000
+cmax_i = 80000
+Dp_i = 5.5e-19
+# Dm_i = 2.94e-19
diff --git a/tests/ref_outputs/test026/sim_output/commit.diff b/tests/ref_outputs/test026/sim_output/commit.diff
new file mode 100644
index 00000000..e04a8c49
--- /dev/null
+++ b/tests/ref_outputs/test026/sim_output/commit.diff
@@ -0,0 +1,60 @@
+diff --git a/configs/params_LFP.cfg b/configs/params_LFP.cfg
+index 1dc1a98..f763c88 100644
+--- a/configs/params_LFP.cfg
++++ b/configs/params_LFP.cfg
+@@ -16,14 +16,14 @@ Omega_a = 1.8560e-20
+ kappa = 5.0148e-10
+ B = 0.1916e9
+ rho_s = 1.3793e28
+-D = 5.3e-19
++D = 1e-18
+ Dfunc = lattice
+-dgammadc = 0e-30
++dgammadc = 1e-29
+ cwet = 0.98
+ 
+ [Reactions]
+ rxnType = BV
+-k0 = 1.6e-1
++k0 = 0.6e-1
+ E_A = 13000
+ alpha = 0.5
+ # Fraggedakis et al. 2020, lambda = 8.3kBT
+diff --git a/configs/params_system.cfg b/configs/params_system.cfg
+index f763a6e..e307f18 100644
+--- a/configs/params_system.cfg
++++ b/configs/params_system.cfg
+@@ -62,7 +62,7 @@ seed = 0
+ # printing internal variable concentrations) files. hdf5 files
+ # are better for cycling, as they store less information and there is less
+ # opening/rewriting of files. Default is mat
+-dataReporter = hdf5
++dataReporter = mat
+ # Series resistance, [Ohm m^2]
+ Rser = 0.
+ # Cathode, anode, and separator numer disc. in x direction (volumes in electrodes)
+@@ -209,13 +209,13 @@ Dm = 2.94e-10
+ # Simulate interface region?
+ # Options: true, false (Default: false)
+ simInterface_c = false
+-simInterface_a = false
++simInterface_a = true
+ c0_int = 40000
+ cmax_i = 80000
+ # for numerical dicretization
+ Nvol_i = 4
+ # interface length
+-L_i = 20e-9
++L_i = 20e-8
+ BruggExp_i = -0.5
+ poros_i = 1.0
+ # Options: dilute, SM
+@@ -234,5 +234,5 @@ interfaceModelType = solid
+ #   solid_elyte: uses the lattice gas model for solid electrolytes.
+ interfaceSMset = solid_elyte
+ # diff coefficinets of the interface, can also be liquid
+-Dp_i = 5.5e-17
++Dp_i = 5.5e-18
+ Dm_i = 2.94e-19
+\ No newline at end of file
+
diff --git a/tests/ref_outputs/test026/sim_output/daetools_config_options.txt b/tests/ref_outputs/test026/sim_output/daetools_config_options.txt
new file mode 100644
index 00000000..bd8d80ba
--- /dev/null
+++ b/tests/ref_outputs/test026/sim_output/daetools_config_options.txt
@@ -0,0 +1,137 @@
+{
+    "daetools": {
+        "core": {
+            "checkForInfiniteNumbers": "false",
+            "eventTolerance": "1E-7",
+            "logIndent": "    ",
+            "pythonIndent": "    ",
+            "checkUnitsConsistency": "true",
+            "resetLAMatrixAfterDiscontinuity": "true",
+            "printInfo": "false",
+            "nodes": {
+                "useNodeMemoryPools": "false",
+                "deleteNodesThreshold": "1000000"
+            },
+            "equations": {
+                "info": [
+                    "If simplifyExpressions is true equation expressions will be simplified.",
+                    "evaluationMode specifies the mode of evaluation: evaluationTree_OpenMP or computeStack_OpenMP.",
+                    "computeStack_External is set by specifying the adComputeStackEvaluator_t object to the simulation.",
+                    "If numThreads is 0 the default number of threads will be used (typically the number of cores in the system)."
+                ],
+                "simplifyExpressions": "false",
+                "evaluationMode": "computeStack_OpenMP",
+                "evaluationTree_OpenMP": {
+                    "numThreads": "0"
+                },
+                "computeStack_OpenMP": {
+                    "numThreads": "0"
+                }
+            }
+        },
+        "activity": {
+            "printHeader": "false",
+            "printStats": "false",
+            "timeHorizon": "100.0",
+            "reportingInterval": "1.0",
+            "reportTimeDerivatives": "false",
+            "reportSensitivities": "false",
+            "stopAtModelDiscontinuity": "true",
+            "reportDataAroundDiscontinuities": "true",
+            "objFunctionAbsoluteTolerance": "1E-8",
+            "constraintsAbsoluteTolerance": "1E-8",
+            "measuredVariableAbsoluteTolerance": "1E-8"
+        },
+        "datareporting": {
+            "tcpipDataReceiverAddress": "127.0.0.1",
+            "tcpipDataReceiverPort": "50000",
+            "tcpipNumberOfRetries": "10",
+            "tcpipRetryAfterMilliSecs": "1000"
+        },
+        "logging": {
+            "tcpipLogAddress": "127.0.0.1",
+            "tcpipLogPort": "51000"
+        },
+        "minlpsolver": {
+            "printInfo": "false"
+        },
+        "IDAS": {
+            "relativeTolerance": "1E-5",
+            "integrationMode": "Normal",
+            "reportDataInOneStepMode": "false",
+            "nextTimeAfterReinitialization": "1E-7",
+            "printInfo": "false",
+            "numberOfSTNRebuildsDuringInitialization": "1000",
+            "SensitivitySolutionMethod": "Staggered",
+            "SensErrCon": "false",
+            "sensRelativeTolerance": "1E-5",
+            "sensAbsoluteTolerance": "1E-5",
+            "MaxOrd": "5",
+            "MaxNumSteps": "1000",
+            "InitStep": "0.0",
+            "MaxStep": "0.0",
+            "MaxErrTestFails": "10",
+            "MaxNonlinIters": "4",
+            "MaxConvFails": "10",
+            "NonlinConvCoef": "0.33",
+            "SuppressAlg": "false",
+            "NoInactiveRootWarn": "false",
+            "NonlinConvCoefIC": "0.0033",
+            "MaxNumStepsIC": "5",
+            "MaxNumJacsIC": "4",
+            "MaxNumItersIC": "10",
+            "LineSearchOffIC": "false",
+            "gmres": {
+                "kspace": "30",
+                "EpsLin": "0.05",
+                "JacTimesVecFn": "DifferenceQuotient",
+                "DQIncrementFactor": "1.0",
+                "MaxRestarts": "5",
+                "GSType": "MODIFIED_GS"
+            }
+        },
+        "superlu": {
+            "factorizationMethod": "SamePattern_SameRowPerm",
+            "useUserSuppliedWorkSpace": "false",
+            "workspaceSizeMultiplier": "3.0",
+            "workspaceMemoryIncrement": "1.5"
+        },
+        "superlu_mt": {
+            "numThreads": "0"
+        },
+        "intel_pardiso": {
+            "numThreads": "0"
+        },
+        "BONMIN": {
+            "IPOPT": {
+                "print_level": "0",
+                "tol": "1E-5",
+                "linear_solver": "mumps",
+                "hessianApproximation": "limited-memory",
+                "mu_strategy": "adaptive"
+            }
+        },
+        "NLOPT": {
+            "printInfo": "false",
+            "xtol_rel": "1E-6",
+            "xtol_abs": "1E-6",
+            "ftol_rel": "1E-6",
+            "ftol_abs": "1E-6",
+            "constr_tol": "1E-6"
+        },
+        "deal_II": {
+            "printInfo": "false",
+            "assembly": {
+                "info": [
+                    "parallelAssembly can be: Sequential or OpenMP.",
+                    "If numThreads is 0 the default number of threads will be used (typically the number of cores in the system).",
+                    "queueSize specifies the size of the internal queue; when this size is reached the local data are copied to the global matrices."
+                ],
+                "parallelAssembly": "OpenMP",
+                "numThreads": "0",
+                "queueSize": "32"
+            }
+        }
+    }
+}
+
diff --git a/tests/ref_outputs/test026/sim_output/input_dict_cathode.p b/tests/ref_outputs/test026/sim_output/input_dict_cathode.p
new file mode 100644
index 00000000..b854b331
Binary files /dev/null and b/tests/ref_outputs/test026/sim_output/input_dict_cathode.p differ
diff --git a/tests/ref_outputs/test026/sim_output/input_dict_derived_values.p b/tests/ref_outputs/test026/sim_output/input_dict_derived_values.p
new file mode 100644
index 00000000..e3ffd631
Binary files /dev/null and b/tests/ref_outputs/test026/sim_output/input_dict_derived_values.p differ
diff --git a/tests/ref_outputs/test026/sim_output/input_dict_system.p b/tests/ref_outputs/test026/sim_output/input_dict_system.p
new file mode 100644
index 00000000..d8a3e635
Binary files /dev/null and b/tests/ref_outputs/test026/sim_output/input_dict_system.p differ
diff --git a/tests/ref_outputs/test026/sim_output/output_data b/tests/ref_outputs/test026/sim_output/output_data
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/ref_outputs/test026/sim_output/output_data.mat b/tests/ref_outputs/test026/sim_output/output_data.mat
new file mode 100644
index 00000000..595ebf25
Binary files /dev/null and b/tests/ref_outputs/test026/sim_output/output_data.mat differ
diff --git a/tests/ref_outputs/test026/sim_output/run_info.txt b/tests/ref_outputs/test026/sim_output/run_info.txt
new file mode 100644
index 00000000..b27fc81f
--- /dev/null
+++ b/tests/ref_outputs/test026/sim_output/run_info.txt
@@ -0,0 +1,17 @@
+mpet version:
+0.1.8
+
+branch name:
+feature/solid
+
+commit hash:
+05c7393
+
+to run, from the root repo directory, copy relevant files there,
+edit input_params_system.cfg to point to correct material
+params files, and:
+$ git checkout [commit hash]
+$ patch -p1 < commit.diff:
+$ python[3] mpetrun.py input_params_system.cfg
+
+Total run time: 1.1276493072509766 s