diff --git a/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py
index a9195f72a7c45853439c87ec7e0c503fdbb83eb0..17b8f55aaf8709b89a4d8181efe3087f0af02eaf 100755
--- a/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py
+++ b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py
@@ -20,17 +20,24 @@ datestr = date.strftime("%Y-%m-%d")
 sym.init_printing()
 
 use_case = "TP-TP-2-patch-pure-dd"
-solver_tol = 6E-7
+# solver_tol = 6E-7
 max_iter_num = 1000
 FEM_Lagrange_degree = 1
 mesh_study = True
-resolutions = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
+resolutions = { 1: 1e-7,
+                2: 2e-5,
+                4: 1e-6,
+                8: 1e-6,
+                16: 5e-7,
+                32: 5e-7,
+                64: 5e-7,
+                128: 5e-7}
 
 ############ GRID #######################
 # mesh_resolution = 20
-timestep_size = 0.0001
-number_of_timesteps = 200
-plot_timestep_every = 1
+timestep_size = 0.000025
+number_of_timesteps = 4000
+plot_timestep_every = 40
 # decide how many timesteps you want analysed. Analysed means, that we write out
 # subsequent errors of the L-iteration within the timestep.
 number_of_timesteps_to_analyse = 0
@@ -46,7 +53,7 @@ include_gravity = False
 debugflag = False
 analyse_condition = False
 
-output_string = "./output/{}-{}_timesteps{}_P{}-solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)
+output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
 
 # toggle what should be written to files
 if mesh_study:
@@ -537,7 +544,7 @@ for subdomain in isRichards.keys():
 #
 # sa
 
-for mesh_resolution in resolutions:
+for mesh_resolution, solver_tol in resolutions.items():
     # initialise LDD simulation class
     simulation = ldd.LDDsimulation(
         tol=1E-14,
diff --git a/TP-TP-layered-soil-case/TP-TP-layered_soil.py b/TP-TP-layered-soil-case/TP-TP-layered_soil.py
index d787da758767cf07f00bc2f5af70a092f824e624..5de8dbc47ef0a00759dc00faff91b2201bbb2f09 100755
--- a/TP-TP-layered-soil-case/TP-TP-layered_soil.py
+++ b/TP-TP-layered-soil-case/TP-TP-layered_soil.py
@@ -17,33 +17,80 @@ import functools as ft
 import domainPatch as dp
 import LDDsimulation as ldd
 import helpers as hlp
+import datetime
+import os
+import pandas as pd
+
+date = datetime.datetime.now()
+datestr = date.strftime("%Y-%m-%d")
 
 # init sympy session
 sym.init_printing()
-
-use_case="TP-TP-layered-soil"
-solver_tol = 1E-6
-
-############ GRID #######################ΓΌ
-mesh_resolution = 30
-timestep_size = 0.0005
+# solver_tol = 6E-7
+use_case = "TP-TP-layered-soil"
+max_iter_num = 1000
+FEM_Lagrange_degree = 1
+mesh_study = False
+resolutions = {
+                # 1: 1e-7,  # h=2
+                # 2: 2e-5,  # h=1.1180
+                # 4: 1e-6,  # h=0.5590
+                # 8: 1e-6,  # h=0.2814
+                # 16: 5e-7, # h=0.1412
+                32: 5e-7,
+                # 64: 5e-7,
+                # 128: 5e-7
+                }
+
+############ GRID #######################
+# mesh_resolution = 20
+timestep_size = 0.000025
 number_of_timesteps = 20
+plot_timestep_every = 100
 # decide how many timesteps you want analysed. Analysed means, that we write out
 # subsequent errors of the L-iteration within the timestep.
-number_of_timesteps_to_analyse = 5
-starttime = 0
+number_of_timesteps_to_analyse = 0
+starttime = 0.0
 
-Lw = 0.25  #/timestep_size
+Lw = 0.25 #/timestep_size
 Lnw=Lw
 
-l_param_w = 40
-l_param_nw = 40
+lambda_w = 40
+lambda_nw = 40
 
 include_gravity = True
 debugflag = False
 analyse_condition = False
 
-output_string = "./output/test-after-bugfix-number_of_timesteps{}_".format(number_of_timesteps)
+if mesh_study:
+    output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
+else:
+    for tol in resolutions.values():
+        solver_tol = tol
+    output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)
+
+
+# toggle what should be written to files
+if mesh_study:
+    write_to_file = {
+        'space_errornorms': True,
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': False,
+        'solutions': False,
+        'absolute_differences': False,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': False
+    }
+else:
+    write_to_file = {
+        'space_errornorms': True,
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': False,
+        'solutions': True,
+        'absolute_differences': True,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': True
+    }
 
 # global domain
 subdomain0_vertices = [df.Point(-1.0,-1.0), #
@@ -240,14 +287,14 @@ L = {
 
 # subdom_num : lambda parameter for the L-scheme
 lambda_param = {
-    1: {'wetting': l_param_w,
-         'nonwetting': l_param_nw},#
-    2: {'wetting': l_param_w,
-         'nonwetting': l_param_nw},#
-    3: {'wetting': l_param_w,
-         'nonwetting': l_param_nw},#
-    4: {'wetting': l_param_w,
-         'nonwetting': l_param_nw},#
+    1: {'wetting': lambda_w,
+         'nonwetting': lambda_nw},#
+    2: {'wetting': lambda_w,
+         'nonwetting': lambda_nw},#
+    3: {'wetting': lambda_w,
+         'nonwetting': lambda_nw},#
+    4: {'wetting': lambda_w,
+         'nonwetting': lambda_nw},#
 }
 
 
@@ -265,12 +312,12 @@ def rel_perm1nw(s):
 ## relative permeabilty functions on subdomain 2
 def rel_perm2w(s):
     # relative permeabilty wetting on subdomain2
-    return s**2
+    return s**3
 
 
 def rel_perm2nw(s):
     # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
-    return (1-s)**2
+    return (1-s)**3
 
 
 _rel_perm1w = ft.partial(rel_perm1w)
@@ -316,11 +363,11 @@ def rel_perm1nw_prime(s):
 # relative permeabilty functions on subdomain 1
 def rel_perm2w_prime(s):
     # relative permeabilty on subdomain1
-    return 2*s
+    return 3*s**2
 
 def rel_perm2nw_prime(s):
     # relative permeabilty on subdomain1
-    return -2*(1-s)
+    return -3*(1-s)**2
 
 _rel_perm1w_prime = ft.partial(rel_perm1w_prime)
 _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -380,22 +427,22 @@ def saturation_sym_prime(pc, n_index, alpha):
 S_pc_sym = {
     1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
     2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-    3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-    4: ft.partial(saturation_sym, n_index=3, alpha=0.001)
+    3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
+    4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
 }
 
 S_pc_sym_prime = {
     1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
     2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-    3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-    4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
+    3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
+    4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
 }
 
 sat_pressure_relationship = {
     1: ft.partial(saturation, n_index=3, alpha=0.001),
     2: ft.partial(saturation, n_index=3, alpha=0.001),
-    3: ft.partial(saturation, n_index=3, alpha=0.001),
-    4: ft.partial(saturation, n_index=3, alpha=0.001)
+    3: ft.partial(saturation, n_index=6, alpha=0.001),
+    4: ft.partial(saturation, n_index=6, alpha=0.001)
 }
 
 
@@ -409,7 +456,7 @@ p_e_sym_2patch = {
     1: {'wetting': -7 - (1+t*t)*(1 + x*x + y*y),
         'nonwetting': -1-t*(1.1 + y + x**2)**2},
     2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x),
-        'nonwetting': -1-t*(1.1 + x**2)**2 - sym.sqrt(2+t**2)*y**2},
+        'nonwetting': -1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2)*y**2},
 }
 
 p_e_sym = {
@@ -496,42 +543,71 @@ for subdomain in isRichards.keys():
                 {outer_boundary_ind: exact_solution[subdomain]}
                 )
 
-write_to_file = {
-    'meshes_and_markers': True,
-    'L_iterations': True
-}
-
-# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=debugflag, LDDsolver_tol=solver_tol)
-simulation.set_parameters(use_case=use_case,
-                          output_dir=output_string,
-                          subdomain_def_points=subdomain_def_points,
-                          isRichards=isRichards,
-                          interface_def_points=interface_def_points,
-                          outer_boundary_def_points=outer_boundary_def_points,
-                          adjacent_subdomains=adjacent_subdomains,
-                          mesh_resolution=mesh_resolution,
-                          viscosity=viscosity,
-                          porosity=porosity,
-                          L=L,
-                          lambda_param=lambda_param,
-                          relative_permeability=relative_permeability,
-                          saturation=sat_pressure_relationship,
-                          starttime=starttime,
-                          number_of_timesteps=number_of_timesteps,
-                          number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
-                          timestep_size=timestep_size,
-                          sources=source_expression,
-                          initial_conditions=initial_condition,
-                          dirichletBC_expression_strings=dirichletBC,
-                          exact_solution=exact_solution,
-                          densities=densities,
-                          include_gravity=include_gravity,
-                          write2file=write_to_file,
-                          )
-
-simulation.initialise()
-# print(simulation.__dict__)
-simulation.run(analyse_condition=analyse_condition)
-# simulation.LDDsolver(time=0, debug=True, analyse_timestep=True)
-# df.info(parameters, True)
+for mesh_resolution, solver_tol in resolutions.items():
+    # initialise LDD simulation class
+    simulation = ldd.LDDsimulation(
+        tol=1E-14,
+        LDDsolver_tol=solver_tol,
+        debug=debugflag,
+        max_iter_num=max_iter_num,
+        FEM_Lagrange_degree=FEM_Lagrange_degree,
+        mesh_study=mesh_study
+        )
+
+    simulation.set_parameters(use_case=use_case,
+                              output_dir=output_string,
+                              subdomain_def_points=subdomain_def_points,
+                              isRichards=isRichards,
+                              interface_def_points=interface_def_points,
+                              outer_boundary_def_points=outer_boundary_def_points,
+                              adjacent_subdomains=adjacent_subdomains,
+                              mesh_resolution=mesh_resolution,
+                              viscosity=viscosity,
+                              porosity=porosity,
+                              L=L,
+                              lambda_param=lambda_param,
+                              relative_permeability=relative_permeability,
+                              saturation=sat_pressure_relationship,
+                              starttime=starttime,
+                              number_of_timesteps=number_of_timesteps,
+                              number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
+                              plot_timestep_every=plot_timestep_every,
+                              timestep_size=timestep_size,
+                              sources=source_expression,
+                              initial_conditions=initial_condition,
+                              dirichletBC_expression_strings=dirichletBC,
+                              exact_solution=exact_solution,
+                              densities=densities,
+                              include_gravity=include_gravity,
+                              write2file=write_to_file,
+                              )
+
+    simulation.initialise()
+    output_dir = simulation.output_dir
+    # simulation.write_exact_solution_to_xdmf()
+    output = simulation.run(analyse_condition=analyse_condition)
+    for subdomain_index, subdomain_output in output.items():
+        mesh_h = subdomain_output['mesh_size']
+        for phase, different_errornorms in subdomain_output['errornorm'].items():
+            filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
+            # for errortype, errornorm in different_errornorms.items():
+
+                # eocfile = open("eoc_filename", "a")
+                # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
+                # eocfile.close()
+                # if subdomain.isRichards:mesh_h
+            data_dict = {
+                'mesh_parameter': mesh_resolution,
+                'mesh_h': mesh_h,
+            }
+            for error_type, errornorms in different_errornorms.items():
+                data_dict.update(
+                    {error_type: errornorms}
+                )
+            errors = pd.DataFrame(data_dict, index=[mesh_resolution])
+            # check if file exists
+            if os.path.isfile(filename) == True:
+                with open(filename, 'a') as f:
+                    errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+            else:
+                errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)