use already calculated wetting pressure for nonwetting calculation

LDDsimulation/LDDsimulation.py

@@ -417,7 +417,7 @@ class LDDsimulation(object):
                 # set subdomain iteration to current iteration
                 subdomain.iteration_number = iteration
                 # solve the problem on subdomain
-                self.Lsolver_step(subdomain_index = sd_index,#
+                subsequent_error = self.Lsolver_step(subdomain_index = sd_index,#
                                   debug = debug,
                                   )
                 # bypass subsequent error calculation for the first iteration, becauses
@@ -425,11 +425,11 @@ class LDDsimulation(object):
                 if iteration > 1:
                     subsequent_iter_error = dict()
                     for phase in subdomain.has_phases:
-                        error = df.Function(subdomain.function_space["pressure"][phase])
-                        error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
-                        error_calculated = df.norm(error, 'L2')
+                        # error = df.Function(subdomain.function_space["pressure"][phase])
+                        # error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
+                        # error_calculated = df.norm(error, 'L2')
                         subsequent_iter_error.update(
-                            {phase: error_calculated}
+                            {phase: subsequent_error[phase]}
                             )
                         if debug:
                             print(f"time = {time}: subsequent error on subdomain {sd_index} for {phase} phase in iteration {iteration} = ", subsequent_iter_error[phase])
@@ -462,10 +462,10 @@ class LDDsimulation(object):
                         write_iter_for_fixed_time = True,
                         )
 
-                for phase in subdomain.has_phases:
-                    subdomain.pressure_prev_iter[phase].assign(
-                        subdomain.pressure[phase]
-                    )
+                # 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.
@@ -579,7 +579,7 @@ class LDDsimulation(object):
         return solution_over_iteration_within_timestep
 
     def Lsolver_step(self, subdomain_index: int,#
-                     debug: bool = False) -> None:
+                     debug: bool = False) -> float:
         """ L-scheme solver iteration step for an object of class subdomain
 
         Lsolver_step implements L-scheme solver iteration step for an subdomain
@@ -596,7 +596,7 @@ class LDDsimulation(object):
         """
         subdomain = self.subdomain[subdomain_index]
         iteration = subdomain.iteration_number
-
+        subsequent_error = 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)
@@ -633,6 +633,16 @@ class LDDsimulation(object):
             if debug and subdomain.mesh.num_cells() < 36:
                 print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
 
+            error = df.Function(subdomain.function_space["pressure"][phase])
+            error.assign(subdomain.pressure[phase] - subdomain.pressure_prev_iter[phase])
+            subsequent_error.update({phase: df.norm(error, 'L2')})
+
+            subdomain.pressure_prev_iter[phase].assign(
+                subdomain.pressure[phase]
+            )
+
+        return subsequent_error
+
 
     def _init_solution_files(self):
         """ set up solution files for saving the solution of the LDD scheme for