From 715daf79bcd2f8eceddbe12bc4f90edf37634da1 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Mon, 1 Apr 2019 16:17:43 +0200
Subject: [PATCH] finish domainPatch.calc_gli_term method

---
 domainPatch.py | 193 ++++++++++++++++++++++++++++++++++---------------
 1 file changed, 135 insertions(+), 58 deletions(-)

diff --git a/domainPatch.py b/domainPatch.py
index b67873c..1f69ce3 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -209,6 +209,7 @@ class DomainPatch(df.SubDomain):
         dx = self.dx
         ds = self.ds
         dt = self.timestep_size
+        porosity = self.porosity
         if self.isRichards:
             if phase == 'nonwetting':
                 print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
@@ -241,7 +242,7 @@ class DomainPatch(df.SubDomain):
             #     gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
             #     rhs_interface_forms.append((dt*gli*phi_w)*ds[index])
             rhs1 = (Lw*pw_prev_iter*phi_w)*dx
-            rhs2 = -((S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
+            rhs2 = -(porosity*(S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
             rhs_without_gli = rhs1 + rhs2 + dt*source_w*phi_w*dx
             form_and_rhs = {#
                 'form': form,#
@@ -278,7 +279,7 @@ class DomainPatch(df.SubDomain):
                 form = form1 + form2 + df.sum(interface_forms)
 
                 rhs1 = (La*pw_prev_iter*phi_a)*dx
-                rhs2 = -((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
+                rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
                 rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
 
                 form_and_rhs = {#
@@ -305,7 +306,7 @@ class DomainPatch(df.SubDomain):
 
                 rhs1 = (La*pnw_prev_iter*phi_a)*dx
                 # the non wetting phase has a + before instead of a - before the bracket
-                rhs2 = ((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
+                rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
                 rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
 
                 form_and_rhs = {#
@@ -327,6 +328,9 @@ class DomainPatch(df.SubDomain):
         ds form to be added to the rhs by the solver. This is due to the need of
         having to communicate dofs over the interface. Additionally, the dofs are
         being saved to an interface dictionary exactly as the pressures.
+
+        Caution: The array self.gli_assembled needs to be added to the rhs with
+        a minus sign by the solver!
         """
         for ind in self.has_interface:
             # self._calc_gli_term_on(interface, iteration)
@@ -351,13 +355,13 @@ class DomainPatch(df.SubDomain):
                 if not neighbour_iter_num == 0:
                     # save the gl term from the neighbour first to use it in the
                     # next iteration
-                    if self.isRichards:
+                    interface.gli_term_prev[subdomain].update(#
+                        {'wetting': interface.gli_term[neighbour]['wetting']}
+                        )
+                    # in case we have two-phase and the neighbour is two-phase, two,
+                    # we need to Additionally save the nonwetting values.
+                    if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
                         interface.gli_term_prev[subdomain].update(#
-                            {'wetting': interface.gli_term[neighbour]['wetting']}
-                            )
-                    else:
-                        interface.gli_term_prev[subdomain].update(#
-                            {'wetting': interface.gli_term[neighbour]['wetting']},#
                             {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
                             )
 
@@ -436,58 +440,131 @@ class DomainPatch(df.SubDomain):
                 interface.current_iteration[subdomain] += 1
 
             else: # iteration > 0
+                # the wetting phase update is always to be calculated even if
+                # only Richards is assumed on the current patch or the neighbour
+                # only assumes the Richards equation.
+
+                # in case neighbour_iter_num == iteration we need to first save
+                # the gli_term of the neighbour to gli_term_prev.
                 if neighbour_iter_num == iteration:
-                    if self.isRichards:
-                        # first, save the gli_terms of the neighbour to gli_prev_term
-                        # of the current subdomain.
+                    # in this case, save the gli_terms of the neighbour to gli_prev_term
+                    # of the current subdomain first. this will then be used by
+                    # the subsequent code to calculate the gli term with
+                    interface.gli_term_prev[subdomain].update(#
+                        {'wetting': interface.gli_term[neighbour]['wetting']}
+                        )
+
+                    if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
                         interface.gli_term_prev[subdomain].update(#
-                            {'wetting': interface.gli_term[neighbour]['wetting']}
+                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
                             )
-                        # then read gli_prev term from interface and write it to
-                        # the subdomain gli_assembled array
-                        interface.write_gli_dofs(to_array = self.gli_assembled['wetting'], #
-                                       interface_dofs = interface_dofs['wetting'],#
-                                       dof_to_vert_map = self.dof2vertex['wetting'],#
-                                       local_to_parent_vertex_map = self.parent_mesh_index,#
-                                       phase = 'wetting',#
-                                       subdomain_ind = subdomain,#
-                                       previous_iter = True
-                                       )
-
-                        # read the neighbouring pressure values and write them to
-                        # the placeholder fucntion self.neighbouring_interface_pressure
-                        interface.write_pressure_dofs(to_function = self.neighbouring_interface_pressure['wetting'], #
-                                       interface_dofs = interface_dofs['wetting'],#
-                                       dof_to_vert_map = self.dof2vertex['wetting'],#
-                                       local_to_parent_vertex_map = self.parent_mesh_index,#
-                                       phase = 'wetting',#
-                                       subdomain_ind = subdomain,#
-                                       previous_iter = True)
-
-
-                        pw_prev = pw_prev_iter['wetting']
-                        #calc gl0 from the flux term
-                        flux = -1/mu_w*kw(S(pw_prev))*df.grad(pw_prev)
-                        gl0 = (df.dot(flux, n) - Lambda['wetting']*pw_prev)
-                        interface.gli_term[subdomain].update({'wetting': gl0})
-                    else:
-                        mu_w = mu['wetting']
-                        mu_nw = mu['nonwetting']
-                        pw_prev = pw_prev_iter['wetting']
-                        pnw_prev = pw_prev_iter['nonwetting']
-                        kw = k['wetting']
-                        knw = k['nonwetting']
-                        #calc gl0 from the flux term
-                        fluxw = -1/mu_w*kw(S(pnw_prev - pw_prev))*df.grad(pw_prev)
-                        fluxnw = -1/mu_nw*knw(1-S(pnw_prev - pw_prev))*df.grad(pnw_prev)
-                        gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
-                        gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
-                        interface.gli_term[subdomain].update({'wetting': gwl0})
-                        interface.gli_term[subdomain].update({'nonwetting': gnwl0})
 
+                # then read gli_prev term from interface and write it to
+                # the gli_assembled array of the subdomain
+                interface.write_gli_dofs(to_array = self.gli_assembled['wetting'], #
+                               interface_dofs = interface_dofs['wetting'],#
+                               dof_to_vert_map = self.dof2vertex['wetting'],#
+                               local_to_parent_vertex_map = self.parent_mesh_index,#
+                               phase = 'wetting',#
+                               subdomain_ind = subdomain,#
+                               previous_iter = True
+                               )
+
+                # read neighbouring pressure values from interface and write them to
+                # the function self.neighbouring_interface_pressure
+                interface.write_pressure_dofs(to_function = self.neighbouring_interface_pressure['wetting'], #
+                               interface_dofs = interface_dofs['wetting'],#
+                               dof_to_vert_map = self.dof2vertex['wetting'],#
+                               local_to_parent_vertex_map = self.parent_mesh_index,#
+                               phase = 'wetting',#
+                               subdomain_ind = neighbour,#
+                               previous_iter = False)
+
+                # use the above to calculate the gli update with
+                # in the paper p^{i-1}_{3-l}
+                pw_prev = self.neighbouring_interface_pressure['wetting']
+                phi_w  = self.testfunction['wetting']
+                gli_pressure_part = df.assemble(-2*dt*Lambda['wetting']*pw_prev*phi_w*ds)
+                gli_p_part = gli_pressure_part.get_local()
+                gli_prev_neighbour = self.gli_assembled['wetting']
+                self.gli_assembled.update(
+                    {'wetting': gli_p_part + gli_prev_neighbour}
+                )
+                # write the newly calculated gli term to the interface dictionary
+                interface.read_gli_dofs(from_array = self.gli_assembled['wetting'], #
+                               interface_dofs = interface_dofs['wetting'],#
+                               dof_to_vert_map = self.dof2vertex['wetting'],#
+                               local_to_parent_vertex_map = self.parent_mesh_index,#
+                               phase = 'wetting',#
+                               subdomain_ind = subdomain,#
+                               previous_iter = False
+                               )
+                # so far, only the wetting phase has been treated. Treat the nonwetting
+                # phase now. If the patch we are currently on assumes two-phase
+                # (self.isRichards == False) then we only need to calculate a gli
+                # updaet for the nonwetting phase if the neighbour also assumes
+                # two-phase equation.
+                if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
+                    # here we need to update the nonwetting phase gli
+
+                    # then read gli_prev term from interface and write it to
+                    # the gli_assembled array of the subdomain
+                    interface.write_gli_dofs(to_array = self.gli_assembled['nonwetting'], #
+                                   interface_dofs = interface_dofs['nonwetting'],#
+                                   dof_to_vert_map = self.dof2vertex['nonwetting'],#
+                                   local_to_parent_vertex_map = self.parent_mesh_index,#
+                                   phase = 'nonwetting',#
+                                   subdomain_ind = subdomain,#
+                                   previous_iter = True
+                                   )
+
+                    # read neighbouring pressure values from interface and write them to
+                    # the function self.neighbouring_interface_pressure
+                    interface.write_pressure_dofs(to_function = self.neighbouring_interface_pressure['nonwetting'], #
+                                   interface_dofs = interface_dofs['nonwetting'],#
+                                   dof_to_vert_map = self.dof2vertex['nonwetting'],#
+                                   local_to_parent_vertex_map = self.parent_mesh_index,#
+                                   phase = 'nonwetting',#
+                                   subdomain_ind = neighbour,#
+                                   previous_iter = False)
+
+                    # in the paper p^{i-1}_{3-l}
+                    pnw_prev = self.neighbouring_interface_pressure['nonwetting']
+                    phi_nw  = self.testfunction['nonwetting']
+                    gli_pressure_part = df.assemble(-2*dt*Lambda['nonwetting']*pnw_prev*phi_nw*ds)
+                    gli_p_part = gli_pressure_part.get_local()
+                    gli_prev_neighbour = self.gli_assembled['nonwetting']
+                    self.gli_assembled.update(
+                        {'nonwetting': gli_p_part + gli_prev_neighbour}
+                    )
+                    # write the newly calculated gli term to the interface dictionary
+                    interface.read_gli_dofs(from_array = self.gli_assembled['nonwetting'], #
+                                   interface_dofs = interface_dofs['nonwetting'],#
+                                   dof_to_vert_map = self.dof2vertex['nonwetting'],#
+                                   local_to_parent_vertex_map = self.parent_mesh_index,#
+                                   phase = 'nonwetting',#
+                                   subdomain_ind = subdomain,#
+                                   previous_iter = False
+                                   )
+
+
+                if neighbour_iter_num -1 == iteration:
+                    # gli_term_prev has been used by the above and can now safely
+                    # be overwritten with the gli_term of the neighbour.
+                    interface.gli_term_prev[subdomain].update(#
+                        {'wetting': interface.gli_term[neighbour]['wetting']}
+                        )
+
+                    if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
+                        interface.gli_term_prev[subdomain].update(#
+                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+                            )
                 else:
-                    print('uiae')
+                    print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
+                    " and iteration = {iteration} on subdomain {subdomain}\n",
+                    "Didn't update any gli terms.")
 
+                # finally we are done and can step up the current iteration count.
                 interface.current_iteration[subdomain] += 1
 
 
@@ -999,10 +1076,10 @@ class Interface(BoundaryPart):
                    subdomain_ind: int,#
                    previous_iter: bool = False,#
                    ) -> None:
-        """ write dofs of from_array with indices interface_dofs to self.gli_term
+        """ Read dofs of from_array with indices interface_dofs and write to self.gli_term
 
-        Write the degrees of freedom of the array from_array with indices
-        interface_dofs, i.e. the dofs of the interface, to one of the dictionaries
+        Read the degrees of freedom of the array from_array with indices
+        interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
         of the interface
             self.gli_term[subdomain_ind][phase], (if gl == True)
             self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
-- 
GitLab