Skip to content
Snippets Groups Projects
Commit 66da3ba8 authored by David Seus's avatar David Seus
Browse files

fix implementation of calc_gli_term

parent 6b0d09a7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
LDDsimulation/domainPatch.py
+ 185
− 162
View file @ 66da3ba8
......@@ -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.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment