rewrite self._eval_sources() to operate on a single domain

LDDsimulation/LDDsimulation.py

@@ -145,9 +145,6 @@ class LDDsimulation(object):
         self._init_interfaces()
         self._init_subdomains()
         self._init_initial_values()
-        self._eval_sources()
-        # 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
@@ -165,6 +162,11 @@ class LDDsimulation(object):
         # and
         for ind, subdomain in self.subdomain.items():
             subdomain_calc_isdone.update({ind, False})
+            # update the time of the source terms before starting the iteration.
+            # The sources and sinks are the same throughout the iteration.
+            self._eval_sources(interpolation_degree = 2,
+                               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.
@@ -177,6 +179,7 @@ class LDDsimulation(object):
                 # 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
@@ -189,12 +192,12 @@ class LDDsimulation(object):
                 # finished.
                 if not subdomain_calc_isdone[sd_index]:
                     # solve the problem on subdomain
-                    self.Lsolver_step(subdomain = subdomain,#
+                    self.Lsolver_step(subdomain_index = sd_index,#
                                       iteration = iteration,#
                                       debug = True
                     )
                     subsequent_iterations_err = self.calc_iteration_error(
-                        subdomain = subdomain,#
+                        subdomain_index = sd_index,#
                         error_type = "L2"
                     )
                     # check if error criterion has been met
@@ -226,6 +229,26 @@ class LDDsimulation(object):
             iteration += 1
             # end iteration while loop.
 
+    def Lsolver_step(subdomain_index: int],#
+                     iteration: int,#
+                     debug: bool = False) -> None:
+        """ L-scheme solver iteration step for an object of class subdomain
+
+        Lsolver_step implements L-scheme solver iteration step for an subdomain
+        with index subdomain_index, i.e. for the model form (already in L-scheme
+        form) stored in the domain patch subdomain. 
+
+        # parameters
+        subdomain_index:    int                         subdomain index that defines
+                                                        which subdomain to perform
+                                                        the iteration on.
+        iteration:          int                         iteration number
+        debug:              bool                        flag to toggle weather or not
+                                                        to write out each iteration.
+        """
+
+
+
     ## Private methods
     def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
                                 mesh_resolution: float = None) -> None:
@@ -444,16 +467,18 @@ class LDDsimulation(object):
                     {'nonwetting': df.DirichletBC(V['nonwetting'], pDCnw, boundary_marker, 1)},#
                 )
 
-    def _eval_sources(self, interpolation_degree: int = 2, time: float = 0):
+    def _eval_sources(self, interpolation_degree: int = 2, #
+                            time: float = 0,
+                            subdomain_index: int):
         """ evaluate time dependent source terms or initialise them if time == 0
         """
+        subdomain = self.subdomain[subdomain_index]
         if np.abs(time) < self.tol:
             # here t = 0 and we have to initialise the sources.
-            for num, subdomain in self.subdomain.items():
             mesh = subdomain.mesh
             V = subdomain.function_space
-                # p0 contains both pw_0 and pnw_0 for subdomain[num]
-                f = self.sources[num]
+            # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
+            f = self.sources[subdomain_index]
             if subdomain.isRichards:
                 # note that the default interpolation degree is 2
                 fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
@@ -462,25 +487,14 @@ class LDDsimulation(object):
             else:
                 fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
                 fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time)
-                    # subdomain.pressure_prev_iter = {'wetting' : df.interpolate(fw, V['wetting']),#
-                    #                                 'nonwetting' : df.interpolate(fnw, V['nonwetting'])}
-                    # subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(fw, V['wetting']),#
-                    #                                     'nonwetting' : df.interpolate(fnw, V['nonwetting'])}
                 subdomain.source = {'wetting' : fw,#
                                     'nonwetting' : fnw}
         else:
-            for num, subdomain in self.subdomain.items():
-                # mesh = subdomain.mesh
-                # V = subdomain.function_space
-                # # p0 contains both pw_0 and pnw_0 for subdomain[num]
-                # f = self.sources[num]
+            # here the sources are already initiated and stored in the subdomain.
+            # they must be updated in time.
             if subdomain.isRichards:
                 # note that the default interpolation degree is 2
-                    # fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
-                    # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
                 subdomain.source['wetting'].t = time
             else:
-                    # fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
-                    # fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time)
                 subdomain.source['wetting'].t = time
                 subdomain.source['nonwetting'].t = time