put the gravity flux term function in unused code

LDDsimulation/unused_code.py

+# this was part of domainPatch.DomainPatch
+
+def calc_gravity_neumann_flux(self,
+                  interface_index: int,
+                  phase: str,
+                  initial: bool = False,
+                  debug: bool = False
+                  ) -> None:
+    """calculate the gravity term of the neumann flux
+
+    This method is written such that it can calculate the gravity neumann
+    flux for both phases in case it is needed, but is actually only used
+    for the flux of the nonwetting gravity term of a Richards domain with
+    TPR coupling.
+    """
+    subdomain = self.subdomain_index
+    dt = self.timestep_size
+    Lambda = self.lambda_param
+    interf_ind = interface_index
+
+    interface = self.interface[interf_ind]
+
+    if initial:
+        pressure_previous = self.pressure_prev_timestep
+    else:
+        pressure_previous = self.pressure_prev_iter
+
+    if self.isRichards:
+        # assuming that we normalised atmospheric pressure to zero,
+        # pc = -pw in the Richards case.
+        pc_prev_iter = -pressure_previous['wetting']
+    else:
+        pw_prev_iter = pressure_previous['wetting']
+        pnw_prev_iter = pressure_previous['nonwetting']
+        pc_prev_iter = pnw_prev - pw_prev
+
+    # n = df.FacetNormal(self.mesh)
+    S = self.saturation
+
+    # viscosity
+    mu_a = self.viscosity[phase]
+    # relative permeability
+    ka = self.relative_permeability[phase]
+    z_a = self.gravity_term[phase]
+
+    if phase == 'wetting':
+        # z_a included the minus sign, see self._calc_gravity_expressions
+        flux = 1/mu_a*ka(S(pc_prev_iter))*df.grad(z_a)
+    else:
+        # phase == 'nonwetting'
+        flux = 1/mu_a*ka(1-S(pc_prev_iter))*df.grad(z_a)
+
+    flux_function = df.project(flux, self.function_space["flux"][phase])
+    flux_x, flux_y = flux_function.split()
+    gravity_neumann_flux = df.Function(self.function_space["gli"][phase])
+    # # get dictionaries of global dof indices and coordinates along
+    # # the interface. These dictionaries have the facet indices of
+    # # facets belonging to the interface as keys (with respect to
+    # # the submesh numbering) and the dictionary
+    # # containing pairs {global_dof_index: dof_coordinate} for all
+    # # dofs along the facet as values.
+    neumann_flux_dof_coords = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
+    # p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
+    # flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
+    # # the attributes local and global for facet numbering mean
+    # # the following: local facet index is the index of the facet
+    # # with respect to the numbering of the submesh belonging to
+    # # the subdomain. global facet index refers to the facet numbering
+    # # of the mesh of the whole computational domain.
+    local2global_facet = interface.local_to_global_facet_map[subdomain]
+    corresponding_dof_index = self.interface_corresponding_dof_index[interf_ind][phase]
+    for local_facet_ind, normal in interface.outer_normals[subdomain].items():
+        neumann_flux_dofs_along_facet = neumann_flux_dof_coords[local_facet_ind]
+        global_facet_ind = local2global_facet[local_facet_ind]
+        for neumann_flux_dof_index, neumann_flux_dof_coord in neumann_flux_dofs_along_facet.items():
+            flux_x_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_x"]
+            flux_y_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_y"]
+
+            gravity_neumann_flux.vector()[neumann_flux_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+                + flux_y.vector()[flux_y_dof_index]*normal[1]
+
+            # save this newly calculated dof to the interface
+            # dictionary neumann_flux_save_dict
+            dof_coord_tupel = (neumann_flux_dof_coord[0], neumann_flux_dof_coord[1])
+            if initial:
+                # because we don't know if the Richards domain or the TP domain
+                # will be up next in the solver, we need to save the initial
+                # gravity neumann flux to both the current and the prev dictionary.
+                interface.gravity_neumann_flux_prev[subdomain][phase][global_facet_ind].update(
+                    {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
+                    )
+                interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
+                    {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
+                    )
+            else:
+                # in the noninitial case we always save to the current dictionary
+                # similarly to the pressures.
+                interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
+                    {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
+                    )
+
+
+def read_gravity_neumann_flux(self,
+                  interface_index: int,
+                  phase: str,
+                  debug: bool = False) -> df.Function: # tp.Dict[str, fl.Form]
+    """in the case of a TPR coupling read the additional gravity neumann flux
+    for the nonwetting phase from the neighbour
+
+    IMPLEMENTATION
+    same logic as in the calc_gli_term method applies.
+
+    PARAMETERS
+    ------------------------------------
+    interface_index    int       interface index for which to calculate the
+                                 gravity_neumann_flux term.
+    phase              str       phase to calculate the gravity_neumann_flux term for.
+    # DEBUG:           bool      toggle the debug output.
+    """
+    # this is the number of the current iteration of the subdomain.
+    iteration = self.iteration_number
+    subdomain = self.subdomain_index
+    interface = self.interface[interface_index]
+    # gravity_neumann_flux should be initialised by zero.
+    gravity_neumann_flux = df.Function(self.function_space["gli"][phase])
+
+    # neighbour index
+    neighbour = interface.neighbour[subdomain]
+    # update the current_iteration number of subdomain stored in the interface
+    # to the current iteration number of the subdomain.
+    interface.current_iteration[subdomain] = iteration
+    # save the iteration number of the neighbour
+    neighbour_iter_num = interface.current_iteration[neighbour]
+
+    # read previous neighbouring gravity neumann flux from interface.
+    # depending on our iteration number and the iteration number of
+    # the neighbour, weed to either read from the previous pressure
+    # dictionary or the current one.
+    if iteration == neighbour_iter_num + 1:
+        # in this case the neighbour has not yet iterated beyond
+        # the iteration we are currently in, so
+        # interface.gravity_neumann_flux[neighbour] holds the "previous neumann gravity flux"
+        # we need to read these values from the interface.
+        gravity_flux_read_dict = interface.gravity_neumann_flux[neighbour][phase]
+
+    else:
+        # iteration == neighbour_iter_num:
+        # this is the only possible other case. In this case the
+        # neighbour has iterated once beyond the iteration we are
+        # currently in, so interface.gravity_neumann_flux_prev holds the "previous neumann gravity flux"
+        # we need to read from the interface.
+        gravity_flux_read_dict = interface.gravity_neumann_flux_prev[neighbour][phase]
+
+    flux_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
+    local2global_facet = interface.local_to_global_facet_map[subdomain]
+    for local_facet_index, flux_facet_dofs2coordinates in flux_interface_dofs.items():
+        global_facet_index = local2global_facet[local_facet_index]
+        # read gravity_neumann_flux_prev term from the neighbour.
+        neighbour_flux_on_facet = gravity_flux_read_dict[global_facet_index]
+
+        for flux_dof_index, flux_dof_coord in flux_facet_dofs2coordinates.items():
+            # same with the gravity_neumann_flux_previous_dof
+            gravity_neumann_flux_neighbour_dof = None
+            for flux_neighbour_coord_tupel, flux_neighbour_dof in neighbour_flux_on_facet.items():
+                flux_neighbour_coord = np.array([flux_neighbour_coord_tupel[0], flux_neighbour_coord_tupel[1]])
+                # due to the warning in the documentation of np.allclose
+                # that np.allclose is not symetric, we test if both
+                # np.allclose(a,b) and np.allclose(b,a) are true.
+                a_close_b = np.allclose(flux_dof_coord, flux_neighbour_coord, rtol=1e-10, atol=1e-14)
+                b_close_a = np.allclose(flux_neighbour_coord, flux_dof_coord, rtol=1e-10, atol=1e-14)
+                if a_close_b and b_close_a:
+                    gravity_neumann_flux_neighbour_dof = flux_neighbour_dof
+                    break
+
+            if gravity_neumann_flux_neighbour_dof is None:
+                raise RuntimeError(f"didn't find gravity_neumann_flux_neighbour_dof corresponding to dict key ({flux_dof_coord})",
+                                    "something is wrong")
+
+            gravity_neumann_flux.vector()[flux_dof_index] = gravity_neumann_flux_neighbour_dof
+
+        # In the following case, the gravity_neumann_flux of the neighbour is saved to
+        # gravity_neumann_flux currently and will be overwritten in the next step,
+        # if we don't save it to gravity_neumann_flux_prev of the neighbour.
+        if iteration == neighbour_iter_num:
+            # print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
+            #       f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
+            interface.gravity_neumann_flux_prev[neighbour][phase].update(#
+                {global_facet_index: interface.gravity_neumann_flux[neighbour][phase][global_facet_index].copy()}
+                )
+
+    return gravity_neumann_flux
+### END read_gravity_neumann_flux