diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 421643d638f656f2b5fd2b37d0969fcba54ef550..c83752a31382eb7a88d228a035a9a6f8b839d053 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -315,8 +315,8 @@ class LDDsimulation(object):
         self._init_exact_solution_expression()
         self._init_DirichletBC_dictionary()
         self._init_simulation_output()
-        if self.exact_solution and self.write2file['solutions']:
-            df.info(colored("writing exact solution for all time steps to xdmf files ...", "yellow"))
+        if self.write2file['solutions']:
+            df.info(colored("writing exact solutions and/or source terms for all time steps to xdmf files ...", "yellow"))
             self.write_exact_solution_to_xdmf()
         self._initialised = True
 
@@ -903,44 +903,48 @@ class LDDsimulation(object):
         for subdom_ind, subdomain in self.subdomain.items():
             file = self.solution_file[subdom_ind]
             for timestep in self.timesteps_to_plot:
-                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"]['wetting'])
                 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"]['wetting'])
+
                 for phase in subdomain.has_phases:
                     f_expr = subdomain.source[phase]
-                    pa_exact = subdomain.pressure_exact[phase]
-                    pa_exact.t = timestep
                     f_expr.t = timestep
-                    exact_pressure.update(
-                        {phase: df.project(pa_exact, subdomain.function_space["pressure"][phase])}
-                        )
                     source.update(
                         {phase: df.project(f_expr, subdomain.function_space["pressure"][phase])}
                     )
-                    exact_pressure[phase].rename("exact_pressure_"+"{}".format(phase), "exact_pressure_"+"{}".format(phase))
-                    file.write(exact_pressure[phase], timestep)
                     source[phase].rename("source_"+"{}".format(phase), "source_"+"{}".format(phase))
                     file.write(source[phase], timestep)
 
-                exact_capillary_pressure = df.Function(subdomain.function_space["pressure"]['wetting'])
-                if subdomain.isRichards:
-                    exact_capillary_pressure.assign(-exact_pressure['wetting'])
-                else:
-                    pc_temp = 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)
-                saturation_w = df.project(S(exact_capillary_pressure), subdomain.function_space["pressure"]["wetting"])
-                # saturation_w.assign(Sat_w)
-                saturation_w.rename("Sw", "Sw")
-                file.write(saturation_w, timestep)
-                # S_nw = 1-S(exact_capillary_pressure).vector().get_local()
-                saturation_nw = df.project(1-S(exact_capillary_pressure), subdomain.function_space["pressure"]["wetting"])
-                # saturation_nw.assign(S_nw)
-                saturation_nw.rename("Snw", "Snw")
-                file.write(saturation_nw, timestep)
+                    if self.exact_solution:
+                        pa_exact = subdomain.pressure_exact[phase]
+                        pa_exact.t = timestep
+                        exact_pressure.update(
+                            {phase: df.project(pa_exact, subdomain.function_space["pressure"][phase])}
+                            )
+                        exact_pressure[phase].rename("exact_pressure_"+"{}".format(phase), "exact_pressure_"+"{}".format(phase))
+                        file.write(exact_pressure[phase], timestep)
+                if self.exact_solution:
+                    exact_capillary_pressure = df.Function(subdomain.function_space["pressure"]['wetting'])
+                    if subdomain.isRichards:
+                        exact_capillary_pressure.assign(-exact_pressure['wetting'])
+                    else:
+                        pc_temp = 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)
+                    saturation_w = df.project(S(exact_capillary_pressure), subdomain.function_space["pressure"]["wetting"])
+                    # saturation_w.assign(Sat_w)
+                    saturation_w.rename("Sw", "Sw")
+                    file.write(saturation_w, timestep)
+                    # S_nw = 1-S(exact_capillary_pressure).vector().get_local()
+                    saturation_nw = df.project(1-S(exact_capillary_pressure), subdomain.function_space["pressure"]["wetting"])
+                    # saturation_nw.assign(S_nw)
+                    saturation_nw.rename("Snw", "Snw")
+                    file.write(saturation_nw, timestep)
 
 
     def write_subsequent_errors_to_csv(self,#