diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 5442a2dc34a915b02ff22db0aae2fb376be68e53..79d657b237bc504f49f9a8e2e820781472fcf0b5 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -167,14 +167,6 @@ class LDDsimulation(object):
         # dictionary holding bool variables for each subdomain indicating wether
         # minimal error has been achieved, i.e. the calculation can be considered
         # done or not.
-        self._calc_isdone_for_subdomain = dict()
-        # dictionary holding for each subdomain a dictionary with bool variables
-        # for each phase indicating wether minimal error has been achieved, i.e.
-        # the calculation can be considered done or not. This dictionary is necessary,
-        # because it can happen, that one phase needs severly less iterations than
-        # the other to achieve the error criterion.
-        self._calc_isdone_for_phase = dict()
-
 
     def set_parameters(self,#
                        output_dir: str,#
@@ -376,6 +368,10 @@ class LDDsimulation(object):
         if time is None:
             time = self.t
 
+        # numpy array holding for each subdomain the maximum of the relative error
+        # of subsequent iterations of both phases on the subdomain.
+        # max(subsequent_error) is used as a measure if the solver is done or not.
+        max_subsequent_error = np.zeros(self._number_of_subdomains, dtype=float)
         error_criterion = 1E-8#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
         # in case analyse_timestep is None, solution_over_iteration_within_timestep is None
         solution_over_iteration_within_timestep = self.prepare_LDDsolver(analyse_timestep = analyse_timestep)
@@ -384,105 +380,77 @@ class LDDsimulation(object):
         # gobal stopping criterion for the iteration.
         all_subdomains_are_done = False
         while iteration <= self._max_iter_num and not all_subdomains_are_done:
+            iteration += 1
             # we need to loop over the subdomains and solve an L-scheme type step
             # on each subdomain. here it should be possible to parallelise in
             # the future.
             for sd_index, subdomain in self.subdomain.items():
-                # check if the calculation on subdomain has already be marked as
-                # finished.
-                if not self._calc_isdone_for_subdomain[sd_index]:
-                    # set subdomain iteration to current iteration
-                    subdomain.iteration_number = iteration
-                    # if debug:
-                    #     print(f"On subdomain{sd_index}: Entering iteration {iteration}\n")
-                    #     print(f"subdomain{sd_index}.isRichards = {subdomain.isRichards}.")
-                    #     print(f"subdomain{sd_index}.has_phases = {subdomain.has_phases}.")
-                    # solve the problem on subdomain
-                    self.Lsolver_step(subdomain_index = sd_index,#
-                                      debug = debug,
-                                      )
-                    # bypass error checking for the first iteration.
-                    if iteration > 0:
-                        subsequent_iter_error = dict()
-                        for phase in subdomain.has_phases:
-                            error = df.Function(subdomain.function_space[phase])
-                            error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
-                            error_calculated = df.norm(error, 'L2')
-                            subsequent_iter_error.update(
-                                {phase: error_calculated}
-                                )
+                # set subdomain iteration to current iteration
+                subdomain.iteration_number = iteration
+                # solve the problem on subdomain
+                self.Lsolver_step(subdomain_index = sd_index,#
+                                  debug = debug,
+                                  )
+                # bypass subsequent error calculation for the first iteration, becauses
+                # it is zero anyways.
+                if iteration > 1:
+                    subsequent_iter_error = dict()
+                    for phase in subdomain.has_phases:
+                        error = df.Function(subdomain.function_space[phase])
+                        error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
+                        error_calculated = df.norm(error, 'L2')
+                        subsequent_iter_error.update(
+                            {phase: error_calculated}
+                            )
+
+                    # calculate the maximum over phase of subsequent errors for the
+                    # stopping criterion.
+                    max_subsequent_error[sd_index-1] = max(subsequent_iter_error.values())
 
-                            # set flag to stop solving for that phase if error_criterion is met.
-                            if subsequent_iter_error[phase] < error_criterion:
-                                self._calc_isdone_for_phase[sd_index].update({phase: True})
-
-                            if debug:
-                                print(f"t = {self.t}: subsequent error in iteration {iteration} for phase {phase} is {subsequent_iter_error[phase]}")
-
-                        if analyse_timestep:
-                            # save the calculated L2 errors to csv file for plotting
-                            # with pgfplots
-                            subsequent_error_filename = self.output_dir\
-                                +self.output_filename_parameter_part[sd_index]\
-                                +"subsequent_iteration_errors" +"_at_time"+\
-                                "{number:.{digits}f}".format(number=time, digits=2) +".csv"
-                            self.write_subsequent_errors_to_csv(
-                                filename = subsequent_error_filename, #
-                                subdomain_index = sd_index,
-                                errors = subsequent_iter_error,
-                                time = time)
-
-                        # both phases should be done before we mark the subdomain as
-                        # finished.
-                        total_subsequent_iter_error = max(subsequent_iter_error.values())
-                        # check if error criterion has been met
-                        # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
-                        if debug:
-                            # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
-                            print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}")
-                            #print(f"the minimal distance between any to points in a cell mesh.hmin() on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
-                            print(f"the error criterion is tol = {error_criterion}")
-                        # wetting_done = self._calc_isdone_for_phase[sd_index]['wetting']
-                        # nonwetting_done = self._calc_isdone_for_phase[sd_index]['nonwetting']
-                        if total_subsequent_iter_error < error_criterion:
-                            # if debug:
-                            print(f"subdomain{sd_index} reached error = {total_subsequent_iter_error} < tol = {error_criterion} after {iteration} iterations.")
-                            self._calc_isdone_for_subdomain[sd_index] = True
-
-                    # prepare next iteration
-                    # write the newly calculated solution to the inteface dictionaries
-                    # for communication
-                    subdomain.write_pressure_to_interfaces()
                     if analyse_timestep:
-                        subdomain.write_solution_to_xdmf(#
-                            file = solution_over_iteration_within_timestep[sd_index], #
-                            time = time,#
-                            write_iter_for_fixed_time = True,
-                            )#subsequent_errors = subsequent_errors_over_iteration[sd_index])
-                    for phase in subdomain.has_phases:
-                        subdomain.pressure_prev_iter[phase].assign(
-                            subdomain.pressure[phase]
+                        # save the calculated L2 errors to csv file for plotting
+                        # with pgfplots
+                        subsequent_error_filename = self.output_dir\
+                            +self.output_filename_parameter_part[sd_index]\
+                            +"subsequent_iteration_errors" +"_at_time"+\
+                            "{number:.{digits}f}".format(number=time, digits=2) +".csv"
+                        self.write_subsequent_errors_to_csv(
+                            filename = subsequent_error_filename, #
+                            subdomain_index = sd_index,
+                            errors = subsequent_iter_error,
+                            time = time)
+
+                # prepare next iteration
+                # write the newly calculated solution to the inteface dictionaries
+                # of the subdomain for communication.
+                subdomain.write_pressure_to_interfaces()
+                if analyse_timestep:
+                    subdomain.write_solution_to_xdmf(#
+                        file = solution_over_iteration_within_timestep[sd_index], #
+                        time = time,#
+                        write_iter_for_fixed_time = True,
                         )
 
-                else:
-                    # one subdomain is done. Check if all subdomains are done to
-                    # stop the while loop if needed.
-                    # For this, set
+                for phase in subdomain.has_phases:
+                    subdomain.pressure_prev_iter[phase].assign(
+                        subdomain.pressure[phase]
+                    )
+            # end loop over subdomains
+            ##### stopping criterion for the solver.
+            # only check if error criterion has been met after at least one iteration.
+            for iteration > 1:
+                if debug:
+                    # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
+                    print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}")
+                    #print(f"the minimal distance between any to points in a cell mesh.hmin() on subdomain{sd_index} is h={subdomain.mesh.hmin()}")
+                    print(f"the error criterion is tol = {error_criterion}")
+
+                total_subsequent_error = np.amax(max_subsequent_error)
+                if total_subsequent_error < error_criterion:
+                    # if debug:
+                    print(f"subdomain{sd_index} reached error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.")
                     all_subdomains_are_done = True
-                    # then check if all domains are done. If not, reset
-                    # all_subdomains_are_done to False.
-                    for subdomain, isdone in self._calc_isdone_for_subdomain.items():
-                        if not isdone:
-                            all_subdomains_are_done = False
-                    #TODO check if maybe I still need to do
-                    # subdomain.iteration_number += 1
-                    # or adapt the calc_gli_term method to reflect that one subdomain
-                    # could be done earlier than others.
-
-
-            # you wouldn't be so stupid as to write an infinite loop, would you?
-            iteration += 1
-            # end iteration while loop.
+        # end iteration while loop.
 
     def prepare_LDDsolver(self, time: float = None, #
                   debug: bool = False, #
@@ -492,7 +460,6 @@ class LDDsimulation(object):
         This method executes a few prelininary steps before actually starting
         the L-iteration. Namely,
             - create solution files for analysing time step if analyse_timestep == True
-            - the dictionries self._calc_isdone_for_subdomain are set to False
             - Dirichlet boundary condition terms are evaluated to the current time
             - all subdomain iteration numbers  on interfaces are set to 0.
             - all subdomain iteration numbers subdomain.iteration_number are set
@@ -500,8 +467,6 @@ class LDDsimulation(object):
             - source terms are evaluated to the current time.
             - the calculated solution pressure from the previous timestep
               is assigned to pressure_prev for each subdomain.
-            - The self._calc_isdone_for_phase dictionary for each phase on each
-              subdomain are set to False
             - gl0 terms are calculated.
 
         PARAMETERS
@@ -533,7 +498,6 @@ class LDDsimulation(object):
         else:
             solution_over_iteration_within_timestep = None
         # before the iteration gets started all iteration numbers must be nulled
-        # and the self._calc_isdone_for_subdomain dictionary must be set to False
         for ind, subdomain in self.subdomain.items():
             # create xdmf files for writing out iterations if analyse_timestep is
             # given.
@@ -545,7 +509,6 @@ class LDDsimulation(object):
                     {ind: SolutionFile(self._mpi_communicator, filename)}
                     )
 
-            self._calc_isdone_for_subdomain.update({ind: False})
             # reset all interface[has_interface].current_iteration[ind] = 0, for
             # index has_interface in subdomain.has_interface
             subdomain.null_all_interface_iteration_numbers()
@@ -569,8 +532,6 @@ class LDDsimulation(object):
                 if debug:
                     print(colored(f"id(subdomain{subdomain.subdomain_index}.pressure[{phase}])", "red"), " =\n", f"{id(subdomain.pressure[phase])}")
                     print(f"subdomain{subdomain.subdomain_index}.pressure_prev_timestep[{phase}].vector()[:]", " =\n", f"{subdomain.pressure_prev_timestep[phase].vector()[:]}")
-                # All isdone flags need to be zero before we start iterating
-                self._calc_isdone_for_phase[ind].update({phase: False})
 
             subdomain.calc_gl0_term()
 
@@ -593,65 +554,61 @@ class LDDsimulation(object):
                                                         to write out each iteration.
         """
         subdomain = self.subdomain[subdomain_index]
-        calc_isdone_for_phase = self._calc_isdone_for_phase[subdomain_index]
         iteration = subdomain.iteration_number
         # calculate the gli terms on the right hand side  and store
         subdomain.calc_gli_term()
 
         for phase in subdomain.has_phases:
-            if not calc_isdone_for_phase[phase]:
-                # extract L-scheme form and rhs (without gli term) from subdomain.
-                governing_problem = subdomain.governing_problem(phase = phase)
-                form = governing_problem['form']
-                rhs_without_gli = governing_problem['rhs_without_gli']
-                # assemble the form and rhs
-                form_assembled = df.assemble(form)
-                rhs_without_gli_assembled = df.assemble(rhs_without_gli)
-                # access the assembled gli term for the rhs.
-                gli_term_assembled = subdomain.gli_assembled[phase]
-                # subdomain.calc_gli_term() asslembles gli but on the left hand side
-                # so gli_term_assembled needs to actually be subtracted from the rhs.
-                if debug and subdomain.mesh.num_cells() < 36:
-                    print("\nSystem before applying outer boundary conditions:")
-                    print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
-                rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
-                rhs_assembled = rhs_without_gli_assembled
-
-                if debug and subdomain.mesh.num_cells() < 36:
-                    # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
-                    # print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
-                    print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
-                    print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
-
-                # apply outer Dirichlet boundary conditions if present.
-                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)
-
-                if debug and subdomain.mesh.num_cells() < 36:
-                    print("\nSystem after applying outer boundary conditions:")
-                    # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
-                    print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
-
-                # # set previous iteration as initial guess for the linear solver
-                # pi_prev = subdomain.pressure_prev_iter[phase].vector().get_local()
-                # subdomain.pressure[phase].vector().set_local(pi_prev)
-
-                if debug and subdomain.mesh.num_cells() < 36:
-                    print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local())
-
-                linear_solver = subdomain.linear_solver
-                if linear_solver is None:
-                    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)
-
-                if debug and subdomain.mesh.num_cells() < 36:
-                    print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
+            # if not calc_isdone_for_phase[phase]:
+            # extract L-scheme form and rhs (without gli term) from subdomain.
+            governing_problem = subdomain.governing_problem(phase = phase)
+            form = governing_problem['form']
+            rhs_without_gli = governing_problem['rhs_without_gli']
+            # assemble the form and rhs
+            form_assembled = df.assemble(form)
+            rhs_without_gli_assembled = df.assemble(rhs_without_gli)
+            # access the assembled gli term for the rhs.
+            gli_term_assembled = subdomain.gli_assembled[phase]
+            # subdomain.calc_gli_term() asslembles gli but on the left hand side
+            # so gli_term_assembled needs to actually be subtracted from the rhs.
+            if debug and subdomain.mesh.num_cells() < 36:
+                print("\nSystem before applying outer boundary conditions:")
+                print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+            rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
+            rhs_assembled = rhs_without_gli_assembled
+
+            if debug and subdomain.mesh.num_cells() < 36:
+                # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+                # print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+                print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
+                print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+
+            # apply outer Dirichlet boundary conditions if present.
+            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)
+
+            if debug and subdomain.mesh.num_cells() < 36:
+                print("\nSystem after applying outer boundary conditions:")
+                # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+                print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+
+            # # set previous iteration as initial guess for the linear solver
+            # pi_prev = subdomain.pressure_prev_iter[phase].vector().get_local()
+            # subdomain.pressure[phase].vector().set_local(pi_prev)
+
+            if debug and subdomain.mesh.num_cells() < 36:
+                print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local())
+
+            linear_solver = subdomain.linear_solver
+            if linear_solver is None:
+                df.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled, self.solver, self.preconditioner)
             else:
-                print(f"calculation for phase {phase} is already done in timestep = {self.timestep_num},  t = {self.t}. skipping ..." )
-                pass
+                linear_solver.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled)
+
+            if debug and subdomain.mesh.num_cells() < 36:
+                print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
 
     def _init_solution_files(self):
         """ set up solution files for saving the solution of the LDD scheme for
@@ -962,10 +919,6 @@ class LDDsimulation(object):
             # for parameter, value in self.linear_solver.parameters.items():
             #     print(f"parameter: {parameter} = {value}")
             # df.info(solver_param, True)
-            # setup self._calc_isdone_for_phase dictionary.
-            self._calc_isdone_for_phase.update({subdom_num: dict()})
-
-
 
 
     def _get_interfaces(self, subdomain_index: int, #