diff --git a/domainPatch.py b/domainPatch.py
index c92816d5d527b4759b1f74c9f2baca076986d4dc..44b7964cb3602dc87dfcab9d3cb3d09bc9fb7f90 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -190,16 +190,21 @@ class DomainPatch(df.SubDomain):
def govering_problem(self, phase: str, iter_num: int = 0) -> tp.Dict[str, fl.Form]:
- """ return the governing form ant right hand side as a dictionary"""
+ """ return the governing form and right hand side for phase phase and iteration
+ iter_num as a dictionary.
+
+ return the governing form ant right hand side for phase phase and iteration
+ iter_num (without the gli terms) as a dictionary.
+ """
# define measures
dx = self.dx
ds = self.ds
dt = self.timestep_size
if self.isRichards:
if phase == 'nonwetting':
- print("CAUTION: You invoked domainPatch.governing_problem() with \n#
- Parameter phase = 'nonwetting'. But isRichards = True for \n#
- current subdomain. \nReturning wetting phase equation only.\n")
+ print("CAUTION: You invoked domainPatch.governing_problem() with\n", #
+ "Parameter phase = 'nonwetting'. But isRichards = True for",
+ "current subdomain. \nReturning wetting phase equation only.\n")
Lw = self.L['wetting']
# this is pw_i in the paper
pw = self.pressure['wetting']
@@ -271,7 +276,7 @@ class DomainPatch(df.SubDomain):
rhs1 = (La*pw_prev_iter*phi_a)*dx
rhs2 = -((S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
- rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a)*dx
+ rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
form_and_rhs = {#
'form': form,#
@@ -283,7 +288,7 @@ class DomainPatch(df.SubDomain):
# gather interface forms
interface_forms = []
for interface in self.has_interface:
- if self.interface[interface].neighbour_is_Richards(self.subdomain_index):
+ if self.interface[interface].neighbour_isRichards[self.subdomain_index]:
# if the neighbour of our subdomain (with index self.subdomain_index)
# assumes a Richards model, we assume constant normalised atmospheric pressure
# and no interface term is practically appearing.
@@ -298,7 +303,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
- rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a)*dx
+ rhs_without_gli = rhs1 + rhs2 + dt*source_a*phi_a*dx
form_and_rhs = {#
'form': form,#
@@ -307,8 +312,9 @@ class DomainPatch(df.SubDomain):
return form_and_rhs
else:
- raise RuntimeError('missmatch for input parameter phase. \n#
- Expected either phase == wetting or phase == nonwetting ')
+ raise RuntimeError('missmatch for input parameter phase.',
+ 'Expected either phase == wetting or phase == nonwetting ')
+ return None
#### PRIVATE METHODS
def _init_function_space(self) -> None:
@@ -736,8 +742,7 @@ class Interface(BoundaryPart):
# this is an array of bools
isRichards: np.array = None,#
**kwds):
- # domain_index should be contained in adjacent_subdomains.
- # self.domain_index = domain_index
+ #
self.isRichards = isRichards
self.adjacent_subdomains = adjacent_subdomains
# pressure_values is a dictionary of dictionaries, one for each of the
@@ -755,6 +760,13 @@ class Interface(BoundaryPart):
self.adjacent_subdomains[0]: 0,#
self.adjacent_subdomains[1]: 0
}
+ # dictionary holding a bool variable for each subdomain index in self.adjacent_subdomains
+ # indicating wether the neighbouring subdomain w.r.t. that index assumes
+ # Richards model or not.
+ self.neighbour_isRichards = {
+ self.adjacent_subdomains[0]: self.isRichards[self.adjacent_subdomains[1]],#
+ self.adjacent_subdomains[1]: self.isRichards[self.adjacent_subdomains[0]]
+ }
# dictionary holding the neighbouring domain index relative to the key index
# this is used for calculating the gl_terms.
self.neighbour = {
@@ -797,7 +809,7 @@ class Interface(BoundaryPart):
def init_interface_values(self,#
interface_marker: df.MeshFunction,#
- interface_marker_value: int) -> dict:
+ interface_marker_value: int) -> None:
""" allocate dictionaries that are used to store the interface dof values
of the pressure, previous pressure, gl decoupling term as well as previous
gl update term respectively as pairs(parent_mesh_index: dof_value)
@@ -835,15 +847,15 @@ class Interface(BoundaryPart):
for phase in ['wetting', 'nonwetting']:
self.pressure_values[subdom_ind].update({phase: dict()})
self.pressure_values_prev[subdom_ind].update({phase: dict()})
- # the gli_term will hold forms not functions, therefor we leave
- # the dictionary empty. it will be populated by domainPatch._calc_gli_term()
- # self.gli_term[subdom_ind].update({phase: dict()})
- # self.gli_term_prev[subdom_ind].update({phase: dict()})
+ # the gli_term will hold dofs of the gli terms for communications
+ # they are initiated with 0
+ self.gli_term[subdom_ind].update({phase: dict()})
+ self.gli_term_prev[subdom_ind].update({phase: dict()})
for vert_ind in vertex_indices:
self.pressure_values[subdom_ind][phase].update({vert_ind: 0.0})
self.pressure_values_prev[subdom_ind][phase].update({vert_ind: 0.0})
- # self.gli_term[subdom_ind][phase].update({vert_ind: 0.0})
- # self.gli_term_prev[subdom_ind][phase].update({vert_ind: 0.0})
+ self.gli_term[subdom_ind][phase].update({vert_ind: 0.0})
+ self.gli_term_prev[subdom_ind][phase].update({vert_ind: 0.0})
# end init_interface_values
def write_dofs(self, from_function: df.Function, #