diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 268c558c0479a503c2e26e8bf29c29e6cfa75839..9c34840afc9f7c0e16cc3b92e686f3bdfe3229cb 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -632,23 +632,24 @@ class LDDsimulation(object):
         OUTPUT
         ------------------------------------
         analysing_solution_file_dict   tp.Dict[ind, SolutionFile]   dictionary
-                                  containing the solution file for analysing the
-                                  timestep if analyse_timestep == True, else None.
+                              containing the solution file for analysing the
+                              timestep if analyse_timestep == True, else None.
         """
         # in case no time is given, use the time of the class.
         if time is None:
-            # remember:  this has been set to the target timestep by the solver.
+            # remember: this has been set to the target timestep by the solver.
             # the timestep we are calculating the solution for.
             time = self.t
 
         if analyse_timestep:
             # dictionary holding filename strings for each subdomain to store
-            # iterations in, used to analyse the behaviour of the iterates within
-            # the current time step.
+            # iterations in, used to analyse the behaviour of the iterates
+            # within the current time step.
             solution_over_iteration_within_timestep = dict()
         else:
             solution_over_iteration_within_timestep = None
-        # before the iteration gets started all iteration numbers must be nulled
+        # before the iteration gets started all iteration numbers must be
+        # nulled
         for subdomain_index in self.subdomain.keys():
             prepare_subdomain = mp.Process(
                         target=self.prepare_subdomain,
@@ -713,8 +714,6 @@ class LDDsimulation(object):
                                   the timestep if analyse_timestep == True,
                                   else None.
         """
-
-        # remember:  this has been set to the target timestep by the solver.
         # the timestep we are calculating the solution for.
         time = self.t
         subdomain = self.subdomain[subdomain_index]
@@ -773,6 +772,7 @@ class LDDsimulation(object):
                     f"pressure_prev_timestep[{phase}].vector()[:]" + " =\n" + \
                     f"{subdomain.pressure_prev_timestep[phase].vector()[:]}"
                 print(debugstr5)
+
         # we need to update the previous pressure dictionaries of all
         # interfaces of the subdomain with the current solution, so that the
         # subdomain.calc_gli_term method works properly.
@@ -784,32 +784,31 @@ class LDDsimulation(object):
             subdomain.write_pressure_to_interfaces(previous_iter=True)
             subdomain.calc_gl0_term()
 
-
-    def Lsolver_step(self, subdomain_index: int,#
+    def Lsolver_step(self, subdomain_index: int,
                      analyse_condition: bool = False,
                      debug: bool = False) -> tp.Dict[str, float]:
-        """ L-scheme solver iteration step for an object of class subdomain
+        """L-scheme solver iteration step for an object of class subdomain.
 
         Lsolver_step implements L-scheme solver iteration step for an subdomain
-        with index subdomain_index, i.e. for the model form (already in L-scheme
-        form) stored in the domain patch subdomain.
+        with index subdomain_index, i.e. for the model form (already in
+        L-scheme form) stored in the domain patch subdomain.
 
         # parameters
-        subdomain_index:    int                         subdomain index that defines
-                                                        which subdomain to perform
-                                                        the iteration on.
-        iteration:          int                         iteration number
-        debug:              bool                        flag to toggle weather or not
-                                                        to write out each iteration.
+        subdomain_index:    int                 subdomain index that defines
+                                                which subdomain to perform
+                                                the iteration on.
+        iteration:          int                 iteration number
+        debug:              bool                flag to toggle weather or not
+                                                to write out each iteration.
         """
         subdomain = self.subdomain[subdomain_index]
-        iteration = subdomain.iteration_number
+        # iteration = subdomain.iteration_number
         condition_number = None
         if analyse_condition:
             condition_number = dict()
         for phase in subdomain.has_phases:
             # extract L-scheme form and rhs (without gli term) from subdomain.
-            governing_problem = subdomain.governing_problem(phase = phase)
+            governing_problem = subdomain.governing_problem(phase=phase)
             form = governing_problem['form']
             rhs = governing_problem['rhs']
             # assemble the form and rhs
@@ -825,13 +824,22 @@ class LDDsimulation(object):
             if subdomain.outer_boundary is not None:
                 # apply boundary dirichlet conditions for all outer boundaries.
                 for index, boundary in enumerate(subdomain.outer_boundary):
-                    self.outerBC[subdomain_index][index][phase].apply(form_assembled, rhs_assembled)
+                    self.outerBC[subdomain_index][index][phase].apply(
+                        form_assembled, rhs_assembled
+                        )
 
             linear_solver = subdomain.linear_solver
             if linear_solver is None:
-                df.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled, self.solver, self.preconditioner)
+                df.solve(form_assembled,
+                         subdomain.pressure[phase].vector(),
+                         rhs_assembled, self.solver,
+                         self.preconditioner
+                         )
             else:
-                linear_solver.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled)
+                linear_solver.solve(form_assembled,
+                                    subdomain.pressure[phase].vector(),
+                                    rhs_assembled
+                                    )
 
         return condition_number
 
@@ -897,7 +905,7 @@ class LDDsimulation(object):
                         error_out[subdom_ind].update({phase: None})
             else:
                 # error_out.update({subdom_ind: None})
-                error_outupdate({subdom_ind: None})
+                error_out.update({subdom_ind: None})
         return error_out
 
 
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
index 4e2fbf9c7f8fb1adbfcd0c764aba9ec76e7c4ebf..589a3d964d76ddc7cfd72c0ad003bc81010f1e84 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -43,13 +43,13 @@ max_iter_num = 500
 FEM_Lagrange_degree = 1
 
 # GRID AND MESH STUDY SPECIFICATIONS  #########################################
-mesh_study = True
+mesh_study = False
 resolutions = {
                 # 1: 2e-6,  # h=2
                 # 2: 2e-6,  # h=1.1180
                 # 4: 2e-6,  # h=0.5590
-                8: 2e-6,  # h=0.2814
-                # 16: 8e-6, # h=0.1412
+                # 8: 2e-6,  # h=0.2814
+                16: 1e-5, # h=0.1412
                 # 32: 5e-6,
                 # 64: 2e-6,
                 # 128: 2e-6
@@ -60,7 +60,7 @@ resolutions = {
 #  for each element t_0 in starttimes.
 starttimes = [0.0]
 timestep_size = 0.001
-number_of_timesteps = 5
+number_of_timesteps = 500
 
 # LDD scheme parameters  ######################################################
 
@@ -110,8 +110,8 @@ lambda56_w = 4
 lambda56_nw = 4
 
 include_gravity = True
-debugflag = True
-analyse_condition = True
+debugflag = False
+analyse_condition = False
 
 # I/O CONFIG  #################################################################
 # when number_of_timesteps is high, it might take a long time to write all