diff --git a/MaterialLib/SolidModels/MFront/PowerLawLinearCreep.mfront b/MaterialLib/SolidModels/MFront/PowerLawLinearCreep.mfront
index 4f945dc4d7296c85f4692e1934ae627ff195389e..e53eab8889f1d0c34c7d7764f1cceff88e65f6d2 100644
--- a/MaterialLib/SolidModels/MFront/PowerLawLinearCreep.mfront
+++ b/MaterialLib/SolidModels/MFront/PowerLawLinearCreep.mfront
@@ -8,7 +8,7 @@ model for salt creep. Bérest et al. (2019). Rock Mech Rock Eng.
 @Algorithm NewtonRaphson;
 @MaximumNumberOfIterations 100;
 
-@Epsilon 1.e-14;
+@Epsilon 1.e-12;
 @Theta 1.0;
 
 @ModellingHypotheses{".+"};
@@ -51,23 +51,20 @@ Ru.setEntryName("UniversalGasConstant");
 @Integrator
 {
     const auto s = deviator(sig);
-    const auto norm_s = sigmaeq(sig) / std::sqrt(3. / 2.);
+    const auto norm_s = std::max(sigmaeq(sig),1e-14*mu) / std::sqrt(3. / 2.);
     constexpr auto Pdev = Stensor4::K();
 
     const auto bPL = std::pow(3. / 2., (m + 1) / 2) * A1 *
                      std::exp(-Q1 / Ru / (T + dT)) / std::pow(sig0, m);
-    const auto bL = 3. / 2. * A2 / std::pow(Dgrain, 3) *
+    const auto bL = 3. / 2. * A2 / (std::pow(Dgrain, 3) * (T + dT)) *
                     std::exp(-Q2 / Ru / (T + dT)) / sig0;
 
-    if (norm_s > 1e-14 * mu)
-    {
-        const auto norm_s_pow = std::pow(norm_s, m - 1);
-        depsPL = dt * bPL * norm_s_pow * s;
-        depsL = dt * bL * s;
-        feel += depsPL + depsL;
-        dfeel_ddeel +=
-            2. * mu * dt *
-            (bPL * norm_s_pow * (Pdev + (((m - 1) / norm_s / norm_s) * s ^ s)) +
-             bL * Pdev);
-    }
+    const auto norm_s_pow = std::pow(norm_s, m - 1);
+    depsPL = dt * bPL * norm_s_pow * s;
+    depsL = dt * bL * s;
+    feel += depsPL + depsL;
+    dfeel_ddeel +=
+        2. * mu * dt *
+        (bPL * norm_s_pow * (Pdev + (((m - 1) / norm_s / norm_s) * s ^ s)) +
+         bL * Pdev);
 }
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index db3edc67e295d044f6ad04c7fe245066f6b7ce68..6a7f1ecf480ed51d4041292603a34819110fef56 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -95,6 +95,7 @@ if (OGS_USE_MFRONT)
     if(NOT OGS_USE_MPI)
         OgsTest(PROJECTFILE Mechanics/GuentherSalzer/model_triaxtest.prj)
     endif()
+    OgsTest(PROJECTFILE Mechanics/PLLC/uniax_compression.prj)
 
     # Linear elastic, no internal state variables, no external state variables.
     AddTest(
@@ -266,4 +267,7 @@ if(NOT OGS_USE_PETSC)
     if(NOT WIN32)
         NotebookTest(NOTEBOOKFILE Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb RUNTIME 40)
     endif()
+    if (OGS_USE_MFRONT)
+        NotebookTest(NOTEBOOKFILE Mechanics/PLLC/PLLC.ipynb RUNTIME 7)
+    endif()
 endif()
diff --git a/Tests/Data/Mechanics/PLLC/PLLC.ipynb b/Tests/Data/Mechanics/PLLC/PLLC.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..b129c90a90480b8794243a42dcfb9843ecc80c9b
--- /dev/null
+++ b/Tests/Data/Mechanics/PLLC/PLLC.ipynb
@@ -0,0 +1,231 @@
+{
+ "cells": [
+  {
+   "cell_type": "raw",
+   "id": "bb0907b4-4e26-4c4e-ab1f-22b5330cb1d2",
+   "metadata": {},
+   "source": [
+    "title = \"Power Law Linear Creep\"\n",
+    "date = \"2023-01-02\"\n",
+    "author = \"Florian Zill\"\n",
+    "notebook = \"Mechanics/PLLC/PLLC.ipynb\"\n",
+    "web_subsection = \"small-deformations\"\n",
+    "<!--eofm-->"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "fa8cd9b5",
+   "metadata": {},
+   "source": [
+    "### Power Law Linear Creep\n",
+    "\n",
+    "This benchmark shows the increased creep rate of salt rock at lower deviatoric stress. A two component power law (which we call Power Law Linear Creep, or short PLLC) provides an easy way to capture the power law behaviour (dislocation creep) and the linear behaviour (pressure solution creep). For more details have a look at (Zill et al., 2022).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7962f42f-fd53-4fc1-b966-a8ba924aca6c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import contextlib\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib as mpl\n",
+    "import numpy as np\n",
+    "import os\n",
+    "import vtuIO\n",
+    "from ogs6py import ogs\n",
+    "\n",
+    "prj_name = \"uniax_compression\"\n",
+    "data_dir = os.environ.get('OGS_DATA_DIR', str(os.getcwd()).split(\"/Data/\")[0] + \"/Data/\")\n",
+    "input_file = f\"{data_dir}/Mechanics/PLLC/{prj_name}.prj\"\n",
+    "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', f'{data_dir}/Mechanics/PLLC/_out')\n",
+    "\n",
+    "if not os.path.exists(out_dir):\n",
+    "    os.makedirs(out_dir)\n",
+    "os.chdir(out_dir)\n",
+    "\n",
+    "prj_file = f\"{out_dir}/{prj_name}_out.prj\"\n",
+    "ogs_model = ogs.OGS(INPUT_FILE=input_file, PROJECT_FILE=prj_file)\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "7250c5b0",
+   "metadata": {},
+   "source": [
+    "### Experimental data\n",
+    "\n",
+    "A nice overview for the strain rates of salt for different temperatures and differential stresses can be found in (Li et al., 2021)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5a20a14e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Unfortunately the source for the WIPP data has gone missing - will be added if it's found again\n",
+    "ExData = {\n",
+    "    \"WIPP CS 25\": (25, \"^\", [[9.87970002, 2.013560846E-05], [11.84642707, 3.178356756E-05], [7.87388785, 1.66059726E-06]]),\n",
+    "    \"WIPP CS 60\": (60, \"^\", [[3.98589289, 5.7824853E-06], [5.94266985, 2.075776623E-05], [7.87388785, 1.953209818E-05], [9.96978837, 5.841438703E-05], [11.84642707, 0.00011762092257], [13.94911482, 0.00026749321794], [17.9857158, 0.00111804208073], [1.9814251, 8.7645834E-07], [3.91418422, 4.01350889E-06], [5.88897108, 3.34371363E-06], [7.87388785, 1.129440706E-05], [9.87970002, 2.99068674E-05], [11.84642707, 7.681792203E-05], [13.82306874, 0.00011067584933], [15.83934389, 0.00052247037957]]),\n",
+    "    \"DeVries 1988 25\": (25, \"s\", [[4.99, 2.10816E-06], [4.99, 2.4192E-06], [5, 1.8144E-06], [9.99, 2.2032E-05], [14.96, 9.2448E-05], [14.98, 0.000216]]),\n",
+    "    \"DeVries 1988 100\": (100, \"s\", [[4.95, 9.6768E-05], [6.77, 0.000292896], [7.46, 0.000324], [8.55, 0.000664416], [8.92, 0.00091584], [8.98, 0.0009936], [9.91, 0.00124416], [10.1, 0.00139968], [10.22, 0.00093312], [10.27, 0.00132192], [12.1, 0.00216], [12.3, 0.00409536], [12.35, 0.00320544], [12.37, 0.00292032], [12.39, 0.00253152], [12.4, 0.0026784], [12.46, 0.0025056], [12.49, 0.00347328], [13.57, 0.00273024], [13.78, 0.00242784], [14.7, 0.00482112], [16.87, 0.0095904], [17.2, 0.0123552], [19.96, 0.030672]]),\n",
+    "    \"DeVries 1988 200\": (200, \"s\", [[3.47, 0.00117504], [4.71, 0.0032832], [6.67, 0.0104544], [6.78, 0.0132192], [9.86, 0.214272]]),\n",
+    "    \"Berest 2015 14.3\": (14.3, \"P\", [[0.09909639, 8.944207E-08], [0.19575886, 1.4118213E-07], [0.29452325, 1.4118213E-07], [0.49411031, 9.799173E-08]]),\n",
+    "    \"Berest 2017 7.8\": (7.8, \"P\", [[0.19575886,2.2285256E-07], [0.19575886,9.505469E-08], [0.19754389,2.5947583E-07], [0.19754389,2.647936E-08], [0.39379426,4.9162047E-07], [0.39738509,6.801413E-08], [0.59247161,4.0957628E-07], [0.59247161,5.7241269E-07], [0.59787408,1.0735864E-07], [1.0591736,1.11804208E-06]])}"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "6d38a65c",
+   "metadata": {},
+   "source": [
+    "### Parameters\n",
+    "\n",
+    "This set of parameters gives a good fit with the experimental data. The grain size is a bit larger than the usual grain size of roughly 1 cm."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8066a6d3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "A1 = 0.18 # d^-1\n",
+    "Q1 = 54e3 # kJ / mol\n",
+    "A2 = 6.5e-5 # m^3 K d^−1\n",
+    "Q2 = 24.5e3 # kJ / mol\n",
+    "dGrain = 5e-2 # m\n",
+    "sref = 1. # MPa\n",
+    "BGRa = lambda sig, T: A1 * np.exp(-Q1/(8.3145*(273.15+T))) * np.power(sig/sref,5.)\n",
+    "PLLC = lambda sig, T: A1 * np.exp(-Q1/(8.3145*(273.15+T))) * np.power(sig/sref,5.) + \\\n",
+    "                      A2 * np.exp(-Q2/(8.3145*(273.15+T))) * sig/sref / np.power(dGrain, 3) / (273.15+T)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "f1e5ec96",
+   "metadata": {},
+   "source": [
+    "### Simulation and plot\n",
+    "\n",
+    "The experimental data is compared against the model results (analytically and numerically)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7e2e294c-e803-4f02-b5ab-9bdfef94b00f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lo_stresses = np.array([0.2e6, 0.6e6])\n",
+    "hi_stresses = np.array([2e6, 10e6])\n",
+    "Exps = {7.8: ('blue', lo_stresses), 14.3: ('orange', lo_stresses),\n",
+    "        25: ('lime', hi_stresses), 60: ('red', hi_stresses),\n",
+    "        100: ('gray', hi_stresses), 200: ('mediumpurple', hi_stresses)}\n",
+    "\n",
+    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
+    "ax.set_xlabel('$\\\\sigma_\\\\mathrm{ax}$ / MPa')\n",
+    "ax.set_ylabel('$\\\\dot{\\\\epsilon}_{zz}$ / d$^{-1}$')\n",
+    "ax.set_xlim(0.15, 30)\n",
+    "ax.set_ylim(1e-15, 1e1)\n",
+    "ax.grid(visible=True, which='both')\n",
+    "points = {'pt0': (1., 1., 1.)}\n",
+    "\n",
+    "sigs = np.logspace(-1, 2, 100)\n",
+    "for temp, (col, stresses) in Exps.items():   \n",
+    "    # plot analytical curves\n",
+    "    if temp >= 25:\n",
+    "        ax.plot(sigs, BGRa(sigs, temp), color=col, ls='--')\n",
+    "        ax.plot(sigs, PLLC(sigs, temp), color=col, ls='-')\n",
+    "\n",
+    "    # simulation in ogs and plot results\n",
+    "    eps_dot = []\n",
+    "    ogs_model.replace_parameter_value(\"T_ref\", str(temp + 273.15))\n",
+    "    for stress in stresses:\n",
+    "        ogs_model.replace_parameter_value(\"sigma_ax\", str(-stress))\n",
+    "        ogs_model.write_input()\n",
+    "        # hide output\n",
+    "        with contextlib.redirect_stdout(None):\n",
+    "            ogs_model.run_model(logfile=f\"{out_dir}/out.txt\", \n",
+    "                                args=\"-m \" + f\"{data_dir}/Mechanics/PLLC/\")\n",
+    "            pvdfile = vtuIO.PVDIO(f\"{prj_name}.pvd\", dim=3)\n",
+    "        eps_zz = pvdfile.read_time_series(\"epsilon\", points)[\"pt0\"][:, 2]\n",
+    "        eps_zz_dot = np.abs(np.diff(eps_zz)) / np.diff(pvdfile.timesteps)\n",
+    "        # omit the first timestep\n",
+    "        eps_dot += [np.mean(eps_zz_dot[1:])]\n",
+    "    ax.loglog(1e-6*stresses, eps_dot, 'o', c=col, markeredgecolor=\"k\")\n",
+    "\n",
+    "# plot experimental data points\n",
+    "for Ex, (temp, m, Data) in ExData.items():\n",
+    "    stresses, eps_dot = np.array(Data).T\n",
+    "    ax.loglog(stresses, eps_dot, m, c=Exps[temp][0])\n",
+    "\n",
+    "# create legend\n",
+    "patches = [mpl.patches.Patch(color=col, label=str(temp) + '°C')\n",
+    "           for temp, (col, _) in Exps.items() if temp >= 25][::-1]\n",
+    "addLeg = lambda **args : patches.append(mpl.lines.Line2D([], [], **args))\n",
+    "addLeg(c='k', label='PLLC')\n",
+    "addLeg(c='k', ls='--', label='BGRa')\n",
+    "addLeg(c='w', ls='None', marker='o', mec=\"k\", label='OGS')\n",
+    "addLeg(c='k', ls='None', marker='s', label='DeVries (1988)')\n",
+    "addLeg(c='k', ls='None', marker='^', label='WIPP CS')\n",
+    "addLeg(c='b', ls='None', marker='P', label='Bérest (2017) 7.8°C')\n",
+    "addLeg(c='orange', ls='None', marker='P', label='Bérest (2015) 14.3°C')\n",
+    "ax.legend(handles=patches, loc='best')\n",
+    "\n",
+    "fig.tight_layout()\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "id": "05f56358",
+   "metadata": {},
+   "source": [
+    "### References\n",
+    "\n",
+    "Zill, Florian, Wenqing Wang, and Thomas Nagel. Influence of THM Process Coupling and Constitutive Models on the Simulated Evolution of Deep Salt Formations during Glaciation. The Mechanical Behavior of Salt X. CRC Press, 2022. https://doi.org/10.1201/9781003295808-33.\n",
+    "\n",
+    "Li, Shiyuan, and Janos Urai. Numerical Studies of the Deformation of Salt Bodies with Embedded Carbonate Stringers. Online, print. Publikationsserver der RWTH Aachen University, 2012. http://publications.rwth-aachen.de/record/211523/files/4415.pdf "
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/Tests/Data/Mechanics/PLLC/cube_1x1x1.gml b/Tests/Data/Mechanics/PLLC/cube_1x1x1.gml
new file mode 100644
index 0000000000000000000000000000000000000000..6ff4684bcbc304c16240037249ddeff6c7b692ea
--- /dev/null
+++ b/Tests/Data/Mechanics/PLLC/cube_1x1x1.gml
@@ -0,0 +1,44 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<?xml-stylesheet type="text/xsl" href="OpenGeoSysGLI.xsl"?>
+
+<OpenGeoSysGLI xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ogs="http://www.opengeosys.org">
+    <name>cube_1x1x1_geometry</name>
+    <points>
+        <point id="0" x="0" y="0" z="0"/>
+        <point id="1" x="0" y="0" z="1"/>
+        <point id="2" x="0" y="1" z="1"/>
+        <point id="3" x="0" y="1" z="0"/>
+
+        <point id="4" x="1" y="0" z="0"/>
+        <point id="5" x="1" y="0" z="1"/>
+        <point id="6" x="1" y="1" z="1"/>
+        <point id="7" x="1" y="1" z="0"/>
+    </points>
+
+    <surfaces>
+        <surface id="0" name="left"><!-- x=0 -->
+            <element p1="0" p2="1" p3="2"/>
+            <element p1="0" p2="3" p3="2"/>
+        </surface>
+        <surface id="1" name="right"><!-- x=1 -->
+            <element p1="4" p2="6" p3="5"/>
+            <element p1="4" p2="6" p3="7"/>
+        </surface>
+        <surface id="2" name="top"><!-- z=1 -->
+            <element p1="1" p2="2" p3="5"/>
+            <element p1="5" p2="2" p3="6"/>
+        </surface>
+        <surface id="3" name="bottom"><!-- z=0 -->
+            <element p1="0" p2="3" p3="4"/>
+            <element p1="4" p2="3" p3="7"/>
+        </surface>
+        <surface id="4" name="front"><!-- y=0 -->
+            <element p1="0" p2="1" p3="4"/>
+            <element p1="4" p2="1" p3="5"/>
+        </surface>
+        <surface id="5" name="back"><!-- y=1 -->
+            <element p1="2" p2="3" p3="6"/>
+            <element p1="6" p2="3" p3="7"/>
+        </surface>
+    </surfaces>
+</OpenGeoSysGLI>
diff --git a/Tests/Data/Mechanics/PLLC/cube_1x1x1_hex_1e0.vtu b/Tests/Data/Mechanics/PLLC/cube_1x1x1_hex_1e0.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..3d9f172f2db01e6e865cc58b038a48a40372ab28
--- /dev/null
+++ b/Tests/Data/Mechanics/PLLC/cube_1x1x1_hex_1e0.vtu
@@ -0,0 +1,33 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="8" NumberOfCells="1">
+      <PointData>
+      </PointData>
+      <CellData>
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="1.7320508076">
+          AQAAAACAAADAAAAAHAAAAA==eJxjYMAHPtjjlcaQh/GJ1YdLPy51mDQAp2EONQ==
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="binary" RangeMin="0" RangeMax="7">
+          AQAAAACAAABAAAAAHgAAAA==eJxjYIAARijNDKWZoDQLlGaF0uxQmg1KAwAC8AAd
+        </DataArray>
+        <DataArray type="Int64" Name="offsets" format="binary" RangeMin="8" RangeMax="8">
+          AQAAAACAAAAIAAAACwAAAA==eJzjYIAAAABIAAk=
+        </DataArray>
+        <DataArray type="UInt8" Name="types" format="binary" RangeMin="12" RangeMax="12">
+          AQAAAACAAAABAAAACQAAAA==eJzjAQAADQAN
+        </DataArray>
+        <DataArray type="Int64" Name="faces" format="binary" RangeMin="0" RangeMax="0">
+          AQAAAACAAAAIAAAACwAAAA==eJxjYIAAAAAIAAE=
+        </DataArray>
+        <DataArray type="Int64" Name="faceoffsets" format="binary" RangeMin="1" RangeMax="1">
+          AQAAAACAAAAIAAAACwAAAA==eJxjZIAAAAAQAAI=
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/Tests/Data/Mechanics/PLLC/uniax_compression.prj b/Tests/Data/Mechanics/PLLC/uniax_compression.prj
new file mode 100644
index 0000000000000000000000000000000000000000..56e422837acb037f69933cb431cf5e34d4d0c57a
--- /dev/null
+++ b/Tests/Data/Mechanics/PLLC/uniax_compression.prj
@@ -0,0 +1,240 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+	<mesh>cube_1x1x1_hex_1e0.vtu</mesh>
+	<geometry>cube_1x1x1.gml</geometry>
+    <processes>
+        <process>
+            <name>SD</name>
+            <type>SMALL_DEFORMATION</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation>
+                <type>MFront</type>
+                <behaviour>PowerLawLinearCreep</behaviour>
+                <material_properties>
+                    <material_property name="YoungModulus" parameter="E"/>
+                    <material_property name="PoissonRatio" parameter="nu"/>
+                    <material_property name="PowerLawFactor" parameter="A"/>
+                    <material_property name="PowerLawEnergy" parameter="Q"/>
+                    <material_property name="PowerLawExponent" parameter="n"/>
+                    <material_property name="ReferenceStress" parameter="sigma_f"/>
+                    <material_property name="LinearLawFactor" parameter="A2"/>
+                    <material_property name="LinearLawEnergy" parameter="Q2"/>
+                    <material_property name="SaltGrainSize" parameter="Dgrain"/>
+                </material_properties>
+            </constitutive_relation>
+            <solid_density>rho_sr</solid_density>
+            <specific_body_force>0 0 0</specific_body_force>
+			<reference_temperature>T_ref</reference_temperature>
+            <process_variables>
+                <process_variable>displacement</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="ElasticStrain" output_name="ElasticStrain"/>
+                <secondary_variable internal_name="sigma" output_name="sigma"/>
+                <secondary_variable internal_name="epsilon" output_name="epsilon"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="SD">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-12</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>2</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>uniax_compression</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>1</each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable>displacement</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+                <variable>ElasticStrain</variable>
+            </variables>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>T_ref</name>
+            <type>Constant</type>
+            <value>298.15</value>
+        </parameter>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>25e9</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>0.25</value>
+        </parameter>
+        <parameter>
+            <name>A</name>
+            <type>Constant</type>
+            <value>0.18</value>
+        </parameter>
+        <parameter>
+            <name>n</name>
+            <type>Constant</type>
+            <value>5.0</value>
+        </parameter>
+        <parameter>
+            <name>sigma_f</name>
+            <type>Constant</type>
+            <value>1e6</value>
+        </parameter>
+        <parameter>
+            <name>Q</name>
+            <type>Constant</type>
+            <value>54000</value>
+        </parameter>
+        <parameter>
+            <name>A2</name>
+            <type>Constant</type>
+            <value>6.5e-5</value>
+        </parameter>
+        <parameter>
+            <name>Q2</name>
+            <type>Constant</type>
+            <value>24500</value>
+        </parameter>
+        <parameter>
+            <name>Dgrain</name>
+            <type>Constant</type>
+            <value>5e-2</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>2570.0</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0 0</values>
+        </parameter>
+        <parameter>
+            <name>sigma_ax</name>
+            <type>Constant</type>
+            <value>-1e6</value>
+        </parameter>
+        <parameter>
+            <name>zero</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>3</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+				<boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>zero</parameter>
+                </boundary_condition>
+				<boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>front</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>zero</parameter>
+                </boundary_condition>
+				<boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>2</component>
+                    <parameter>zero</parameter>
+                </boundary_condition>
+				<boundary_condition>
+                    <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Neumann</type>
+                    <component>2</component>
+                    <parameter>sigma_ax</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>100</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <eigen>
+                <solver_type>SparseLU</solver_type>
+            <scaling>true</scaling>
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+    <test_definition>
+        <vtkdiff>
+            <file>uniax_compression_ts_2_t_2.000000.vtu</file>
+            <field>displacement</field>
+            <absolute_tolerance>1e-13</absolute_tolerance>
+            <relative_tolerance>1e-14</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>uniax_compression_ts_2_t_2.000000.vtu</file>
+            <field>NodalForces</field>
+            <absolute_tolerance>1e-9</absolute_tolerance>
+            <relative_tolerance>1e-11</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>uniax_compression_ts_2_t_2.000000.vtu</file>
+            <field>sigma</field>
+            <absolute_tolerance>1e-8</absolute_tolerance>
+            <relative_tolerance>1e-12</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>uniax_compression_ts_2_t_2.000000.vtu</file>
+            <field>epsilon</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>1e-15</relative_tolerance>
+        </vtkdiff>
+        <vtkdiff>
+            <file>uniax_compression_ts_2_t_2.000000.vtu</file>
+            <field>ElasticStrain</field>
+            <absolute_tolerance>1e-14</absolute_tolerance>
+            <relative_tolerance>1e-15</relative_tolerance>
+        </vtkdiff>
+    </test_definition>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Mechanics/PLLC/uniax_compression_ts_2_t_2.000000.vtu b/Tests/Data/Mechanics/PLLC/uniax_compression_ts_2_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..f17b12516ea03ae5775173ecbf292a548a8303e5
--- /dev/null
+++ b/Tests/Data/Mechanics/PLLC/uniax_compression_ts_2_t_2.000000.vtu
@@ -0,0 +1,37 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+  <UnstructuredGrid>
+    <FieldData>
+      <DataArray type="Int8" Name="IntegrationPointMetaData" NumberOfTuples="97" format="appended" RangeMin="34"                   RangeMax="125"                  offset="0"                   />
+      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="25" format="appended" RangeMin="45"                   RangeMax="121"                  offset="160"                 />
+      <DataArray type="Float64" Name="sigma_ip" NumberOfComponents="6" NumberOfTuples="8" format="appended" RangeMin="1000000"              RangeMax="1000000"              offset="248"                 />
+    </FieldData>
+    <Piece NumberOfPoints="8"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="Float64" Name="ElasticStrain" NumberOfComponents="6" format="appended" RangeMin="4.2426406871e-05"     RangeMax="4.2426406871e-05"     offset="516"                 />
+        <DataArray type="Float64" Name="MaterialForces" NumberOfComponents="3" format="appended" RangeMin="3.1622776602e+149"    RangeMax="-nan"                 offset="936"                 />
+        <DataArray type="Float64" Name="NodalForces" NumberOfComponents="3" format="appended" RangeMin="250000"               RangeMax="250000"               offset="1004"                />
+        <DataArray type="Float64" Name="displacement" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="4.2636355248e-05"     offset="1280"                />
+        <DataArray type="Float64" Name="epsilon" NumberOfComponents="6" format="appended" RangeMin="4.2636355248e-05"     RangeMax="4.2636355248e-05"     offset="1404"                />
+        <DataArray type="Float64" Name="sigma" NumberOfComponents="6" format="appended" RangeMin="1000000"              RangeMax="1000000"              offset="1792"                />
+      </PointData>
+      <CellData>
+        <DataArray type="Float64" Name="principal_stress_values" NumberOfComponents="3" format="appended" RangeMin="1000000"              RangeMax="1000000"              offset="2348"                />
+        <DataArray type="Float64" Name="principal_stress_vector_1" NumberOfComponents="3" format="appended" RangeMin="1"                    RangeMax="1"                    offset="2436"                />
+        <DataArray type="Float64" Name="principal_stress_vector_2" NumberOfComponents="3" format="appended" RangeMin="1"                    RangeMax="1"                    offset="2512"                />
+        <DataArray type="Float64" Name="principal_stress_vector_3" NumberOfComponents="3" format="appended" RangeMin="1"                    RangeMax="1"                    offset="2580"                />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.7320508076"         offset="2656"                />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="2740"                />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="2824"                />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="2884"                />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _AQAAAAAAAAAAgAAAAAAAAGEAAAAAAAAAVwAAAAAAAAA=eF5VjDEOgDAIAP/C7OTg0K8YQ1DRMAANrYNp+ne7Ot5dcg3EKt9BVdww+yCkCHoLpLX9osfJAWmewEgZEhS5lVAyDPPozoF+4eGa3djqGCx96x8LVCQuAQAAAAAAAAAAgAAAAAAAABkAAAAAAAAAIQAAAAAAAAA=eF4z0zPRM9Y1s7TUTTcwNDSzNE420kvJLCqpBABIQQaGAQAAAAAAAAAAgAAAAAAAAIABAAAAAAAApgAAAAAAAAA=eF5jYGBg2FGx0hZIOSS17Nj7HwjqW/QOHgwx/tvmusoWQt/aC6FP7WXAq74Lqn4WlN4FUv/g3N7DIPUM4vqXbf9iVX9r78T6Z7aqWxDqYPpg6u9s2bm2NmKpLUTdZqj6xSB1E9zqr8BpdPczQAGS+7Gqh5oH9ecqOA3S+/LCFbA/ihiP72UEMhqQ3A8zH8n9BTsbLu+F0Uwo6rNsYebC3A8AtESthw==AQAAAAAAAAAAgAAAAAAAAIABAAAAAAAAGAEAAAAAAAA=eF77mvG4Y+uPJ3ZfoPQPMM2y38ptYeGv+6nWBfm8R34mlew2Oy7qbiQ0d/c7qLo3KDTL/gkJmz7PmTHbushpzzqxokbrecFHzX8llFrD1MPo11D1na0TVv6tXbxb9IYf37N3Xdb1CW9eP4pasxtm7lsUfSz7P01f7cF/PtH6uviaow94F+9WUt7zZqX8butXONzz57DeBI6Vwbs/cpu1/KrcYx3JWjwt6NAG60dQdc+g9COo+mfHxByv9tfvflZa8uBvz27rzcG96ekXT++GhctXKP0Jqj5qeXnHXdlF1j/1RXLMeQ7ubjf8dtd7zibrb1B1n6H0N6j6/5MvGT83mGm9mUEg7c/ifbt/b1p686br+d0Ar2ro5w==AQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAAEQAAAAAAAAA=eF5jYACBH/UMQ5QGAN48Iyk=AQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAArAAAAAAAAAA=eF67dmKuyVK7SXtjrN4173xQbcsABA0tfI6/xFyOXahv3XvwwSHh22ZLbf/9//+/Hij+o/XhOuYTPXsPr9j0+PikqXDx2AT95JvTN9jeemOffOrFKri4tLSDody+ZluRrff9FzQv3gs1/+CCdwUT9xzSs/35KaHF8sf6vX8g6g/eO+Ch5r19xV7blOuhsacb9kLNOai9wnGudugE2y0LevftvdkFFwcAU81rOA==AQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAAPAAAAAAAAAA=eF5jYMAOwn/e9mnVeGqHQ5ohFE0exo8ioE954/2suQKs+2PQ1GlAxWF8mP1qUPEIKD8SSqtCxQFZUzCCAQAAAAAAAAAAgAAAAAAAAIABAAAAAAAAAQEAAAAAAAA=eF6L+Hnbp1XjqV0UlFbbeD9rrgDr/lu/7gR3vAm1ruCQqJ1QkG59wro+bd6Cst2hUHUhUFoJqv4m24yn745P373d9Vn0w1X91lV6LrbTZKZao6tXhKoviPxj/aQhzjr3j8Px1dem7z4ENR+mLgjN/NQXmz7W7Z1orbT1YeCse627K6Dmw9QFQmk5qPqlxbGTmy0KdzN57upfplKx+9+VDZObHCKsfaDqfKG0BFT9xtl8/HOOlO5+0r3o4IMtfdanoO6JhqqLgNIaUPUz1i7eb3NqqXXC/CPvEg2CrM9AzYeFYwyU1oaqr5rDHMP7uNnaRmmN8oK2adYHoeYDAMN4w/k=AQAAAAAAAAAAgAAAAAAAAIABAAAAAAAAgAEAAAAAAAA=eF47+ov1wXy/N3u9ilgWaYm82cvNwMDQ0KJ3cK6pk/5pxibbPJ1n81zfbd8b0+P1UOPEhL1M955EVGS+shVefcy3X/Su7bf////XA9XvPrPocfH2bbZdy9/vt5850fbLMqWZp3fMsb1z++Kj2zNe2MYFX+go2vTK9gdUfZCJwaeV3Af2PkriexcffnzvpHPP7fZrL7PluHDi22e1S3vZHp5m4ntze+8/qPo1ec1/bgfW2nK19+58HHzN9pjhytv74/bvfTd99pm21Fd7S/hnPXqr93jvH6h6+yMvayX0C/Y+kDI6mj39tO3eYw6yPwqu23p5xKxyu3PddtEa/df3FV7ZvoCqL3jPXaf5e+beN905K42fPdvrciDn7qO112w/KEqbnxN9YZu/SqbUOeCFLRM0fLbfOSw26/Je23cmTUcOLD5nq1l/dJuw4J29B0Re2KedOLJ3Ve9vx1MCD/dyQNU3LF6yfc6eLbZK8cIfnfa93PtztVw+i+XNvQAadNRWAQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAIQAAAAAAAAA=eF77//////oWvYMMDAwTunk22Zqe+DzrU+RZWwCrcgzkAQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAGAAAAAAAAAA=eF5jYACDBgbe1EsrtmfYQLgf7AExkwT3AQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAEAAAAAAAAAA=eF5jYACBD/YMaAAAFScBMA==AQAAAAAAAAAAgAAAAAAAABgAAAAAAAAAFgAAAAAAAAA=eF5jYEAGH+wZeFMvrdieYQMAGFMEdw==AQAAAAAAAAAAgAAAAAAAAMAAAAAAAAAAHAAAAAAAAAA=eF5jYMAHPtjjlcaQh/GJ1YdLPy51mDQAp2EONQ==AQAAAAAAAAAAgAAAAAAAAEAAAAAAAAAAHgAAAAAAAAA=eF5jYIAARijNDKWZoDQLlGaF0uxQmg1KAwAC8AAdAQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAACwAAAAAAAAA=eF7jYIAAAABIAAk=AQAAAAAAAAAAgAAAAAAAAAEAAAAAAAAACQAAAAAAAAA=eF7jAQAADQAN
+  </AppendedData>
+</VTKFile>