diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 0d9b8fa87b03eda44733b4abb3af563d40120f63..e640d20a03367faea7880067cde01fc8f8454ba7 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -1037,17 +1037,9 @@ class LDDsimulation(object):
         for subdom_ind, subdomain in self.subdomain.items():
             file = self.solution_file[subdom_ind]
             for timestep in self.timesteps_to_plot:
+                # write source terms to files
                 source = dict()
-                if self.exact_solution:
-                    exact_pressure = dict()
-                    S = self.saturation[subdom_ind]
-                    saturation_w = df.Function(
-                        subdomain.function_space["pressure"]['wetting']
-                        )
-                    saturation_nw = df.Function(
-                        subdomain.function_space["pressure"]['nonwetting']
-                        )
-
+                exact_pressure = dict()
                 for phase in subdomain.has_phases:
                     f_expr = subdomain.source[phase]
                     f_expr.t = timestep
@@ -1062,6 +1054,8 @@ class LDDsimulation(object):
                         "source_"+"{}".format(phase)
                         )
                     file.write(source[phase], timestep)
+
+                    # write exact solutions if they exist
                     if self.exact_solution:
                         pa_exact = subdomain.pressure_exact[phase]
                         pa_exact.t = timestep
@@ -1076,6 +1070,8 @@ class LDDsimulation(object):
                             "exact_pressure_"+"{}".format(phase)
                             )
                         file.write(exact_pressure[phase], timestep)
+
+                # write pc and S to file if exact solution exists
                 if self.exact_solution:
                     exact_capillary_pressure = df.Function(
                         subdomain.function_space["pressure"]['wetting']
@@ -1089,8 +1085,19 @@ class LDDsimulation(object):
                             exact_pressure["nonwetting"].vector().get_local()\
                             - exact_pressure["wetting"].vector().get_local()
                         exact_capillary_pressure.vector().set_local(pc_temp)
+
                     exact_capillary_pressure.rename("pc_exact", "pc_exact")
                     file.write(exact_capillary_pressure, timestep)
+
+                    S = self.saturation[subdom_ind]
+                    # even if we assume Richards model, we have a nw saturation
+                    # and will write that to file.
+                    saturation_w = df.Function(
+                        subdomain.function_space["pressure"]['wetting']
+                        )
+                    saturation_nw = df.Function(
+                        subdomain.function_space["pressure"]['wetting']
+                        )
                     saturation_w = df.project(
                         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 daaa208a3908825041c8254eec6f8cc8a1c1178e..f6d15ca6edc8e6614d34b05465c3f0f2620a4eed 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-layeredSoilIP-same-intrinsic-control-TPTP-nw-zero"
+use_case = "TP-R-layeredSoilIP-same-intrinsic-after-fixing-exact-writeout"
 # 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 = 50
+number_of_timesteps = 5000
 
 # LDD scheme parameters  ######################################################
 Lw1 = 0.01  # /timestep_size
@@ -114,10 +114,10 @@ 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 = 2
+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
+number_of_timesteps_to_analyse = 6
 
 # 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.
@@ -169,8 +169,8 @@ outer_boundary_def_points = substructuring.outer_boundary_def_points
 
 # MODEL CONFIGURATION #########################################################
 isRichards = {
-    1: False,
-    2: False,
+    1: True,
+    2: True,
     3: False,
     4: False,
     5: False,