From a3a9d4027bcf375a1687fce158aa5ce0e5413d35 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Tue, 2 Apr 2019 18:00:50 +0200
Subject: [PATCH] write LDDsolver method.

---
 LDDsimulation/LDDsimulation.py | 81 ++++++++++++++++++++++++++++++++++
 LDDsimulation/__init__.py      |  3 ++
 LDDsimulation/domainPatch.py   |  5 ++-
 3 files changed, 88 insertions(+), 1 deletion(-)

diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index f137716..f6ed481 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -80,6 +80,10 @@ class LDDsimulation(object):
         self.adjacent_subdomains = None
 
         ## Private variables
+        # maximal number of L-iterations that the LDD solver uses.
+        self._max_iter_num = 100
+        # TODO rewrite this with regard to the mesh sizes
+        self.calc_tol = 10^-6
         # The number of subdomains are counted by self.init_meshes_and_markers()
         self._number_of_subdomains = 0
         # variable to check if self.set_parameters() has been called
@@ -145,6 +149,83 @@ class LDDsimulation(object):
         # init Dirichlet Boundary values for the first time
         self.set_DirichletBC(time = 0, first_init = True)
 
+    def LDDsolver(self, time: float):
+        """ calulate the solution on all patches for the next timestep
+
+        The heart of the simulation. The function LDDsolver implements the LDD
+        iteration, that is solves the problem given the initial data at timestep
+        time.
+        """
+        iteration = 0
+        # 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
+        for ind, subdomain in self.subdomain.items():
+            subdomain_calc_isdone.update({ind, False})
+            # 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.
+            if time > 0:
+                # this is the beginning of the solver, so iteration must be set
+                # to 0.
+                subdomain.iteration_number = 0
+                # set the solution of the previous timestep as new prev_timestep pressure
+                subdomain.pressure_prev_timestep = subdomain.pressure
+                # TODO recheck if
+                subdomain.pressure_prev_iter = subdomain.pressure
+                # needs to be set also
+        ### actual iteration starts here
+        # gobal stopping criterion for the iteration.
+        all_subdomains_done = False
+        while iteration < self._max_iter_num and not all_subdomains_done:
+            # we need to loop over the subdomains and solve an L-scheme type step
+            # on each subdomain. here it should be possible to parallelise in
+            # the future.
+            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]:
+                    # solve the problem on subdomain
+                    self.Lsolver_step(subdomain = subdomain,#
+                                      iteration = iteration,#
+                                      debug = True
+                    )
+                    subsequent_iterations_err = self.calc_iteration_error(
+                        subdomain = subdomain,#
+                        error_type = "L2"
+                    )
+                    # check if error criterion has been met
+                    if subsequent_iterations_err < self.calc_tol:
+                        subdomain_calc_isdone[sd_index] = True
+
+                    # prepare next iteration
+                    # write the newly calculated solution to the inteface dictionaries
+                    # for communication
+                    subdomain.write_pressure_to_interfaces()
+                    subdomain.pressure_prev_iter = subdomain.pressure
+                else:
+                    # one subdomain is done. Check if all subdomains are done to
+                    # stop the while loop if needed.
+                    # For this, set
+                    all_subdomains_done = True
+                    # then check if all domains are done. If not, reset
+                    # all_subdomains_done to False.
+                    for subdomain, isdone in subdomain_calc_isdone.items():
+                        if not isdone:
+                            all_subdomains_done = False
+                    #TODO check if maybe I still need to do
+                    # subdomain.iteration_number += 1
+                    # or adapt the calc_gli_term method to reflect that one subdomain
+                    # could be done earlier than others.
+
+
+            # you wouldn't be so stupid as to write an infinite loop, would you?
+            iteration += 1
+            # end iteration while loop.
+
     ## Private methods
     def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
                                 mesh_resolution: float = None) -> None:
diff --git a/LDDsimulation/__init__.py b/LDDsimulation/__init__.py
index e69de29..a0c87a4 100644
--- a/LDDsimulation/__init__.py
+++ b/LDDsimulation/__init__.py
@@ -0,0 +1,3 @@
+from . import domainPatch
+from . import helpers
+from . import LDDsimulation
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 47d72f1..35ddc9d 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -180,7 +180,10 @@ class DomainPatch(df.SubDomain):
 
     #### PUBLIC METHODS
     def write_pressure_to_interfaces(self):
-        """ save the interface values of self.pressure to all neighbouring interfaces"""
+        """ save the interface values of self.pressure to all neighbouring interfaces
+
+        This method is used by the LDDsimulation.LDDsolver method.
+        """
         if self.isRichards:
             for ind in self.has_interface:
                 interface[ind].read_pressure_dofs(from_function = self.pressure['wetting'], #
-- 
GitLab