diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index c33fcf0d2e2a4b0e51565d4022b14b6b78a97df2..984736a9ad3786f154f378f1d32302d694e8fdf3 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -875,21 +875,6 @@ class LDDsimulation(object):
                     pa_exact = subdomain.pressure_exact[phase]
                     pa_exact.t = self.t
                     error_calculated = df.errornorm(pa_exact, subdomain.pressure[phase], 'L2', degree_rise=4)
-                    pressure_exact.update(
-                        {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
-                        )
-                    absolute_difference = df.Function(subdomain.function_space["pressure"][phase])
-                    pressure_difference = pressure_exact[phase].vector()[:] - subdomain.pressure[phase].vector()[:]
-                    abs_diff_tmp = np.fabs(pressure_difference)
-                    absolute_difference.vector()[:] = abs_diff_tmp
-                    # dx = subdomain.dx
-                    # error_calculated_L1 = df.assemble(absolute_difference*dx)
-                    # error_calculated_L2 = np.sqrt(df.assemble(absolute_difference**2*dx))
-                    # error_calculated_L2_2 = df.norm(absolute_difference, norm_type='L2', mesh=subdomain.mesh)
-                    # print(f"Errornorm dolfin: {error_calculated}")
-                    # print(f"Errornorm manually calculated L1: {error_calculated_L1}")
-                    # print(f"Errornorm manually calculated L2: {error_calculated_L2}")
-                    # print(f"Errornorm manually calculated L2 with df.norm: {error_calculated_L2_2}")
                     self.output[subdom_ind]['errornorm'][phase]['L1'] += self.timestep_size*error_calculated
                     self.output[subdom_ind]['errornorm'][phase]['L2'] += self.timestep_size*error_calculated**2
                     # print(f"Linf error on subdomain {subdom_ind} and phase {phase} before checking: {self.output[subdom_ind]['errornorm'][phase]['Linf']}")
@@ -1067,7 +1052,7 @@ class LDDsimulation(object):
                     f_expr = subdomain.source[phase]
                     f_expr.t = timestep
                     source.update(
-                        {phase: df.project(
+                        {phase: df.interpolate(
                             f_expr, subdomain.function_space["pressure"][phase]
                             )
                         }
@@ -1082,7 +1067,7 @@ class LDDsimulation(object):
                         pa_exact = subdomain.pressure_exact[phase]
                         pa_exact.t = timestep
                         exact_pressure.update(
-                            {phase: df.project(
+                            {phase: df.interpolate(
                                 pa_exact,
                                 subdomain.function_space["pressure"][phase]
                                 )}
@@ -1107,7 +1092,7 @@ class LDDsimulation(object):
                         exact_capillary_pressure.vector().set_local(pc_temp)
                     exact_capillary_pressure.rename("pc_exact", "pc_exact")
                     file.write(exact_capillary_pressure, timestep)
-                    saturation_w = df.project(
+                    saturation_w = df.interpolate(
                         S(exact_capillary_pressure),
                         subdomain.function_space["pressure"]["wetting"]
                         )
@@ -1115,7 +1100,7 @@ class LDDsimulation(object):
                     saturation_w.rename("Sw_exact", "Sw_exact")
                     file.write(saturation_w, timestep)
                     # S_nw = 1-S(exact_capillary_pressure).vector().get_local()
-                    saturation_nw = df.project(
+                    saturation_nw = df.interpolate(
                         1-S(exact_capillary_pressure),
                         subdomain.function_space["pressure"]["wetting"]
                         )
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 c7e5397553f8e240f885a1848c97bb22e9707442..a59d77b805dbc2c19accdb1b66d7e5bf72668033 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
@@ -30,7 +30,7 @@ date = datetime.datetime.now()
 datestr = date.strftime("%Y-%m-%d")
 
 # Name of the usecase that will be printed during simulation.
-use_case = "TP-R-layered-soil-with-inner-patch-realistic-same-intrinsic-longterm"
+use_case = "TP-R-layeredSoilIP-same-intrinsic-plot-every-timestep-test"
 # The name of this very file. Needed for creating log output.
 thisfile = "TP-R-layered_soil_with_inner_patch-realistic.py"
 
@@ -58,7 +58,7 @@ resolutions = {
 starttimes = {0: 0.0}
 # starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
 timestep_size = 0.001
-number_of_timesteps = 5000
+number_of_timesteps = 50
 
 # LDD scheme parameters  ######################################################
 Lw1 = 0.01  # /timestep_size
@@ -114,7 +114,7 @@ analyse_condition = False
 # 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 = 10
+plot_timestep_every = 1
 # 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