diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 850cfb5980ed89eb41911256631aa5d91d6af769..7e38bcd13f4f64c94b394a9eac3cf09feab2d48d 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -99,7 +99,7 @@ class LDDsimulation(object):
 
         ## Private variables
         # maximal number of L-iterations that the LDD solver uses.
-        self._max_iter_num = 300
+        self._max_iter_num = 5000
         # TODO rewrite this with regard to the mesh sizes
         self.calc_tol = self.tol
         # The number of subdomains are counted by self.init_meshes_and_markers()
@@ -572,7 +572,7 @@ class LDDsimulation(object):
                 #"th. initial mass": df.assemble(rho*df.dx), \
                 #"energy": self.total_energy(), \
                 "wetting": subsequent_errors["wetting"],
-                "nonwetting": subsequent_errors["nonwetting"]})
+                "nonwetting": subsequent_errors["nonwetting"]}, index=[iteration])
         # if self._rank == 0:
         # check if file exists
         if os.path.isfile(filename) == True:
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 6bb78942124064b50cdb358b65018f08d9dfbf41..deae3286a6d7348064bce6ce9852fd64a5c58975 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -777,8 +777,8 @@ class DomainPatch(df.SubDomain):
                 )
             else:
                 self._interface_has_common_dof_indices[interf_ind].update(
-                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w].copy()},
-                    {'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw].copy()}
+                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w].copy(),
+                    'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw].copy()}
                 )
 
 
diff --git a/TP-TP-patch-test-case/TP-TP-2-patch-test.py b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
index c61dff9c4b907843a6272aa7b9711271548d1675..1e9a7d5ba16aa07b41b614412d6a965f9b7230c2 100755
--- a/TP-TP-patch-test-case/TP-TP-2-patch-test.py
+++ b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
@@ -9,6 +9,9 @@ import LDDsimulation as ldd
 import functools as ft
 #import ufl as ufl
 
+# init sympy session
+sym.init_printing()
+
 ##### Domain and Interface ####
 # global simulation domain domain
 sub_domain0_vertices = [df.Point(0.0,0.0), #
@@ -80,19 +83,21 @@ outer_boundary_def_points = {
 # interface i (i.e. given by interface_def_points[i]).
 adjacent_subdomains = [[1,2]]
 isRichards = {
-    1: True, #
-    2: True
+    1: False, #
+    2: False
     }
 
 Tmax = 1
-timestep_size = 0.005
+timestep_size = 4*0.001
 
 time_interval = [0, Tmax]
 
 viscosity = {#
 # subdom_num : viscosity
-    1 : {'wetting' :1}, #
-    2 : {'wetting' :1}
+    1 : {'wetting' :1,
+         'nonwetting': 1/50}, #
+    2 : {'wetting' :1,
+         'nonwetting': 1/50}
 }
 
 porosity = {#
@@ -103,36 +108,49 @@ porosity = {#
 
 L = {#
 # subdom_num : subdomain L for L-scheme
-    1 : {'wetting' :0.25},#
-    2 : {'wetting' :0.25}
+    1 : {'wetting' :0.25,
+         'nonwetting': 0.25},#
+    2 : {'wetting' :0.25,
+         'nonwetting': 0.25}
 }
 
 lambda_param = {#
 # subdom_num : lambda parameter for the L-scheme
-    1 : {'wetting' :100},#
-    2 : {'wetting' :100}
+    1 : {'wetting' :100,
+         'nonwetting': 100},#
+    2 : {'wetting' :100,
+         'nonwetting': 100}
 }
 
 ## relative permeabilty functions on subdomain 1
-def rel_perm1(s):
-    # relative permeabilty on subdomain1
+def rel_perm1w(s):
+    # relative permeabilty wetting on subdomain1
     return s**2
 
-_rel_perm1 = ft.partial(rel_perm1)
+def rel_perm1nw(s):
+    # relative permeabilty nonwetting on subdomain1
+    return (1-s)**3
 
+_rel_perm1w = ft.partial(rel_perm1w)
+_rel_perm1nw = ft.partial(rel_perm1nw)
 subdomain1_rel_perm = {
-    'wetting': _rel_perm1,#
-    'nonwetting': None
+    'wetting': _rel_perm1w,#
+    'nonwetting': _rel_perm1nw
 }
 ## relative permeabilty functions on subdomain 2
-def rel_perm2(s):
-    # relative permeabilty on subdomain2
+def rel_perm2w(s):
+    # relative permeabilty wetting on subdomain2
     return s**3
-_rel_perm2 = ft.partial(rel_perm2)
+def rel_perm2nw(s):
+    # relative permeabilty nonwetting on subdomain2
+    return (1-s)**2
+
+_rel_perm2w = ft.partial(rel_perm2w)
+_rel_perm2nw = ft.partial(rel_perm2nw)
 
 subdomain2_rel_perm = {
-    'wetting': _rel_perm2,#
-    'nonwetting': None
+    'wetting': _rel_perm2w,#
+    'nonwetting': _rel_perm2nw
 }
 
 ## dictionary of relative permeabilties on all domains.
@@ -141,28 +159,129 @@ relative_permeability = {#
     2: subdomain2_rel_perm
 }
 
-def saturation(pressure, subdomain_index):
+# S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
+# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
+def saturation(capillary_pressure, n_index, alpha):
+    # inverse capillary pressure-saturation-relationship
+    return df.conditional(capillary_pressure < 0, 1/((1 + (alpha*capillary_pressure)**n_index)**((n_index - 1)/n_index)), 1)
+
+# S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
+# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
+def saturation_sym(capillary_pressure, n_index, alpha):
     # inverse capillary pressure-saturation-relationship
-    return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+    #df.conditional(capillary_pressure > 0,
+    return 1/((1 + (alpha*capillary_pressure)**n_index)**((n_index - 1)/n_index))
 
-sat_pressure_relationship = {#
-    1: ft.partial(saturation, subdomain_index = 1),#
-    2: ft.partial(saturation, subdomain_index = 2)
+S_pc_rel = {#
+    1: ft.partial(saturation_sym, n_index = 3, alpha=0.001),# n= 3 stands for non-uniform porous media
+    2: ft.partial(saturation_sym, n_index = 6, alpha=0.001) # n=6 stands for uniform porous media matrix (siehe Helmig)
 }
 
+S_pc_rel_sym = {#
+    1: ft.partial(saturation_sym, n_index = sym.Symbol('n'), alpha = sym.Symbol('a')),# n= 3 stands for non-uniform porous media
+    2: ft.partial(saturation_sym, n_index = sym.Symbol('n'), alpha = sym.Symbol('a')) # n=6 stands for uniform porous media matrix (siehe Helmig)
+}
+
+#### Manufacture source expressions with sympy
+###############################################################################
+## subdomain1
+x, y = sym.symbols('x[0], x[1]') # needed by UFL
+t = sym.symbols('t', positive=True)
+#f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)  # -Laplace(u)
+#f = sym.simplify(f)                         # simplify f
+p1_w = 1 - (1+t)*(1 + x**2 + y**2)
+p1_nw = t*(1-y)**2 - sym.sqrt(2+t**2)*(1-y)
+
+#dtS1_w = sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)
+#dtS1_nw = -sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)
+dtS1_w = porosity[1]*sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)
+dtS1_nw = -porosity[1]*sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)
+print("dtS1_w = ", dtS1_w, "\n")
+print("dtS1_nw = ", dtS1_nw, "\n")
+
+#dxdxflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)
+#dydyflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)
+dxdxflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)
+dydyflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)
+
+rhs1_w = dtS1_w + dxdxflux1_w + dydyflux1_w
+rhs1_w = sym.printing.ccode(rhs1_w)
+print("rhs_w = ", rhs1_w, "\n")
+#rhs_w = sym.expand(rhs_w)
+#print("rhs_w", rhs_w, "\n")
+#rhs_w = sym.collect(rhs_w, x)
+#print("rhs_w", rhs_w, "\n")
+
+#dxdxflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)
+#dydyflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)
+dxdxflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)
+dydyflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)
+
+rhs1_nw = dtS1_nw + dxdxflux1_nw + dydyflux1_nw
+rhs1_nw = sym.printing.ccode(rhs1_nw)
+print("rhs_nw = ", rhs1_nw, "\n")
+
+## subdomain2
+p2_w = 1 - (1+t)*(1 + x**2 + 0.25*y)
+p2_nw = t*(1-y)**2 - sym.sqrt(2+t**2)*(1-y) +(0.5-y)*(t/(0.1+x))
+
+#dtS2_w = sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)
+#dtS2_nw = -sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)
+dtS2_w = porosity[2]*sym.diff(S_pc_rel[2](p2_nw - p2_w), t, 1)
+dtS2_nw = -porosity[2]*sym.diff(S_pc_rel[2](p2_nw - p2_w), t, 1)
+print("dtS2_w = ", dtS2_w, "\n")
+print("dtS2_nw = ", dtS2_nw, "\n")
+
+#dxdxflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, x, 1), x, 1)
+#dydyflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, y, 1), y, 1)
+dxdxflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_w, x, 1), x, 1)
+dydyflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_w, y, 1), y, 1)
+
+rhs2_w = dtS2_w + dxdxflux2_w + dydyflux2_w
+rhs2_w = sym.printing.ccode(rhs2_w)
+print("rhs2_w = ", rhs2_w, "\n")
+#rhs_w = sym.expand(rhs_w)
+#print("rhs_w", rhs_w, "\n")
+#rhs_w = sym.collect(rhs_w, x)
+#print("rhs_w", rhs_w, "\n")
+
+#dxdxflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, x, 1), x, 1)
+#dydyflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, y, 1), y, 1)
+dxdxflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_nw, x, 1), x, 1)
+dydyflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_nw, y, 1), y, 1)
+
+rhs2_nw = dtS2_nw + dxdxflux2_nw + dydyflux2_nw
+rhs2_nw = sym.printing.ccode(rhs2_nw)
+print("rhs2_nw = ", rhs2_nw, "\n")
+
+
+###############################################################################
+
 source_expression = {
-    1: {'wetting': '4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )'},
-    2: {'wetting': '2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))'}
+    1: {'wetting': rhs1_w,
+        'nonwetting': rhs1_nw},
+    2: {'wetting': rhs2_w,
+        'nonwetting': rhs2_nw}
 }
 
+p1_w_00 = p1_w.subs(t, 0)
+p1_nw_00 = p1_nw.subs(t, 0)
+p2_w_00 = p2_w.subs(t, 0)
+p2_nw_00 = p2_nw.subs(t, 0)
+# p1_w_00 = sym.printing.ccode(p1_w_00)
+
 initial_condition = {
-    1: {'wetting': '-(x[0]*x[0] + x[1]*x[1])'},#
-    2: {'wetting': '-x[1]*x[1]'}
+    1: {'wetting': sym.printing.ccode(p1_w_00),
+        'nonwetting': sym.printing.ccode(p1_nw_00)},#
+    2: {'wetting': sym.printing.ccode(p2_w_00),
+        'nonwetting': sym.printing.ccode(p2_nw_00)}
 }
 
 exact_solution = {
-    1: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'},#
-    2: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}
+    1: {'wetting': sym.printing.ccode(p1_w),
+        'nonwetting': sym.printing.ccode(p1_nw)},#
+    2: {'wetting': sym.printing.ccode(p2_w),
+        'nonwetting': sym.printing.ccode(p2_nw)}
 }
 
 # similary to the outer boundary dictionary, if a patch has no outer boundary
@@ -175,8 +294,10 @@ exact_solution = {
 # the actual expression needed for the dirichlet condition for both phases if present.
 dirichletBC = {
 #subdomain index: {outer boudary part index: {phase: expression}}
-    1: { 0: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'}},
-    2: { 0: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}}
+    1: { 0: {'wetting': sym.printing.ccode(p1_w),
+             'nonwetting': sym.printing.ccode(p1_nw)}},
+    2: { 0: {'wetting': sym.printing.ccode(p2_w),
+             'nonwetting': sym.printing.ccode(p2_nw)}}
 }
 
 # def saturation(pressure, subdomain_index):
@@ -190,7 +311,7 @@ write_to_file = {
     'L_iterations': True
 }
 
-mesh_resolution = 30
+mesh_resolution = 10
 
 # initialise LDD simulation class
 simulation = ldd.LDDsimulation(tol = 1E-12)
@@ -206,7 +327,7 @@ simulation.set_parameters(output_dir = "./output/",#
     L = L,#
     lambda_param = lambda_param,#
     relative_permeability = relative_permeability,#
-    saturation = sat_pressure_relationship,#
+    saturation = S_pc_rel,#
     time_interval = time_interval,#
     timestep_size = timestep_size,#
     sources = source_expression,#
diff --git a/TP-TP-patch-test-case/exact_solution_symbolic_computation.ipynb b/TP-TP-patch-test-case/exact_solution_symbolic_computation.ipynb
index 2ff9bfb5c4d58fefd74d6b616b28127104669ba7..745a7b0e0f016422b2b577e3f8480f370712a3b1 100644
--- a/TP-TP-patch-test-case/exact_solution_symbolic_computation.ipynb
+++ b/TP-TP-patch-test-case/exact_solution_symbolic_computation.ipynb
@@ -9,7 +9,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -28,7 +28,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -112,25 +112,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "('dtS1_w', 0, '\\n')\n",
-      "('dtS1_nw', 0, '\\n')\n",
-      "('rhs_w', 4*t + 4, '\\n')\n",
-      "('rhs_w', 4*t + 4, '\\n')\n",
-      "('dxdxflux1_w = ', '2*t + 2', '\\n')\n",
-      "('dydyflux1_w = ', '2*t + 2', '\\n')\n",
-      "('dxdxflux1_nw = ', '0', '\\n')\n",
-      "('dydyflux1_nw = ', '0', '\\n')\n"
+      "dtS1_w =  -0.666666666666667*((0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + x[1]**2 + 1) - 0.001)**3 + 1)**(-1.66666666666667)*(0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + x[1]**2 + 1) - 0.001)**2*(-0.003*t*(1 - x[1])/sqrt(t**2 + 2) + 0.003*x[0]**2 + 0.003*x[1]**2 + 0.003*(1 - x[1])**2 + 0.003) \n",
+      "\n",
+      "dtS1_nw =  0.666666666666667*((0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + x[1]**2 + 1) - 0.001)**3 + 1)**(-1.66666666666667)*(0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + x[1]**2 + 1) - 0.001)**2*(-0.003*t*(1 - x[1])/sqrt(t**2 + 2) + 0.003*x[0]**2 + 0.003*x[1]**2 + 0.003*(1 - x[1])**2 + 0.003) \n",
+      "\n",
+      "rhs_w =  16.0*pow(x[0], 2)*(-t - 1)*(0.001*t + 0.001)*pow(pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 3) + 1, -2.333333333333333)*pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 2) + 2.6666666666666665*x[1]*(-t - 1)*pow(pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 3) + 1, -2.333333333333333)*(0.0030000000000000001*t*(2*x[1] - 2) + 6*x[1]*(0.001*t + 0.001) + 0.0030000000000000001*sqrt(pow(t, 2) + 2))*pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 2) - 4*(-t - 1)*pow(pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 3) + 1, -1.3333333333333333) - 0.66666666666666663*pow(pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 3) + 1, -1.6666666666666665)*pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 2)*(-0.0030000000000000001*t*(1 - x[1])/sqrt(pow(t, 2) + 2) + 0.0030000000000000001*pow(x[0], 2) + 0.0030000000000000001*pow(x[1], 2) + 0.0030000000000000001*pow(1 - x[1], 2) + 0.0030000000000000001) \n",
+      "\n",
+      "rhs_nw =  -2*t*pow(pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 3) + 1, -2.0) + 2.0*(t*(2*x[1] - 2) + sqrt(pow(t, 2) + 2))*pow(pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 3) + 1, -3.0)*(0.0030000000000000001*t*(2*x[1] - 2) + 6*x[1]*(0.001*t + 0.001) + 0.0030000000000000001*sqrt(pow(t, 2) + 2))*pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 2) + 0.66666666666666663*pow(pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 3) + 1, -1.6666666666666665)*pow(0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + pow(x[1], 2) + 1) - 0.001, 2)*(-0.0030000000000000001*t*(1 - x[1])/sqrt(pow(t, 2) + 2) + 0.0030000000000000001*pow(x[0], 2) + 0.0030000000000000001*pow(x[1], 2) + 0.0030000000000000001*pow(1 - x[1], 2) + 0.0030000000000000001) \n",
+      "\n"
      ]
     }
    ],
    "source": [
+    "### domain 1\n",
     "x, y = sym.symbols('x[0], x[1]') # needed by UFL\n",
     "t = sym.symbols('t', positive=True)\n",
     "#f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)  # -Laplace(u)\n",
@@ -140,36 +141,107 @@
     "\n",
     "#dtS1_w = sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)\n",
     "#dtS1_nw = -sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)\n",
-    "dtS1_w = sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)\n",
-    "dtS1_nw = -sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)\n",
-    "print(\"dtS1_w\", dtS1_w, \"\\n\")\n",
-    "print(\"dtS1_nw\", dtS1_nw, \"\\n\")\n",
+    "dtS1_w = porosity[1]*sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)\n",
+    "dtS1_nw = -porosity[1]*sym.diff(S_pc_rel[1](p1_nw - p1_w), t, 1)\n",
+    "print(\"dtS1_w = \", dtS1_w, \"\\n\")\n",
+    "print(\"dtS1_nw = \", dtS1_nw, \"\\n\")\n",
     "\n",
     "#dxdxflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)\n",
     "#dydyflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)\n",
-    "dxdxflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)\n",
-    "dydyflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)\n",
+    "dxdxflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)\n",
+    "dydyflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)\n",
     "\n",
-    "rhs_w = dtS1_w + dxdxflux1_w + dydyflux1_w\n",
-    "print(\"rhs_w\", rhs_w, \"\\n\")\n",
-    "rhs_w = sym.expand(rhs_w)\n",
-    "print(\"rhs_w\", rhs_w, \"\\n\")\n",
+    "rhs1_w = dtS1_w + dxdxflux1_w + dydyflux1_w\n",
+    "rhs1_w = sym.printing.ccode(rhs1_w)\n",
+    "print(\"rhs_w = \", rhs1_w, \"\\n\")\n",
+    "#rhs_w = sym.expand(rhs_w)\n",
+    "#print(\"rhs_w\", rhs_w, \"\\n\")\n",
+    "#rhs_w = sym.collect(rhs_w, x)\n",
+    "#print(\"rhs_w\", rhs_w, \"\\n\")\n",
     "\n",
     "#dxdxflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)\n",
     "#dydyflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)\n",
-    "dxdxflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)\n",
-    "dydyflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)\n",
+    "dxdxflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)\n",
+    "dydyflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-S_pc_rel[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)\n",
     "\n",
+    "rhs1_nw = dtS1_nw + dxdxflux1_nw + dydyflux1_nw\n",
+    "rhs1_nw = sym.printing.ccode(rhs1_nw)\n",
+    "print(\"rhs_nw = \", rhs1_nw, \"\\n\")\n",
     "\n",
     "dxdxflux1_w = sym.printing.ccode(dxdxflux1_w)\n",
     "dxdxflux1_nw = sym.printing.ccode(dxdxflux1_nw)\n",
     "dydyflux1_w = sym.printing.ccode(dydyflux1_w)\n",
     "dydyflux1_nw = sym.printing.ccode(dydyflux1_nw)\n",
     "\n",
-    "print(\"dxdxflux1_w = \", dxdxflux1_w, \"\\n\")\n",
-    "print(\"dydyflux1_w = \", dydyflux1_w, \"\\n\")\n",
-    "print(\"dxdxflux1_nw = \", dxdxflux1_nw, \"\\n\")\n",
-    "print(\"dydyflux1_nw = \", dydyflux1_nw, \"\\n\")"
+    "#print(\"dxdxflux1_w = \", dxdxflux1_w, \"\\n\")\n",
+    "#print(\"dydyflux1_w = \", dydyflux1_w, \"\\n\")\n",
+    "#print(\"dxdxflux1_nw = \", dxdxflux1_nw, \"\\n\")\n",
+    "#print(\"dydyflux1_nw = \", dydyflux1_nw, \"\\n\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "dtS2_w =  -0.833333333333333*((0.001*t*(0.5 - x[1])/(x[0] + 0.1) + 0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + 0.25*x[1] + 1) - 0.001)**6 + 1)**(-1.83333333333333)*(0.001*t*(0.5 - x[1])/(x[0] + 0.1) + 0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + 0.25*x[1] + 1) - 0.001)**5*(-0.006*t*(1 - x[1])/sqrt(t**2 + 2) + 0.006*x[0]**2 + 0.0015*x[1] + 0.006*(0.5 - x[1])/(x[0] + 0.1) + 0.006*(1 - x[1])**2 + 0.006) \n",
+      "\n",
+      "dtS2_nw =  0.833333333333333*((0.001*t*(0.5 - x[1])/(x[0] + 0.1) + 0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + 0.25*x[1] + 1) - 0.001)**6 + 1)**(-1.83333333333333)*(0.001*t*(0.5 - x[1])/(x[0] + 0.1) + 0.001*t*(1 - x[1])**2 - 0.001*(1 - x[1])*sqrt(t**2 + 2) + 0.001*(t + 1)*(x[0]**2 + 0.25*x[1] + 1) - 0.001)**5*(-0.006*t*(1 - x[1])/sqrt(t**2 + 2) + 0.006*x[0]**2 + 0.0015*x[1] + 0.006*(0.5 - x[1])/(x[0] + 0.1) + 0.006*(1 - x[1])**2 + 0.006) \n",
+      "\n",
+      "rhs2_w =  5.0*x[0]*(-t - 1)*(-0.0060000000000000001*t*(0.5 - x[1])/pow(x[0] + 0.10000000000000001, 2) + 12*x[0]*(0.001*t + 0.001))*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -3.5)*pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 5) - 2*(-t - 1)*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -2.5) + 2.5*(-0.25*t - 0.25)*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -3.5)*(0.0060000000000000001*t*(2*x[1] - 2) + 0.0015*t - 0.0060000000000000001*t/(x[0] + 0.10000000000000001) + 0.0060000000000000001*sqrt(pow(t, 2) + 2) + 0.0015)*pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 5) - 0.83333333333333337*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -1.8333333333333335)*pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 5)*(-0.0060000000000000001*t*(1 - x[1])/sqrt(pow(t, 2) + 2) + 0.0060000000000000001*pow(x[0], 2) + 0.0015*x[1] + 0.0060000000000000001*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.0060000000000000001*pow(1 - x[1], 2) + 0.0060000000000000001) \n",
+      "\n",
+      "rhs2_nw =  -1.6666666666666667*t*(0.5 - x[1])*(-0.0060000000000000001*t*(0.5 - x[1])/pow(x[0] + 0.10000000000000001, 2) + 12*x[0]*(0.001*t + 0.001))*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -2.666666666666667)*pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 5)/pow(x[0] + 0.10000000000000001, 2) - 2*t*(0.5 - x[1])*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -1.6666666666666667)/pow(x[0] + 0.10000000000000001, 3) - 2*t*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -1.6666666666666667) + 1.6666666666666667*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -2.666666666666667)*(t*(2*x[1] - 2) - t/(x[0] + 0.10000000000000001) + sqrt(pow(t, 2) + 2))*(0.0060000000000000001*t*(2*x[1] - 2) + 0.0015*t - 0.0060000000000000001*t/(x[0] + 0.10000000000000001) + 0.0060000000000000001*sqrt(pow(t, 2) + 2) + 0.0015)*pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 5) + 0.83333333333333337*pow(pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 6) + 1, -1.8333333333333335)*pow(0.001*t*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.001*t*pow(1 - x[1], 2) - 0.001*(1 - x[1])*sqrt(pow(t, 2) + 2) + 0.001*(t + 1)*(pow(x[0], 2) + 0.25*x[1] + 1) - 0.001, 5)*(-0.0060000000000000001*t*(1 - x[1])/sqrt(pow(t, 2) + 2) + 0.0060000000000000001*pow(x[0], 2) + 0.0015*x[1] + 0.0060000000000000001*(0.5 - x[1])/(x[0] + 0.10000000000000001) + 0.0060000000000000001*pow(1 - x[1], 2) + 0.0060000000000000001) \n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "### domain2\n",
+    "\n",
+    "x, y = sym.symbols('x[0], x[1]') # needed by UFL\n",
+    "t = sym.symbols('t', positive=True)\n",
+    "#f = -sym.diff(u, x, 2) - sym.diff(u, y, 2)  # -Laplace(u)\n",
+    "#f = sym.simplify(f)                         # simplify f\n",
+    "p2_w = 1 - (1+t)*(1 + x**2 + 0.25*y)\n",
+    "p2_nw = t*(1-y)**2 - sym.sqrt(2+t**2)*(1-y) +(0.5-y)*(t/(0.1+x))\n",
+    "\n",
+    "#dtS2_w = sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)\n",
+    "#dtS2_nw = -sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)\n",
+    "dtS2_w = porosity[2]*sym.diff(S_pc_rel[2](p2_nw - p2_w), t, 1)\n",
+    "dtS2_nw = -porosity[2]*sym.diff(S_pc_rel[2](p2_nw - p2_w), t, 1)\n",
+    "print(\"dtS2_w = \", dtS2_w, \"\\n\")\n",
+    "print(\"dtS2_nw = \", dtS2_nw, \"\\n\")\n",
+    "\n",
+    "#dxdxflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, x, 1), x, 1)\n",
+    "#dydyflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, y, 1), y, 1)\n",
+    "dxdxflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_w, x, 1), x, 1)\n",
+    "dydyflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_w, y, 1), y, 1)\n",
+    "\n",
+    "rhs2_w = dtS2_w + dxdxflux2_w + dydyflux2_w\n",
+    "rhs2_w = sym.printing.ccode(rhs2_w)\n",
+    "print(\"rhs2_w = \", rhs2_w, \"\\n\")\n",
+    "#rhs_w = sym.expand(rhs_w)\n",
+    "#print(\"rhs_w\", rhs_w, \"\\n\")\n",
+    "#rhs_w = sym.collect(rhs_w, x)\n",
+    "#print(\"rhs_w\", rhs_w, \"\\n\")\n",
+    "\n",
+    "#dxdxflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, x, 1), x, 1)\n",
+    "#dydyflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, y, 1), y, 1)\n",
+    "dxdxflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_nw, x, 1), x, 1)\n",
+    "dydyflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-S_pc_rel[2](p2_nw - p2_w))*sym.diff(p2_nw, y, 1), y, 1)\n",
+    "\n",
+    "rhs2_nw = dtS2_nw + dxdxflux2_nw + dydyflux2_nw\n",
+    "rhs2_nw = sym.printing.ccode(rhs2_nw)\n",
+    "print(\"rhs2_nw = \", rhs2_nw, \"\\n\")\n",
+    "\n",
+    "#dxdxflux2_w = sym.printing.ccode(dxdxflux2_w)\n",
+    "#dxdxflux2_nw = sym.printing.ccode(dxdxflux2_nw)\n",
+    "#dydyflux2_w = sym.printing.ccode(dydyflux2_w)\n",
+    "#dydyflux2_nw = sym.printing.ccode(dydyflux2_nw)\n"
    ]
   },
   {
@@ -189,21 +261,21 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 2",
+   "display_name": "Python 3",
    "language": "python",
-   "name": "python2"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.16"
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3rc1"
   }
  },
  "nbformat": 4,
diff --git a/TP-TP-patch-test-case/plots/subsequent_errors.tex b/TP-TP-patch-test-case/plots/subsequent_errors.tex
index 2d98427949256aeedd7267e3b01644fce5673c46..5c68a2ce1b8946c0dca47207685e92780859e008 100644
--- a/TP-TP-patch-test-case/plots/subsequent_errors.tex
+++ b/TP-TP-patch-test-case/plots/subsequent_errors.tex
@@ -16,32 +16,32 @@
     \begin{tikzpicture}
 % %     \path[draw,dashed,thick] (3.13,0) -- (3.13,5.2);
         \begin{semilogyaxis}[%
-            width=\textwidth,
+            width=0.5\textwidth,
             title={ Subsequent errors for $t = 0$,  $ h \approx 0.02$, $\tau = 2\cdot 10^{-4}$ },
         % 	    axis lines=left,
         % 	legend style = {draw=none},
-            legend cell align = left,
+%             legend cell align = left,
             xlabel= {iterations},
             ylabel= {subsequent errors},
             xmin= 0,
             xmax= 53,
 %             ymin= 0,
 %             ymax= 0.0003,
-        	grid= both, %major or minor
-            axis line style={-Latex[round]},
+        	grid= major, %both or minor
+% %             axis line style={-Latex[round]},
             legend style={
 %                 anchor=north east,
 %                 at={(1,1)},
-                font=\tiny
+                font=\tiny,
             },
-            legend entries={$\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_1)}$,
-                            $\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_2)}$},
+%             legend entries={$\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_1)}$,
+%                             $\bigl\|p_l^i - p_l^{i-1}\bigr\|_{L^2(\dom_2)}$},
             ] % end of axis options %col sep=comma,
-            \addplot table [x=iteration, y=wetting] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
-            \addplot table [x=iteration, y=wetting] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
-%             \addplot table [col sep=comma] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
-% %             \addplot table [col sep=comma] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%        
+            \addplot table [x=iteration, y=wetting] {../output/subdomain1_subsequent_iteration_error_for_phase_nonwetting_at_time0.csv};%
+            \addplot table [x=iteration, y=wetting] {../output/subdomain2_subsequent_iteration_error_for_phase_nonwetting_at_time0.csv};%
+%             \addplot table [x=iteration, y=nonwetting] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};
+% %             \addplot table [col sep=comma] {../output/subdomain1_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%
+% % %             \addplot table [col sep=comma] {../output/subdomain2_subsequent_iteration_error_for_phase_wetting_at_time0.csv};%        
         \end{semilogyaxis}
     \end{tikzpicture}  
-
 \end{document}