add methods to enforce gravity neumann flux

LDDsimulation/LDDsimulation.py

@@ -535,6 +535,13 @@ class LDDsimulation(object):
                     subdomain.pressure_prev_iter[phase].assign(
                         subdomain.pressure[phase]
                     )
+                # onyl becomes active in TPR case
+                if self.include_gravity and subdomain.isRichards and subdomain.TPR_occures:
+                        for interface in subdomain.has_interface_with_TP_neighbour:
+                            subdomain.calc_gravity_neumann_flux(
+                                interface_index=interface,
+                                phase="nonwetting"
+                                )
             # end loop over subdomains
             ##### stopping criterion for the solver.
             total_subsequent_error = np.amax(max_subsequent_error)
@@ -645,6 +652,17 @@ class LDDsimulation(object):
                 subdomain.write_pressure_to_interfaces()
                 subdomain.write_pressure_to_interfaces(previous_iter=True)
                 subdomain.calc_gl0_term()
+                # If gravity is included and a Richards subdomain has a TPR coupling,
+                # we need to calculate the Neumann flux of the nonwetting phase
+                # that appears due to gravity
+                if self.include_gravity:
+                    if subdomain.isRichards and subdomain.TPR_occures:
+                        for interface in subdomain.has_interface_with_TP_neighbour:
+                            subdomain.calc_gravity_neumann_flux(
+                                initial=True,
+                                interface_index=interface,
+                                phase="nonwetting"
+                                )
 
         return solution_over_iteration_within_timestep
 

LDDsimulation/boundary_and_interface.py

@@ -266,6 +266,8 @@ class Interface(BoundaryPart):
         # dictionaries holding the forms for the decoupling terms
         self.gli_term = dict()
         self.gli_term_prev = dict()
+        self.gravity_neumann_flux = dict()
+        self.gravity_neumann_flux_prev = dict()
         # dictionary holding the current iteration number i for each of the adjacent
         # subdomains with index in adjacent_subdomains
         self.current_iteration = {#
@@ -613,6 +615,8 @@ class Interface(BoundaryPart):
             else:  #function == "gli":
                 self.gli_term.update({subdomain_index: dict()})
                 self.gli_term_prev.update({subdomain_index: dict()})
+                self.gravity_neumann_flux.update({subdomain_index: dict()})
+                self.gravity_neumann_flux_prev.update({subdomain_index: dict()})
 
             for phase in has_phases:
                 if function == "pressure":
@@ -621,7 +625,8 @@ class Interface(BoundaryPart):
                 else:  #function == "gli":
                     self.gli_term[subdomain_index].update({phase: dict()})
                     self.gli_term_prev[subdomain_index].update({phase: dict()})
-
+                    self.gravity_neumann_flux[subdomain_index].update({phase: dict()})
+                    self.gravity_neumann_flux_prev[subdomain_index].update({phase: dict()})
                 space = function_space[function][phase]
                 dof_coordinates = self.dofs_and_coordinates(space, interface_marker, interface_marker_value)
                 # the dof dictionaries are grouped by the global (with respect to the
@@ -635,6 +640,8 @@ class Interface(BoundaryPart):
                     else:  #function == "gli":
                         self.gli_term[subdomain_index][phase].update({global_mesh_facet_index: dict()})
                         self.gli_term_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+                        self.gravity_neumann_flux[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+                        self.gravity_neumann_flux_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
 
                     for dof_coord in facet_dof_coord_dict.values():
                         # we want to use the coordinates as key in a dictionary.
@@ -648,6 +655,8 @@ class Interface(BoundaryPart):
                         else:  #function == "gli":
                             self.gli_term[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                             self.gli_term_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+                            self.gravity_neumann_flux[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+                            self.gravity_neumann_flux_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                     # end init_interface_values