diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
index 1feb5db69e127fd9984744671b3c4532dd462a38..6afe2dc92a9745702ef31231bf3c53686ea76a5f 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
@@ -1,28 +1,45 @@
 #!/usr/bin/python3
+"""TPR 2 patch soil simulation.
+
+This program sets up an LDD simulation
+"""
+
 import dolfin as df
-import mshr
-import numpy as np
 import sympy as sym
-import typing as tp
-import domainPatch as dp
-import LDDsimulation as ldd
 import functools as ft
+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")
-#import ufl as ufl
-
 # init sympy session
 sym.init_printing()
 
+# PREREQUISITS  ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
+if not os.path.exists('./output'):
+    os.mkdir('./output')
+    print("Directory ", './output',  " created ")
+else:
+    print("Directory ", './output',  " already exists. Will use as output \
+    directory")
+
+date = datetime.datetime.now()
+datestr = date.strftime("%Y-%m-%d")
+
+# Name of the usecase that will be printed during simulation.
 use_case = "TPR-2-patch-realistic-testrun"
-# solver_tol = 6E-7
+# The name of this very file. Needed for creating log output.
+thisfile = "TP-R-2-patch-mesh-study.py"
+
+# GENERAL SOLVER CONFIG  ######################################################
+# maximal iteration per timestep
 max_iter_num = 500
 FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS  #########################################
 mesh_study = True
 resolutions = {
                 1: 1e-6,
@@ -36,18 +53,19 @@ resolutions = {
                 256: 1e-6,
                 }
 
-############ GRID #######################
-# mesh_resolution = 20
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+#  for each element t_0 in starttimes.
+starttimes = [0.0]
 timestep_size = 0.001
 number_of_timesteps = 1000
-plot_timestep_every = 4
-# 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
-starttimes = [0.0]
 
-Lw = 0.025 #/timestep_size
-Lnw= 0.025
+# LDD scheme parameters  ######################################################
+Lw1 = 0.025 #/timestep_size
+Lnw1= 0.025
+
+Lw2 = 0.025 #/timestep_size
+Lnw2= 0.025
 
 lambda_w = 40
 lambda_nw = 40
@@ -56,22 +74,38 @@ include_gravity = False
 debugflag = False
 analyse_condition = True
 
-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)
+# I/O CONFIG  #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
 
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
 if mesh_study:
     write_to_file = {
+        # output the relative errornorm (integration in space) w.r.t. an exact
+        # solution for each timestep into a csv file.
         'space_errornorms': True,
+        # save the mesh and marker functions to disk
         'meshes_and_markers': True,
-        'L_iterations_per_timestep': True,
+        # save xdmf/h5 data for each LDD iteration for timesteps determined by
+        # number_of_timesteps_to_analyse. I/O intensive!
+        'L_iterations_per_timestep': False,
+        # save solution to xdmf/h5.
         'solutions': True,
+        # save absolute differences w.r.t an exact solution to xdmf/h5 file
+        # to monitor where on the domains errors happen
         'absolute_differences': True,
+        # analyise condition numbers for timesteps determined by
+        # number_of_timesteps_to_analyse and save them over time to csv.
         'condition_numbers': analyse_condition,
+        # output subsequent iteration errors measured in L^2  to csv for
+        # timesteps determined by number_of_timesteps_to_analyse.
+        # Usefull to monitor convergence of the acutal LDD solver.
         'subsequent_errors': True
     }
 else:
@@ -85,9 +119,20 @@ else:
         'subsequent_errors': True
     }
 
+# OUTPUT FILE STRING  #########################################################
+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
+        )
 
 
-##### Domain and Interface ####
+# DOMAIN AND INTERFACE  #######################################################
 # global simulation domain domain
 sub_domain0_vertices = [df.Point(-1.0, -1.0),
                         df.Point(1.0, -1.0),
@@ -181,28 +226,19 @@ gravity_acceleration = 9.81
 
 L = {#
 # subdom_num : subdomain L for L-scheme
-    1 : {'wetting' :Lw,
-         'nonwetting': Lnw},#
-    2 : {'wetting' :Lw,
-         'nonwetting': Lnw}
+    1 : {'wetting' :Lw1,
+         'nonwetting': Lnw1},#
+    2 : {'wetting' :Lw2,
+         'nonwetting': Lnw2}
 }
 
 
 lambda_param = {#
 # subdom_num : lambda parameter for the L-scheme
-    1 : {'wetting' :lambda_w,
+    0 : {'wetting' :lambda_w,
          'nonwetting': lambda_nw},#
-    2 : {'wetting' :lambda_w,
-         'nonwetting': lambda_nw}
 }
 
-# intrinsic_permeability = {
-#     1: {"wetting": 1,
-#         "nonwetting": 1},
-#     2: {"wetting": 1,
-#         "nonwetting": 1},
-# }
-
 intrinsic_permeability = {
     1: 1,
     2: 1,
@@ -229,7 +265,7 @@ def rel_perm2w(s):
     # relative permeabilty wetting on subdomain2
     return intrinsic_permeability[2]*s**3
 def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
+    # relative permeabilty nonwetting on subdomain2
     return intrinsic_permeability[2]*(1-s)**3
 
 _rel_perm2w = ft.partial(rel_perm2w)
@@ -442,7 +478,7 @@ dirichletBC = dict()
 
 # subdomain index: {outer boudary part index: {phase: expression}}
 for subdomain in isRichards.keys():
-    # if subdomain has no outer boundary, outer_boundary_def_points[subdomain] is None
+    # subdomain can have no outer boundary
     if outer_boundary_def_points[subdomain] is None:
         dirichletBC.update({subdomain: None})
     else:
@@ -455,13 +491,9 @@ for subdomain in isRichards.keys():
                 )
 
 
-# def saturation(pressure, subdomain_index):
-#     # inverse capillary pressure-saturation-relationship
-#     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
-#
-# sa
-
-f = open('TP-R-2-patch-mesh-study.py', 'r')
+# read this file and print it to std out. This way the simulation can produce a
+# log file with ./TP-R-layered_soil.py | tee simulation.log
+f = open(thisfile, 'r')
 print(f.read())
 f.close()
 
@@ -478,33 +510,34 @@ for starttime in starttimes:
             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.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
@@ -512,26 +545,39 @@ for starttime in starttimes:
         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
+            for phase, error_dict in subdomain_output['errornorm'].items():
+                filename = output_dir \
+                    + "subdomain{}".format(subdomain_index)\
+                    + "-space-time-errornorm-{}-phase.csv".format(phase)
+                # for errortype, errornorm in error_dict.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():
+                for norm_type, errornorm in error_dict.items():
                     data_dict.update(
-                        {error_type: errornorms}
+                        {norm_type: errornorm}
                     )
                 errors = pd.DataFrame(data_dict, index=[mesh_resolution])
                 # check if file exists
-                if os.path.isfile(filename) == True:
+                if os.path.isfile(filename) is True:
                     with open(filename, 'a') as f:
-                        errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
+                        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)
+                    errors.to_csv(
+                        filename,
+                        sep='\t',
+                        encoding='utf-8',
+                        index=False
+                        )