diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 8d66b564807a7dd0af207a6c581986b08f8c00ae..1aad5c959b2df321a2b9bec2a65e8009d8149de0 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -349,9 +349,20 @@ class DomainPatch(df.SubDomain):
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
+ # vector (numpy array) that we will use to add all dofs of the gli terms
+ # into one array.
+ gli_assembled_tmp = dict()
+ if self.isRichards:
+ gli_assembled_tmp = {
+ 'wetting': np.zeros(subdomain.pressure['wetting'].vector().size, dtype=float)
+ }
+ else:
+ gli_assembled_tmp = {
+ 'wetting': np.zeros(subdomain.pressure['wetting'].vector().size, dtype=float),
+ 'nonwetting': np.zeros(subdomain.pressure['nonwetting'].vector().size, dtype=float)
+ }
for ind in self.has_interface:
- # self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
# neighbour index
neighbour = interface.neighbour[subdomain]
@@ -359,6 +370,11 @@ class DomainPatch(df.SubDomain):
ds = self.ds(ind)
# needed for the read_gli_dofs() functions
interface_dofs = self._dof_indices_of_interface[ind],#
+ # for each interface, this 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[ind]
if iteration == 0:
n = df.FacetNormal(self.mesh)
@@ -381,77 +397,80 @@ class DomainPatch(df.SubDomain):
)
if self.isRichards:
- mu_w = mu['wetting']
- kw = k['wetting']
- pw_prev = pw_prev_iter['wetting']
- phi_w = self.testfunction['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)
- self.gli_assembled = {
- 'wetting': df.assemble(dt*gl0*phi_w*ds).get_local()
- }
- # write glo to the current iteration dictionary of the interface
- 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
- )
+ subdomain_has_phases = ['wetting']
+ # assuming that we normalised atmospheric pressure to zero,
+ # pc = pw in the Richards case.
+ pc_prev = pw_prev_iter['wetting']
else:
- phi_w = self.testfunction['wetting']
- mu_w = mu['wetting']
+ subdomain_has_phases = ['wetting', 'nonwetting']
pw_prev = pw_prev_iter['wetting']
pnw_prev = pw_prev_iter['nonwetting']
pc_prev = pnw_prev - pw_prev
- kw = k['wetting']
- #calc gl0 from the flux term
- if 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.
- fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
- gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
- self.gli_assembled = {
- 'wetting': df.assemble(dt*gwl0*phi_w*ds).get_local(),
- # assemble the nonwetting phase gl0 term as zero functional
- 'nonwetting': df.assemble(df.Constant(0)*ds).get_local()
- }
+
+ for phase in subdomain_has_phases:
+ # viscosity
+ mu_a = mu['phase']
+ # relative permeability
+ ka = k['phase']
+ # previous pressure iteration
+ pa_prev = pw_prev_iter['phase']
+ # testfunction
+ phi_a = self.testfunction['phase']
+ 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)
+ 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:
+ 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['alpha'] += gl0_assembled
else:
- # here we need to assemble a gli term for the nonwetting phase
- phi_nw = self.testfunction['nonwetting']
- mu_nw = mu['nonwetting']
- knw = k['nonwetting']
-
- fluxw = -1/mu_w*kw(S(pc_prev))*df.grad(pw_prev)
- fluxnw = -1/mu_nw*knw(1-S(pc_prev))*df.grad(pnw_prev)
- gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
- gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
- self.gli_assembled = {
- 'wetting': df.assemble(dt*gwl0*phi_w*ds).get_local(),
- 'nonwetting': df.assemble(dt*gnwl0*phi_nw*ds).get_local()
- }
- # write glo to the current iteration dictionary of the interface
- # wetting
- 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
- )
- # nonwetting
- 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
- )
+ # phase == 'nonwetting'
+ # if we are in the two-phase case but our neighbour assumes
+ # Richards equation, we hase a zero communication term.
+ if 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.
+ gli_assembled_tmp['phase'] += df.assemble(df.Constant(0)*ds).get_local()
+ else:
+ # in this case the neighbour assumes two-phase flow
+ # and we need to calculte a gl0 torm for the nonwetting
+ # phase.
+ # 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:
+ for common_dof in dofs_in_common_with_other_interfaces['phase']:
+ # see previous case.
+ 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 loop over phases
interface.current_iteration[subdomain] += 1
else: # iteration > 0
@@ -462,9 +481,10 @@ class DomainPatch(df.SubDomain):
# 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:
- # 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
+ # in this case, first save the gli_terms of the neighbour to
+ # gli_prev_term of the current subdomain. This is because
+ # we assume, that the neighbouring patch has done a previous
+ # iteration and we need the value to calculate the gli term with
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting']}
)
@@ -473,95 +493,56 @@ class DomainPatch(df.SubDomain):
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
-
- # 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 self.isRichards:
+ subdomain_has_phases = ['wetting']
+ else:
+ subdomain_has_phases = ['wetting', 'nonwetting']
+
+ for phase in subdomain_has_phases:
+ # in case phase == 'nonwetting' only proceed if the neighbouring
+ # subdomain does not assume a Richards model.
+ if (phase == 'wetting') or (not interface.neighbour_isRichards[subdomain]):
+ # then read gli_prev term from interface and write it to
+ # the gli_assembled array of the subdomain
+ # to be able to sum the gli terms of all interfaces together,
+ # we need a dummy vector to hold the dofs read from the interface.
+ gli_readout_tmp = np.zeros(gli_assembled_tmp[phase].size, dtype = float)
+ interface.write_gli_dofs(to_array = gli_readout_tmp, #
+ 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,#
+ # this is set to True because we need gli_prev
+ 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[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 = neighbour,#
+ previous_iter = False)
+
+ # use the above to calculate the gli update with
+ # in the paper p^{i-1}_{3-l}
+ pa_prev = self.neighbouring_interface_pressure[phase]
+ phi_a = self.testfunction[phase]
+ gli_pressure_part = df.assemble(-2*dt*Lambda[phase]*pa_prev*phi_a*ds)
+ gli_p_part = gli_pressure_part.get_local()
+ gli_prev_neighbour = gli_readout_tmp
+ gli = gli_p_part + gli_prev_neighbour
+ # check if there are shared dofs with another interface.
+ if dofs_in_common_with_other_interfaces['phase'].size != 0:
+ for common_dof in dofs_in_common_with_other_interfaces['phase']:
+ # see previous case.
+ gli[common_dof] = 0.5*gli[common_dof]
+ # now, add the just assembled and corrected term to the Global
+ # gli_assembled_tmp vector
+ gli_assembled_tmp['phase'] += gli
if neighbour_iter_num -1 == iteration:
# gli_term_prev has been used by the above and can now safely
@@ -576,12 +557,54 @@ class DomainPatch(df.SubDomain):
)
else:
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.")
+ f"and iteration = {iteration} on subdomain {subdomain}\n",
+ "It is expected that this case should only be happening if the \n",
+ f"neighbour (index {neighbour}) has already finished its iteration, i.e. reached\n",
+ f"the error criterion. In this case we should have {iteration} > {neighbour_iter_num}.\n",
+ "If this is not the case revisite DomainPatch.calc_gli_term().")
# finally we are done and can step up the current iteration count.
interface.current_iteration[subdomain] += 1
+ ### END for ind in self.has_interface:
+ # save the result of the calculation to the subdomain and to the interface
+ # dictionaries.
+ if self.isRichards:
+ self.gli_assembled.update({
+ 'wetting': gli_assembled_tmp['wetting']
+ })
+ for ind in self.has_interface:
+ # self._calc_gli_term_on(interface, iteration)
+ interface = self.interface[ind]
+ # write gli dofs to interface dicts for communiction
+ 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
+ )
+ else:
+ self.gli_assembled.update({
+ 'wetting': gli_assembled_tmp['wetting'],
+ 'nonwetting': gli_assembled_tmp['nonwetting']
+ })
+ for ind in self.has_interface:
+ # self._calc_gli_term_on(interface, iteration)
+ interface = self.interface[ind]
+ for phase in ['wetting', 'nonwetting']:
+ # write gli dofs to interface dicts for communiction
+ interface.read_gli_dofs(from_array = 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 = False
+ )
+ ### END calc_gli_term
+
def null_all_interface_iteration_numbers(self) -> None:
""" reset interface.current_iteration[subdomain] = 0 for all interfaces
of the subdomain.