diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 2b432b6b461150840bcdbdd12d0515d1c0d66ed8..c746a95b50ed46b9eeea4f698f825d9f62c88907 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -410,6 +410,8 @@ class LDDsimulation(object):
# All isdone flags need to be zero before we start iterating
self._calc_isdone_for_phase[ind].update({phase: False})
+ subdomain.calc_gl0_term()
+
### actual iteration starts here
# gobal stopping criterion for the iteration.
iteration = 0
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 9fff6770af7056a936500525ff3c9b656e79cd6f..42c97c032dff7c27372c07a0539bf62fe0767617 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -1,4 +1,4 @@
-"""Class definitions for subdomains representing soil patches for LDD-TPR simulation
+[interf_ind]"""Class definitions for subdomains representing soil patches for LDD-TPR simulation
This file contains class definitions for the class DomainPatch which defines the
subdomains used for the LDD simulation. Each subdomain of this class holds
@@ -333,8 +333,148 @@ class DomainPatch(df.SubDomain):
return form_and_rhs
+ def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
+ """calculate the gl0 terms for all interfaces of the subdomain.
+
+ assemble the gl0 terms for all interfaces of the subdomain and save them
+ to the dictionries of the interfaces. This method gets called by the LDDsolver
+ before entering the L-iterations in each time step.
+ This method assumes that self.pressure_prev_timestep has assigned the values
+ of self.pressure for all phases.
+ """
+ subdomain = self.subdomain_index
+ dt = self.timestep_size
+ Lambda = self.lambda_param
+ # since we loop over all interfaces in self.has_interface we need a temporary
+ # vectors (numpy arrays) for each phase that we will use to add all dofs of the gli terms
+ # into one array.
+ gli_assembled_tmp = dict()
+ for phase in self.has_phases:
+ # print("number of dofs: ", self.pressure_prev_iter['wetting'].vector()[:].size)
+ gli_assembled_tmp.update(
+ {phase: np.zeros(self.pressure[phase].vector()[:].size, dtype=float)}
+ )
+
+ for interf_ind in self.has_interface:
+ interface = self.interface[interf_ind]
+ # neighbour index
+ neighbour = interface.neighbour[subdomain]
+ neighbour_iter_num = interface.current_iteration[neighbour]
+ ds = self.ds(interface.marker_value)
+ # needed for the read_gli_dofs() functions
+ interface_dofs = self._dof_indices_of_interface[interf_ind]
+ # for each interface, the following is a list of dofs indices belonging to another
+ # interface. The dofs corresponding to these indices must be multiplied
+ # by 1/2 to to get an average of the gli terms for dofs belonging to
+ # two different interfaces.
+ dofs_in_common_with_other_interfaces = self._interface_has_common_dof_indices[interf_ind]
+ if debug:
+ print(f"On subdomain{subdomain}, we have:\n",
+ f"Interface{interf_ind}, Neighbour index = {neighbour}\n",
+ f"dofs on interface:\n{interface_dofs['wetting']}\n",
+ f"dof2global_vertex_map:\n{self.parent_mesh_index[self.dof2vertex['wetting'][interface_dofs['wetting']]]}\n"
+ f"and dofs common with other interfaces:\n{dofs_in_common_with_other_interfaces['wetting']}")
+
+ n = df.FacetNormal(self.mesh)
+ pa_prev_timestep = self.pressure_prev_timestep
+ k = self.relative_permeability
+ S = self.saturation
+ mu = self.viscosity
+
+ if self.isRichards:
+ # assuming that we normalised atmospheric pressure to zero,
+ # pc = -pw in the Richards case.
+ pc_prev = -pa_prev_timestep['wetting']
+ else:
+ pw_prev = pa_prev_timestep['wetting']
+ pnw_prev = pa_prev_timestep['nonwetting']
+ pc_prev = pnw_prev - pw_prev
+
+ for phase in self.has_phases:
+ # viscosity
+ mu_a = mu[phase]
+ # relative permeability
+ ka = k[phase]
+ # previous timestep pressure
+ pa_prev = pa_prev_timestep[phase]
+ # testfunction
+ phi_a = self.testfunction[phase]
+ if phase == 'nonwetting' and interface.neighbour_isRichards[subdomain]:
+ # if the neighbour of our subdomain (with index self.subdomain_index)
+ # assumes a Richards model, we don't have gli terms for the
+ # nonwetting phase and assemble the terms as zero form. We don't
+ # need to do anything, since gli_assembled_tmp[phase] has been
+ # initialised to zero.
+ pass
+ # gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
+ else:
+ if phase == 'wetting':
+ # the wetting phase is always present and we always need
+ # to calculate a gl0 term.
+ # TODO: add intrinsic permeabilty!!!!
+ # TODO: add gravitiy term!!!!!
+ flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
+ else:
+ # phase == 'nonwetting'
+ # we need anothre flux in this case
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+
+ gl0 = df.dot(flux, n) - Lambda[phase]*pa_prev
+ gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
+
+ # if dofs_in_common_with_other_interfaces[phase].size == 0
+ # there are no dofs that lie on another interface and we can
+ # skip the following step.
+ if dofs_in_common_with_other_interfaces[phase].size != 0:
+ if debug:
+ print(f"the following dofs (phase = {phase}) of interface{interf_ind} are in common with \n",
+ f" other interfaces: {dofs_in_common_with_other_interfaces[phase]}")
+ for common_dof in dofs_in_common_with_other_interfaces[phase]:
+ # from the standpoint of the subdomain we are currently on,
+ # each dof can belong maximally to two interfaces. On these
+ # vertices, common to two interfaces, the dofs of the
+ # corresponding gli terms of both interfaces are averaged
+ # to give one value for this vertex. Because we achieve this
+ # by adding the assembled gli terms together we have to
+ # multiply the dofs that appear in more than one interface
+ # by 1/2.
+ gl0_assembled[common_dof] = 0.5*gl0_assembled[common_dof]
+ # now, add the just assembled and corrected term to the Global
+ # gli_assembled_tmp vector
+ gli_assembled_tmp[phase] += gl0_assembled
+
+ # writing out glo to the current iteration dictionary of the
+ # interface takes place after the loop, see below.
+
+ ### END for ind in self.has_interface:
+ # save the result of the calculation to the subdomain and to the interface
+ # dictionaries.
+ for phase in self.has_phases:
+ # self.gli_assembled.update(
+ # {phase: gli_assembled_tmp[phase].copy()}
+ # )
+ for interf_ind in self.has_interface:
+ # self._calc_gli_term_on(interface, iteration)
+ interface = self.interface[interf_ind]
+ if debug:
+ print(f"Interface{interf_ind}.gli_term_prev[{subdomain}][{phase}]=\n",interface.gli_term_prev[subdomain][phase])
+ print(f"before writing to interface[{interf_ind}] the values gli_assembled[{phase}]=\n",gli_assembled_tmp[phase])
+ # write gli dofs to interface dicts for communiction.
+ interface.read_gli_dofs(from_array = gli_assembled_tmp[phase],#self.gli_assembled[phase], #
+ interface_dofs = interface_dofs[phase],#
+ dof_to_vert_map = self.dof2vertex[phase],#
+ local_to_parent_vertex_map = self.parent_mesh_index,#
+ phase = phase,#
+ subdomain_ind = subdomain,#
+ previous_iter = True
+ )
+ if debug:
+ print(f"gl0_terms after writing to interface[{interf_ind}] Interface[{interf_ind}].gli_term_prev[{subdomain}][{phase}]=\n",interface.gli_term_prev[subdomain][phase])
+ ### END calc_gl0_term
+
+
def calc_gli_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
- """calculate the gl terms for the iteration number of the subdomain
+ """calculate the gli terms for the iteration number of the subdomain
calculate the gl terms for the iteration number of the subdomain and save them
to self.gli_assembled. The term gets saved as an already assembled sparse