From b1af8421a8a83f4f8f02d3205c71e559c1d11956 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Thu, 2 May 2019 09:22:56 +0200
Subject: [PATCH] rename _eval_sources() -> _init_sources(). Change LDDsolver
 accordingly to update source expressions.

---
 LDDsimulation/LDDsimulation.py | 58 ++++++++++++++++++----------------
 1 file changed, 31 insertions(+), 27 deletions(-)

diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 8161e99..1c2bf30 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -100,7 +100,7 @@ class LDDsimulation(object):
 
         ## Private variables
         # maximal number of L-iterations that the LDD solver uses.
-        self._max_iter_num = 1
+        self._max_iter_num = 1000
         # TODO rewrite this with regard to the mesh sizes
         # self.calc_tol = self.tol
         # list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps
@@ -161,6 +161,8 @@ class LDDsimulation(object):
             'maximum_iterations': 1000,
             'report': False
         }
+
+        self.debug_solver = 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
@@ -195,7 +197,7 @@ class LDDsimulation(object):
                        number_of_timesteps_to_analyse: int,#
                        sources: tp.Dict[int, tp.Dict[str, str]],#
                        initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
-                       dirichletBC_Expressions: tp.Dict[int, tp.Dict[str, str]],#
+                       dirichletBC_expression_strings: tp.Dict[int, tp.Dict[str, str]],#
                        write2file: tp.Dict[str, bool],#
                        exact_solution: tp.Dict[int, tp.Dict[str, str]] = None, #
                        )-> None:
@@ -223,7 +225,7 @@ class LDDsimulation(object):
         self.timestep_size = timestep_size
         self.sources = sources
         self.initial_conditions = initial_conditions
-        self.dirichletBC_expression_strings = dirichletBC_Expressions
+        self.dirichletBC_expression_strings = dirichletBC_expression_strings
         self.exact_solution = exact_solution
         self.write2file = write2file
         self._parameters_set = True
@@ -237,22 +239,25 @@ class LDDsimulation(object):
             raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)')
         self._add_parameters_to_output_dirname()
         self._init_analyse_timesteps()
-        print("initialise meshes and marker functions ...\n")
+        df.info(colored("initialise meshes and marker functions ...\n", "orange"))
         self._init_meshes_and_markers()
-        print("initialise interfaces ...\n")
+        df.info(colored("initialise interfaces ...\n", "orange"))
         self._init_interfaces()
-        print("initialise subdomains ...\n")
+        df.info(colored("initialise subdomains ...\n", "orange"))
         self._init_subdomains()
-        print("creating xdmf files for each subdomain ...\n")
+        df.info(colored("creating xdmf files for each subdomain ...\n", "orange"))
         self._init_solution_files()
         # initial values get initialised here to be able to call the solver at different
         # timesteps without having to call self.run().
+        df.info(colored("create initial values ...", "orange"))
         self._init_initial_values()
+        df.info(colored("create source Expressions ...", "orange"))
+        self._init_sources()
         # initialise exact solution expression dictionary if we have an exact
         # solution case.
         self._init_exact_solution_expression()
         if self.exact_solution:
-            print("writing exact solution for all time steps to xdmf files ...")
+            df.info(colored("writing exact solution for all time steps to xdmf files ...", "orange"))
             self.write_exact_solution_to_xdmf()
         self._initialised = True
 
@@ -289,6 +294,8 @@ class LDDsimulation(object):
         df.info(colored("Start time stepping","green"))
         # write initial values to file.
         self.timestep_num = int(0)
+        # note that at this point of the code self.t = self.starttime as set in
+        # the beginning.
         while self.timestep_num < self.number_of_timesteps:
             print(f"entering timestep ",
                   "{}".format(self.timestep_num),
@@ -296,9 +303,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 = True, analyse_timestep = True)
+                self.LDDsolver(debug = self.debug_solver, analyse_timestep = True)
             else:
-                self.LDDsolver(debug = True, analyse_timestep = False)
+                self.LDDsolver(debug = self.debug_solver, analyse_timestep = False)
 
             self.t += self.timestep_size
             self.timestep_num += 1
@@ -390,8 +397,8 @@ class LDDsimulation(object):
             subdomain.null_all_interface_iteration_numbers()
             # update the time of the source terms before starting the iteration.
             # The sources and sinks are the same throughout the iteration.
-            self._eval_sources(time = time,
-                               subdomain_index = ind)
+            # self._eval_sources(time = time,
+            #                    subdomain_index = ind)
             self.update_DirichletBC_dictionary(time = time, #
                 subdomain_index = ind,
                 )
@@ -399,6 +406,8 @@ class LDDsimulation(object):
             #     # to 0.
             subdomain.iteration_number = 0
             for phase in subdomain.has_phases:
+                # evaluate the sources to the current time.
+                subdomain.source[phase].t = time
                 # set the solution of the previous timestep as new prev_timestep pressure
                 subdomain.pressure_prev_timestep[phase].assign(
                     subdomain.pressure[phase]
@@ -475,9 +484,9 @@ class LDDsimulation(object):
                             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}")
-                        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:
+                        # wetting_done = self._calc_isdone_for_phase[sd_index]['wetting']
+                        # nonwetting_done = self._calc_isdone_for_phase[sd_index]['nonwetting']
+                        if total_subsequent_iter_error < error_criterion:
                             if debug:
                                 print(f"subdomain{sd_index} reachd error tol = {error_criterion} in iteration {iteration}.")
                             self._calc_isdone_for_subdomain[sd_index] = True
@@ -973,7 +982,7 @@ class LDDsimulation(object):
             subdomain.write_pressure_to_interfaces()
             subdomain.write_solution_to_xdmf(#
                 file = self.solution_file[num], #
-                time = self.t,#
+                time = self.starttime,#
                 write_iter_for_fixed_time = False,
                 )
 
@@ -1060,27 +1069,22 @@ class LDDsimulation(object):
         else:
             pass
 
-    def _eval_sources(self, subdomain_index: int, #
-                            interpolation_degree: int = None, #
-                            time: float = 0,#
-                            ):
-        """ evaluate time dependent source terms or initialise them if time == 0
+    def _init_sources(self, interpolation_degree: int = None):
+        """ create source expressions at self.starttime and populate subdomain
+        source  dictionaries.
         """
         if interpolation_degree is None:
             interpolation_degree = self.interpolation_degree
 
-        subdomain = self.subdomain[subdomain_index]
-        for phase in subdomain.has_phases:
-            if np.abs(time-self.starttime) < self.tol:
+        for subdomain_index, subdomain in self.subdomain.items():
+            for phase in subdomain.has_phases:
                 # here t = 0 and we have to initialise the sources.
                 f = self.sources[subdomain_index][phase]
                 subdomain.source.update(
                     {phase: df.Expression(f, domain = subdomain.mesh,#
                                              degree = interpolation_degree, #
-                                             t = time)}
+                                             t = self.starttime)}
                     )
-            else:
-                subdomain.source[phase].t = time
 
     def _add_parameters_to_output_dirname(self):
         """ function to create a new output directory for given set of input
-- 
GitLab