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

finish domainPatch.calc_gli_term method

parent a24455cf
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 135 additions and 58 deletions
domainPatch.py
+ 135
− 58
View file @ 715daf79
......@@ -209,6 +209,7 @@ class DomainPatch(df.SubDomain):
dx = self.dx
ds = self.ds
dt = self.timestep_size
porosity = self.porosity
if self.isRichards:
if phase == 'nonwetting':
print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
......@@ -241,7 +242,7 @@ class DomainPatch(df.SubDomain):
# gli = self.interface[index].gli_term[self.subdomain_index]['wetting']
# rhs_interface_forms.append((dt*gli*phi_w)*ds[index])
rhs1 = (Lw*pw_prev_iter*phi_w)*dx
rhs2 = -((S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
rhs2 = -(porosity*(S(pw_prev_iter) - S(pw_prev_timestep))*phi_w)*dx
rhs_without_gli = rhs1 + rhs2 + dt*source_w*phi_w*dx
form_and_rhs = {#
'form': form,#
......@@ -278,7 +279,7 @@ class DomainPatch(df.SubDomain):
form = form1 + form2 + df.sum(interface_forms)
rhs1 = (La*pw_prev_iter*phi_a)*dx
rhs2 = -((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
form_and_rhs = {#
......@@ -305,7 +306,7 @@ class DomainPatch(df.SubDomain):
rhs1 = (La*pnw_prev_iter*phi_a)*dx
# the non wetting phase has a + before instead of a - before the bracket
rhs2 = ((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
form_and_rhs = {#
......@@ -327,6 +328,9 @@ class DomainPatch(df.SubDomain):
ds form to be added to the rhs by the solver. This is due to the need of
having to communicate dofs over the interface. Additionally, the dofs are
being saved to an interface dictionary exactly as the pressures.
Caution: The array self.gli_assembled needs to be added to the rhs with
a minus sign by the solver!
"""
for ind in self.has_interface:
# self._calc_gli_term_on(interface, iteration)
......@@ -351,13 +355,13 @@ class DomainPatch(df.SubDomain):
if not neighbour_iter_num == 0:
# save the gl term from the neighbour first to use it in the
# next iteration
if self.isRichards:
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['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(#
{'wetting': interface.gli_term[neighbour]['wetting']}
)
else:
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting']},#
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
......@@ -436,58 +440,131 @@ class DomainPatch(df.SubDomain):
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:
if self.isRichards:
# first, save the gli_terms of the neighbour to gli_prev_term
# of the current subdomain.
# 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
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting']}
)
if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
interface.gli_term_prev[subdomain].update(#
{'wetting': interface.gli_term[neighbour]['wetting']}
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
# then read gli_prev term from interface and write it to
# the subdomain gli_assembled array
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 the neighbouring pressure values and write them to
# the placeholder fucntion 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 = subdomain,#
previous_iter = True)
pw_prev = pw_prev_iter['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)
interface.gli_term[subdomain].update({'wetting': gl0})
else:
mu_w = mu['wetting']
mu_nw = mu['nonwetting']
pw_prev = pw_prev_iter['wetting']
pnw_prev = pw_prev_iter['nonwetting']
kw = k['wetting']
knw = k['nonwetting']
#calc gl0 from the flux term
fluxw = -1/mu_w*kw(S(pnw_prev - pw_prev))*df.grad(pw_prev)
fluxnw = -1/mu_nw*knw(1-S(pnw_prev - pw_prev))*df.grad(pnw_prev)
gwl0 = (df.dot(fluxw, n) - Lambda['wetting']*pw_prev)
gnwl0 = (df.dot(fluxnw, n) - Lambda['nonwetting']*pnw_prev)
interface.gli_term[subdomain].update({'wetting': gwl0})
interface.gli_term[subdomain].update({'nonwetting': gnwl0})
# 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 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']}
)
if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
interface.gli_term_prev[subdomain].update(#
{'nonwetting': interface.gli_term[neighbour]['nonwetting']}
)
else:
print('uiae')
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.")
# finally we are done and can step up the current iteration count.
interface.current_iteration[subdomain] += 1
......@@ -999,10 +1076,10 @@ class Interface(BoundaryPart):
subdomain_ind: int,#
previous_iter: bool = False,#
) -> None:
""" write dofs of from_array with indices interface_dofs to self.gli_term
""" Read dofs of from_array with indices interface_dofs and write to self.gli_term
Write the degrees of freedom of the array from_array with indices
interface_dofs, i.e. the dofs of the interface, to one of the dictionaries
Read the degrees of freedom of the array from_array with indices
interface_dofs, i.e. the dofs of the interface, and write them to one of the dictionaries
of the interface
self.gli_term[subdomain_ind][phase], (if gl == True)
self.gli_term_prev[subdomain_ind][phase] (if previous_iter == True)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment