From 2dc9b5a4583c99c7132c5e9fa2d0f70c48bac78e Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Fri, 3 May 2019 20:15:27 +0200
Subject: [PATCH] rewrite calc_gli_term()
---
LDDsimulation/LDDsimulation.py | 5 +-
LDDsimulation/domainPatch.py | 358 ++++++++++++---------------------
2 files changed, 129 insertions(+), 234 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index c746a95..764ac95 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -969,8 +969,11 @@ class LDDsimulation(object):
subdomain.pressure_prev_iter[phase].assign(pa0_interpolated)
subdomain.pressure_prev_timestep[phase].assign(pa0_interpolated)
- # populate interface dictionaries with pressure values.
+ # populate interface dictionaries with pressure values. The calc_gli_method
+ # of each subdomain reads from the interface.pressure_values_prev
+ # and from interface.pressure_values dictionary so both need to be populated.
subdomain.write_pressure_to_interfaces()
+ subdomain.write_pressure_to_interfaces(previous_iter = True)
subdomain.write_solution_to_xdmf(#
file = self.solution_file[num], #
time = self.starttime,#
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 42c97c0..d2e1bb6 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -215,24 +215,17 @@ class DomainPatch(df.SubDomain):
# END constructor
#### PUBLIC METHODS
- def write_pressure_to_interfaces(self):
+ def write_pressure_to_interfaces(self, previous_iter: int = False):
""" save the interface values of self.pressure to all neighbouring interfaces
This method is used by the LDDsimulation.LDDsolver method.
# parameters
- initial_values bool Flag to toggle wether to write self.pressure
- or self.pressure_prev_iter to the interface.
- This is needed for the first ever iteration
- at time = 0.
+ previous_iter bool Flag to toggle wether to write self.pressure
+ to interface.pressure_values or to
+ interface.pressure_values_prev
"""
- # if initial_values:
- # # this assumes, that self.pressure_prev_iter has been set by
- # # LDDsimulation._init_initial_values().
- # from_func = self.pressure_prev_iter
- # else:
from_func = self.pressure
-
for ind in self.has_interface:
for phase in self.has_phases:
self.interface[ind].read_pressure_dofs(from_function = from_func[phase], #
@@ -240,8 +233,8 @@ class DomainPatch(df.SubDomain):
dof_to_vert_map = self.dof2vertex[phase],#
local_to_parent_vertex_map = self.parent_mesh_index,#
phase = phase,#
- subdomain_ind = self.subdomain_index)
-
+ subdomain_ind = self.subdomain_index,
+ previous_iter = previous_iter)
def governing_problem(self, phase: str) -> tp.Dict[str, fl.Form]:
""" return the governing form and right hand side for phase phase as a dictionary.
@@ -518,218 +511,85 @@ class DomainPatch(df.SubDomain):
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']}")
- if iteration == 0:
- n = df.FacetNormal(self.mesh)
- p_prev_iter = self.pressure_prev_timestep
- k = self.relative_permeability
- S = self.saturation
- mu = self.viscosity
-
- if not neighbour_iter_num == 0:
- # copy the gl term from the neighbour first to use it in the
- # next iteration
- if debug:
- print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
- interface.gli_term_prev[subdomain].update(#
- {'wetting': interface.gli_term[neighbour]['wetting'].copy()}
- )
- # for glob_dof_ind, dof in interface.gli_term[neighbour]['wetting'].items():
- # interface.gli_term_prev[subdomain]['wetting'][glob_dof_ind] = dof
+ for phase in self.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]):
+ # read gli_prev term from the neighbour 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_prev_of_neighbour = np.zeros(gli_assembled_tmp[phase].size, dtype = float)
+ interface.write_gli_dofs(to_array = gli_prev_of_neighbour, #
+ 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,#
+ # this is set to True because we need gli_prev
+ previous_iter = True
+ )
if debug:
- print(f"Interface{ind}.gli_term[{neighbour}]['wetting']=\n",interface.gli_term[neighbour]['wetting'])
- print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
- # in case we have two-phase and the neighbour is two-phase, two,
- # we need to Additionally save the nonwetting values.
- if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
- interface.gli_term_prev[subdomain].update(#
- {'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
- )
-
- if self.isRichards:
- # assuming that we normalised atmospheric pressure to zero,
- # pc = -pw in the Richards case.
- pc_prev = -p_prev_iter['wetting']
- else:
- pw_prev = p_prev_iter['wetting']
- pnw_prev = p_prev_iter['nonwetting']
- pc_prev = pnw_prev - pw_prev
-
- for phase in self.has_phases:
- # viscosity
- mu_a = mu[phase]
- # relative permeability
- ka = k[phase]
- # previous pressure iteration
- pa_prev = p_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:
- if debug:
- print(f"the following dofs (phase = {phase}) of interface{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
+ print(f"read gli_prev_iter[{subdomain}]:\n", gli_prev_of_neighbour)
+
+ # read previous neighbouring pressure values from interface and write them to
+ # the function pa_prev. We need to make a distinction first.
+ if neighbour_iter_num <= iteration:
+ # in case neighbour_iter_num == iteration the neighbour
+ # has not yet iterated beyond
+ # the iteration we are currently in, so
+ # interface.pressure_values holds the "previous pressure"
+ # we need to read from the interface to calculate gli.
+ # In case neighbour_iter_num < iteration the neighbour has
+ # already reached the error tolerance and has stopped iterating.
+ # We therefor read again the current and not the previous pressure
+ take_previous_pressure_from_interface = False
+ if neighbour_iter_num + 1 < iteration:
+ print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
+ 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().")
else:
- # 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.
- pass
- # 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
- # the wetting phase update is always to be calculated even if
- # only Richards is assumed on the current patch or the neighbour
- # only assumes the Richards equation.
-
- # 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, 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
- if debug:
- print(f"in iteration {iteration} and neighbour iteration = {neighbour_iter_num}:\n")
- print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
-
- interface.gli_term_prev[subdomain].update(#
- {'wetting': interface.gli_term[neighbour]['wetting'].copy()}
- )
- # print(colored(f"id(interface.gli_term_prev[{subdomain}][wetting])= {id(interface.gli_term_prev[subdomain]['wetting'])}","red"))
- # print(colored(f"id(interface.gli_term[{neighbour}][wetting])= {id(interface.gli_term[neighbour]['wetting'])}","red"))
- if debug:
- print(f"Interface{ind}.gli_term[{neighbour}]['wetting']=\n",interface.gli_term[neighbour]['wetting'])
- print(f"Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
-
- if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
- interface.gli_term_prev[subdomain].update(#
- {'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
- )
-
- for phase in self.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.
- if debug:
- print(f"Dict: Interface{ind}.gli_term[{neighbour}]['wetting']=\n",interface.gli_term[neighbour]['wetting'])
- print("should be the same as \n")
- print(f"From Dict Interface{ind}.gli_term_prev[{subdomain}]['wetting']=\n",interface.gli_term_prev[subdomain]['wetting'])
- 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
- )
- if debug:
- print(f"read gli_prev_iter[{subdomain}]:\n", gli_readout_tmp)
-
- # read neighbouring pressure values from interface and write them to
- # the function self.neighbouring_interface_pressure
- pa_prev = df.interpolate(df.Constant(0), self.function_space[phase])
- interface.write_pressure_dofs(to_function = pa_prev,#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
- # be overwritten with the gli_term of the neighbour.
- interface.gli_term_prev[subdomain].update(#
- {'wetting': interface.gli_term[neighbour]['wetting'].copy()}
- )
-
- if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
- interface.gli_term_prev[subdomain].update(#
- {'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
- )
-
- if neighbour_iter_num + 1 < iteration:
- print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
- 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
+ # the case we think here is neighbour_iter_num = iteration + 1.
+ # in this case the neighbour has iterated once beyond
+ # the iteration we are currently in, so
+ # interface.pressure_values_prev holds the "previous pressure"
+ # we need to read from the interface to calculate gli.
+ # neighbour_iter_num > iteration + 1 should not appear since
+ # we only stop to iterate if we reach error tolerance and no
+ # subdomain is left out while iterating.
+ take_previous_pressure_from_interface = True
+ if neighbour_iter_num > iteration + 1:
+ print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
+ f"and iteration = {iteration} on subdomain {subdomain}\n",
+ "This case should not be happening.")
+ pa_prev = df.interpolate(df.Constant(0), self.function_space[phase])
+ interface.write_pressure_dofs(to_function = pa_prev,#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 = take_previous_pressure_from_interface)
+
+ # 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_prev_part = gli_pressure_part.get_local()
+ gli = gli_p_prev_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
### END for ind in self.has_interface:
# save the result of the calculation to the subdomain and to the interface
@@ -738,23 +598,55 @@ class DomainPatch(df.SubDomain):
self.gli_assembled.update(
{phase: gli_assembled_tmp[phase].copy()}
)
+ # communicate the newly calculated terms
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
interface = self.interface[ind]
if debug:
print(f"Interface{ind}.gli_term[{subdomain}][{phase}]=\n",interface.gli_term[subdomain][phase])
print(f"before writing to interface[{ind}] gli_assembled[{phase}]=\n",self.gli_assembled[phase])
- # 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
- )
- if debug:
- print(f"gli_terms after writing to interface[{ind}] Interface[{ind}].gli_term[{subdomain}][{phase}]=\n",interface.gli_term[subdomain][phase])
+ # write gli dofs to interface dicts for communiction. We have to
+ # distinguish different cases.
+ if neighbour_iter_num == iteration:
+ # in this case, the neighbour has not done the next iteration
+ # and will need both pressure and gli_term_prev of the current
+ # subdomain. Therefor, we save the the calculated gli_terms to
+ # the gli_term dictionary of the current subdomain, not the
+ # gli_term_prev dictionary.
+ 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
+ )
+ if neighbour_iter_num -1 == iteration:
+ # in this case, the neighbour has done the next iteration already
+ # and will not need both pressure_prev and gli_term_prev of the current
+ # subdomain. Therefor, we save the the calculated gli_terms to
+ # the gli_term_prev dictionary of the current subdomain, not the
+ # gli_term dictionary.
+ 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 = True
+ )
+ # Additionally, the gli_term of the neighbour is saved to
+ # gli_term currently and will be overwritten in the next step,
+ # if we don't save it to gli_term_prev of the neighbour.
+ # similarly for the pressure values.
+ if phase == 'wetting' or not interface.neighbour_isRichards[subdomain]:
+ interface.gli_term_prev[neighbour].update(#
+ {phase: interface.gli_term[neighbour][phase].copy()}
+ )
+ interface.pressure_values_prev[neighbour][phase].update(
+ {phase: interface.pressure_values[neighbour][phase].copy()}
+ )
+ interface.current_iteration[subdomain] += 1
### END calc_gli_term
def null_all_interface_iteration_numbers(self) -> None:
--
GitLab