From 0cd3814164f5f4a79316903ac879edd0194e5882 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Sat, 26 Oct 2019 16:51:13 +0200
Subject: [PATCH] put the gravity flux term function in unused code
---
LDDsimulation/unused_code.py | 192 +++++++++++++++++++++++++++++++++++
1 file changed, 192 insertions(+)
create mode 100644 LDDsimulation/unused_code.py
diff --git a/LDDsimulation/unused_code.py b/LDDsimulation/unused_code.py
new file mode 100644
index 0000000..0e0b007
--- /dev/null
+++ b/LDDsimulation/unused_code.py
@@ -0,0 +1,192 @@
+# 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
--
GitLab