clean _eval_sources(). check ids

check that source terms have different python ids.

LDDsimulation/LDDsimulation.py

@@ -100,9 +100,9 @@ class LDDsimulation(object):
 
         ## Private variables
         # maximal number of L-iterations that the LDD solver uses.
-        self._max_iter_num = 5000
+        self._max_iter_num = 1
         # TODO rewrite this with regard to the mesh sizes
-        self.calc_tol = self.tol
+        # self.calc_tol = self.tol
         # list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps
         self.analyse_timesteps = None
         # The number of subdomains are counted by self.init_meshes_and_markers()
@@ -143,14 +143,14 @@ class LDDsimulation(object):
         ########## SOLVER DEFINITION AND PARAMETERS ########
         # print("\nLinear Algebra Backends:")
         # df.list_linear_algebra_backends()
-        # print("\nLinear Solver Methods:")
+        # # print("\nLinear Solver Methods:")
         # df.list_linear_solver_methods()
-        # print("\nPeconditioners for Krylov Solvers:")
+        # # print("\nPeconditioners for Krylov Solvers:")
         # df.list_krylov_solver_preconditioners()
 
         ### Define the linear solver to be used.
         self.solver = 'bicgstab' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
-        self.preconditioner = 'hypre_amg' #'ilu'#'hypre_amg' # incomplete LU factorization
+        self.preconditioner = 'jacobi'#'hypre_amg' #'ilu'#'hypre_amg' # incomplete LU factorization
         # dictionary of solver parametrs. This is passed to self._init_subdomains,
         # where for each subdomain a sovler object of type self.solver is created
         # with these parameters.
@@ -161,6 +161,17 @@ class LDDsimulation(object):
             'maximum_iterations': 1000,
             'report': False
         }
+        # Dictionaries needed by the LDD Solver
+        # 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,#
@@ -285,9 +296,9 @@ class LDDsimulation(object):
             # check if the solver should beanalised or not.
             if analyse_solver_at_timesteps is not None and np.isin(self.timestep_num, analyse_solver_at_timesteps, assume_unique=True):
                 print("analyising timestep ...")
-                self.LDDsolver(debug = False, analyse_timestep = True)
+                self.LDDsolver(debug = True, analyse_timestep = True)
             else:
-                self.LDDsolver(debug = False, analyse_timestep = False)
+                self.LDDsolver(debug = True, analyse_timestep = False)
 
             self.t += self.timestep_size
             self.timestep_num += 1
@@ -347,18 +358,16 @@ class LDDsimulation(object):
         if time is None:
             time = self.t
 
+        error_criterion = 1E-8#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
+
         if analyse_timestep:
             # dictionary holding filename strings for each subdomain to store
             # iterations in, used to analyse the behaviour of the iterates within
             # the current time step.
             solution_over_iteration_within_timestep = dict()
             subsequent_errors_over_iteration = dict()
-        # dictionary holding bool variables for each subdomain indicating wether
-        # minimal error has been achieved, i.e. the calculation can be considered
-        # done or not.
-        subdomain_calc_isdone = dict()
         # before the iteration gets started all iteration numbers must be nulled
-        # and the subdomain_calc_isdone dictionary must be set to False
+        # 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.
@@ -375,7 +384,7 @@ class LDDsimulation(object):
                 #     subsequent_errors_over_iteration[ind].update( {phase: dict()} )
                 #
 
-            subdomain_calc_isdone.update({ind: False})
+            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()
@@ -386,18 +395,16 @@ class LDDsimulation(object):
             self.update_DirichletBC_dictionary(time = time, #
                 subdomain_index = ind,
                 )
-            # the following only needs to be done when time is not equal 0 since
-            # our initialisation functions already took care of setting what follows
-            # for the initial iteration step at time = 0.
-            if time > self.calc_tol:
-                # this is the beginning of the solver, so iteration must be set
-                # to 0.
+            #     # this is the beginning of the solver, so iteration must be set
+            #     # to 0.
             subdomain.iteration_number = 0
             for phase in subdomain.has_phases:
                 # set the solution of the previous timestep as new prev_timestep pressure
                 subdomain.pressure_prev_timestep[phase].assign(
                     subdomain.pressure[phase]
                 )
+                # All isdone flags need to be zero before we start iterating
+                self._calc_isdone_for_phase[ind].update({phase: False})
 
         ### actual iteration starts here
         # gobal stopping criterion for the iteration.
@@ -410,7 +417,7 @@ class LDDsimulation(object):
             for sd_index, subdomain in self.subdomain.items():
                 # check if the calculation on subdomain has already be marked as
                 # finished.
-                if not subdomain_calc_isdone[sd_index]:
+                if not self._calc_isdone_for_subdomain[sd_index]:
                     # set subdomain iteration to current iteration
                     subdomain.iteration_number = iteration
                     if debug:
@@ -430,6 +437,10 @@ class LDDsimulation(object):
                             subsequent_iter_error.update(
                                 {phase: error_calculated}
                                 )
+                            # 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"the subsequent error in iteration {iteration} for phase {phase} is {subsequent_iter_error[phase]}")
 
@@ -458,16 +469,18 @@ class LDDsimulation(object):
                         print("max(subsequent_iter_error.values()): ", max(subsequent_iter_error.values()))
                         total_subsequent_iter_error = max(subsequent_iter_error.values())
                         # check if error criterion has been met
-                        error_criterion = 1E-9#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
+
                         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}")
-                        if total_subsequent_iter_error < error_criterion:
+                        wetting_done = self._calc_isdone_for_phase[sd_index]['wetting']
+                        nonwetting_done = self._calc_isdone_for_phase[sd_index]['nonwetting']
+                        if nonwetting_done and wetting_done:
                             if debug:
                                 print(f"subdomain{sd_index} reachd error tol = {error_criterion} in iteration {iteration}.")
-                            subdomain_calc_isdone[sd_index] = True
+                            self._calc_isdone_for_subdomain[sd_index] = True
 
                     # prepare next iteration
                     # write the newly calculated solution to the inteface dictionaries
@@ -491,7 +504,7 @@ class LDDsimulation(object):
                     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 subdomain_calc_isdone.items():
+                    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
@@ -522,11 +535,13 @@ 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']
@@ -576,7 +591,9 @@ class LDDsimulation(object):
 
                 if debug and subdomain.mesh.num_cells() < 36:
                     print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
-
+            else:
+                print(f"calculation for phase {phase} is already done in timestep = {self.timestep_num},  t = {self.t}. skipping ..." )
+                pass
 
     def _init_solution_files(self):
         """ set up solution files for saving the solution of the LDD scheme for
@@ -588,7 +605,7 @@ class LDDsimulation(object):
         for subdom_ind, subdomain in self.subdomain.items():
             tau = subdomain.timestep_size
             L = subdomain.L
-            h = subdomain.mesh.hmin()
+            h = subdomain.mesh.hmax()
             if subdomain.isRichards:
                 Lam = subdomain.lambda_param['wetting']
                 name1 ="subdomain" + "{}".format(subdom_ind)
@@ -836,7 +853,7 @@ class LDDsimulation(object):
             # marker value gets internally set to global_index+1.
             # The reason is that we want to mark outer boundaries with the value 0.
             self.interface[num].mark(self.interface_marker, self.interface[num].marker_value)
-            print("Global interface vertices:")
+            # print("Global interface vertices:")
             self.interface[num].init_interface_values(self.interface_marker, self.interface[num].marker_value)
 
         if self.write2file['meshes_and_markers']:
@@ -887,6 +904,8 @@ 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()})
 
 
 
@@ -1059,21 +1078,15 @@ class LDDsimulation(object):
 
         subdomain = self.subdomain[subdomain_index]
         for phase in subdomain.has_phases:
-            if np.abs(time-self.starttime) < self.calc_tol:
+            if np.abs(time-self.starttime) < self.tol:
                 # here t = 0 and we have to initialise the sources.
-                mesh = subdomain.mesh
-                # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
                 f = self.sources[subdomain_index][phase]
-                # note that the default interpolation degree is 2
-                f_expression = {
-                    phase: df.Expression(f, domain = mesh, degree = interpolation_degree, t = time)
-                    }
-
                 subdomain.source.update(
-                    {phase: f_expression[phase]}#
+                    {phase: df.Expression(f, domain = subdomain.mesh,#
+                                             degree = interpolation_degree, #
+                                             t = time)}
                     )
             else:
-                # note that the default interpolation degree is 2
                 subdomain.source[phase].t = time
 
     def _add_parameters_to_output_dirname(self):