diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index d46d817cbd170b1009bbd33535fbce355e4f2cae..4699e0e493d6450c73ebca7f2515b2fcffbb89f5 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -535,6 +535,13 @@ class LDDsimulation(object):
subdomain.pressure_prev_iter[phase].assign(
subdomain.pressure[phase]
)
+ # onyl becomes active in TPR case
+ if self.include_gravity and subdomain.isRichards and subdomain.TPR_occures:
+ for interface in subdomain.has_interface_with_TP_neighbour:
+ subdomain.calc_gravity_neumann_flux(
+ interface_index=interface,
+ phase="nonwetting"
+ )
# end loop over subdomains
##### stopping criterion for the solver.
total_subsequent_error = np.amax(max_subsequent_error)
@@ -645,6 +652,17 @@ class LDDsimulation(object):
subdomain.write_pressure_to_interfaces()
subdomain.write_pressure_to_interfaces(previous_iter=True)
subdomain.calc_gl0_term()
+ # If gravity is included and a Richards subdomain has a TPR coupling,
+ # we need to calculate the Neumann flux of the nonwetting phase
+ # that appears due to gravity
+ if self.include_gravity:
+ if subdomain.isRichards and subdomain.TPR_occures:
+ for interface in subdomain.has_interface_with_TP_neighbour:
+ subdomain.calc_gravity_neumann_flux(
+ initial=True,
+ interface_index=interface,
+ phase="nonwetting"
+ )
return solution_over_iteration_within_timestep
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 386361676af49a94a96c89b724b0bc6cd5a104d3..7d1a962d09fe354abd30d4814108479472acd430 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -266,6 +266,8 @@ class Interface(BoundaryPart):
# dictionaries holding the forms for the decoupling terms
self.gli_term = dict()
self.gli_term_prev = dict()
+ self.gravity_neumann_flux = dict()
+ self.gravity_neumann_flux_prev = dict()
# dictionary holding the current iteration number i for each of the adjacent
# subdomains with index in adjacent_subdomains
self.current_iteration = {#
@@ -613,6 +615,8 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term.update({subdomain_index: dict()})
self.gli_term_prev.update({subdomain_index: dict()})
+ self.gravity_neumann_flux.update({subdomain_index: dict()})
+ self.gravity_neumann_flux_prev.update({subdomain_index: dict()})
for phase in has_phases:
if function == "pressure":
@@ -621,7 +625,8 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term[subdomain_index].update({phase: dict()})
self.gli_term_prev[subdomain_index].update({phase: dict()})
-
+ self.gravity_neumann_flux[subdomain_index].update({phase: dict()})
+ self.gravity_neumann_flux_prev[subdomain_index].update({phase: dict()})
space = function_space[function][phase]
dof_coordinates = self.dofs_and_coordinates(space, interface_marker, interface_marker_value)
# the dof dictionaries are grouped by the global (with respect to the
@@ -635,6 +640,8 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term[subdomain_index][phase].update({global_mesh_facet_index: dict()})
self.gli_term_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+ self.gravity_neumann_flux[subdomain_index][phase].update({global_mesh_facet_index: dict()})
+ self.gravity_neumann_flux_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
for dof_coord in facet_dof_coord_dict.values():
# we want to use the coordinates as key in a dictionary.
@@ -648,6 +655,8 @@ class Interface(BoundaryPart):
else: #function == "gli":
self.gli_term[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
self.gli_term_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+ self.gravity_neumann_flux[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
+ self.gravity_neumann_flux_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
# end init_interface_values
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 061f36e4d4408effd17566a9ead25fb5e66b0e4d..dc8be6a1e1316ff44349f61bf53773a9b19336d1 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -169,20 +169,34 @@ class DomainPatch(df.SubDomain):
### Class variables set by methods
# sometimes we need to loop over the phases and the length of the loop is
# determined by the model we are assuming
- if self.isRichards:
- self.has_phases = ['wetting']
- # in case we are Richards, determin, whether or not we have a TP
- # neighbour
- self.has_TP_neighbour = False
- for interf in self.has_interface:
- neighbour_assumes_R = self.interface[interf].neighbour_isRichards[self.subdomain_index]
- if not neighbour_assumes_R:
- self.has_TP_neighbour = True
- # if only one neighbour assumes TP we ste up the nonwetting
- # space.
- break
- else:
- self.has_phases = ['wetting', 'nonwetting'] #['nonwetting', 'wetting']
+ # flag that tells whether or not the current subdomain has a two-phase
+ # Richards coupling going on.
+ self.TPR_occures = False
+ if self.has_interface is not None:
+ if self.isRichards:
+ self.has_interface_with_Richards_neighbour = None
+ self.has_interface_with_TP_neighbour = []
+ self.has_phases = ['wetting']
+ # in case we are Richards, determine, whether or not we have a TP
+ # neighbour
+ for interf in self.has_interface:
+ neighbour_assumes_R = self.interface[interf].neighbour_isRichards[self.subdomain_index]
+ if not neighbour_assumes_R:
+ self.TPR_occures = True
+ self.has_interface_with_TP_neighbour.update(interf)
+ else:
+ self.has_interface_with_TP_neighbour = None
+ self.has_phases = ['wetting', 'nonwetting'] #['nonwetting', 'wetting']
+ # in case we assume TP we need to know, which interfaces have a
+ # Richards subdomain as neighbour.
+ self.has_interface_with_Richards_neighbour = []
+ for interface_index in self.has_interface:
+ interface = self.interface[interface_index]
+ if interface.neighbour_isRichards[self.subdomain_index]:
+ self.TPR_occures = True
+ self.has_interface_with_Richards_neighbour.append(interface_index)
+
+
# dictionary holding the dof indices corresponding to an interface of
# given interface. see self._calc_dof_indices_of_interfaces()
self._dof_indices_of_interface = dict()
@@ -230,17 +244,18 @@ class DomainPatch(df.SubDomain):
self.gravity_term = None
if self.has_interface is not None:
+
self._calc_interface_dof_indices_and_coordinates()
self._calc_corresponding_dof_indices()
self._init_interface_values()
- # Dictionary of FEM function holding the same dofs gli_assembled above
- # but as a df.Function. This is needed for writing out the gli terms when
- # analyising the convergence of a timestep.
- self.gli_function = dict()
- for phase in self.has_phases:
- self.gli_function.update(
- {phase: df.Function(self.function_space["gli"][phase])}
- )
+ # # Dictionary of FEM function holding the same dofs gli_assembled above
+ # # but as a df.Function. This is needed for writing out the gli terms when
+ # # analyising the convergence of a timestep.
+ # self.gli_function = dict()
+ # for phase in self.has_phases:
+ # self.gli_function.update(
+ # {phase: df.Function(self.function_space["gli"][phase])}
+ # )
# END constructor
@@ -358,6 +373,14 @@ class DomainPatch(df.SubDomain):
# gather interface forms
interface_forms = []
gli_forms = []
+ if self.include_gravity and phase=="nonwetting" and self.TPR_occures:
+ rhs_gravity_term_interface = []
+ for interface in self.has_interface_with_Richards_neighbour:
+ marker_value = self.interface[interface].marker_value
+ gravity_neumann_flux = self.read_gravity_neumann_flux(interface_index=interface, phase="nonwetting")
+ rhs_gravity_term_interface.append(
+ -dt*gravity_neumann_flux*phi_a*ds(marker_value)
+ )
# the marker values of the interfacese are global interfaces index + 1
for interface in self.has_interface:
# print(f"subdomain{self.subdomain_index} has interfaces: {interface}\n")
@@ -366,6 +389,7 @@ class DomainPatch(df.SubDomain):
interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
gli = self.calc_gli_term(interface_index=interface, phase=phase)
gli_forms.append((dt*gli*phi_a)*ds(marker_value))
+
# if our subdomain assumes Richards and our neighbour TP, we need to
# call self.calc_gli_term for the nonwetting phase even though we
# are not using it to determin the nonwetting pressure on this patch.
@@ -373,7 +397,7 @@ class DomainPatch(df.SubDomain):
# communication to the interface dictionries and
# the neighbour needs to read this updated gli term,
# to approximate the zero pressure.
- if self.isRichards:
+ if self.isRichards and self.TPR_occures:
# note, that in this case we know that phase == "wetting",
# since the solver runs only over phases beloning to our
# subdomain. But if self.isRichard, self.has_phases = ["wetting"]
@@ -401,8 +425,12 @@ class DomainPatch(df.SubDomain):
if self.has_interface is not None:
form = form1 + form2 + sum(interface_forms)
if self.include_gravity:
- # z_a included the minus sign already, see self._calc_gravity_expressions
- rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
+ if phase=="nonwetting" and self.TPR_occures:
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) \
+ + rhs_gravity + sum(rhs_gravity_term_interface)
+ else:
+ # z_a included the minus sign already, see self._calc_gravity_expressions
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
else:
rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
else:
@@ -419,6 +447,198 @@ class DomainPatch(df.SubDomain):
return form_and_rhs
+ def calc_gravity_neumann_flux(self,
+ interface_index: int,
+ phase: str,
+ initial: bool = False,
+ debug: bool = False
+ ) -> None:
+ """calculate the gravity term of the neumann flux
+
+ This method is written such that it can calculate the gravity neumann
+ flux for both phases in case it is needed, but is actually only used
+ for the flux of the nonwetting gravity term of a Richards domain with
+ TPR coupling.
+ """
+ subdomain = self.subdomain_index
+ dt = self.timestep_size
+ Lambda = self.lambda_param
+ interf_ind = interface_index
+
+ interface = self.interface[interf_ind]
+
+ if initial:
+ pressure_previous = self.pressure_prev_timestep
+ else:
+ pressure_previous = self.pressure_prev_iter
+
+ if self.isRichards:
+ # assuming that we normalised atmospheric pressure to zero,
+ # pc = -pw in the Richards case.
+ pc_prev_iter = -pressure_previous['wetting']
+ else:
+ pw_prev_iter = pressure_previous['wetting']
+ pnw_prev_iter = pressure_previous['nonwetting']
+ pc_prev_iter = pnw_prev - pw_prev
+
+ # n = df.FacetNormal(self.mesh)
+ S = self.saturation
+
+ # viscosity
+ mu_a = self.viscosity[phase]
+ # relative permeability
+ ka = self.relative_permeability[phase]
+ z_a = self.gravity_term[phase]
+
+ if phase == 'wetting':
+ # z_a included the minus sign, see self._calc_gravity_expressions
+ flux = 1/mu_a*ka(S(pc_prev_iter))*df.grad(z_a)
+ else:
+ # phase == 'nonwetting'
+ flux = 1/mu_a*ka(1-S(pc_prev_iter))*df.grad(z_a)
+
+ flux_function = df.project(flux, self.function_space["flux"][phase])
+ flux_x, flux_y = flux_function.split()
+ gravity_neummann_flux = df.Function(self.function_space["gli"][phase])
+ # # get dictionaries of global dof indices and coordinates along
+ # # the interface. These dictionaries have the facet indices of
+ # # facets belonging to the interface as keys (with respect to
+ # # the submesh numbering) and the dictionary
+ # # containing pairs {global_dof_index: dof_coordinate} for all
+ # # dofs along the facet as values.
+ neumann_flux_dof_coords = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
+ # p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
+ # flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
+ # # the attributes local and global for facet numbering mean
+ # # the following: local facet index is the index of the facet
+ # # with respect to the numbering of the submesh belonging to
+ # # the subdomain. global facet index refers to the facet numbering
+ # # of the mesh of the whole computational domain.
+ local2global_facet = interface.local_to_global_facet_map[subdomain]
+ corresponding_dof_index = self.interface_corresponding_dof_index[interf_ind][phase]
+ for local_facet_ind, normal in interface.outer_normals[subdomain].items():
+ neumann_flux_dofs_along_facet = neumann_flux_dof_coords[local_facet_ind]
+ global_facet_ind = local2global_facet[local_facet_ind]
+ for neumann_flux_dof_index, neumann_flux_dof_coord in neumann_flux_dofs_along_facet.items():
+ flux_x_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_x"]
+ flux_y_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_y"]
+
+ gravity_neummann_flux.vector()[neumann_flux_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+ + flux_y.vector()[flux_y_dof_index]*normal[1]
+
+ # save this newly calculated dof to the interface
+ # dictionary neumann_flux_save_dict
+ dof_coord_tupel = (neumann_flux_dof_coord[0], neumann_flux_dof_coord[1])
+ if initial:
+ # because we don't know if the Richards domain or the TP domain
+ # will be up next in the solver, we need to save the initial
+ # gravity neumann flux to both the current and the prev dictionary.
+ interface.gravity_neumann_flux_prev[subdomain][phase][global_facet_ind].update(
+ {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
+ )
+ interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
+ {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
+ )
+ else:
+ # in the noninitial case we always save to the current dictionary
+ # similarly to the pressures.
+ interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
+ {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
+ )
+
+
+ def read_gravity_neumann_flux(self,
+ interface_index: int,
+ phase: str,
+ debug: bool = False) -> df.Function: # tp.Dict[str, fl.Form]
+ """in the case of a TPR coupling read the additional gravity neumann flux
+ for the nonwetting phase from the neighbour
+
+ IMPLEMENTATION
+ same logic as in the calc_gli_term method applies.
+
+ PARAMETERS
+ ------------------------------------
+ interface_index int interface index for which to calculate the
+ gravity_neumann_flux term.
+ phase str phase to calculate the gravity_neumann_flux term for.
+ # DEBUG: bool toggle the debug output.
+ """
+ # this is the number of the current iteration of the subdomain.
+ iteration = self.iteration_number
+ subdomain = self.subdomain_index
+ interface = self.interface[interface_index]
+ # gravity_neumann_flux should be initialised by zero.
+ gravity_neumann_flux = df.Function(self.function_space["gli"][phase])
+
+ # neighbour index
+ neighbour = interface.neighbour[subdomain]
+ # update the current_iteration number of subdomain stored in the interface
+ # to the current iteration number of the subdomain.
+ interface.current_iteration[subdomain] = iteration
+ # save the iteration number of the neighbour
+ neighbour_iter_num = interface.current_iteration[neighbour]
+
+ # read previous neighbouring gravity neumann flux from interface.
+ # depending on our iteration number and the iteration number of
+ # the neighbour, weed to either read from the previous pressure
+ # dictionary or the current one.
+ if iteration == neighbour_iter_num + 1:
+ # in this case the neighbour has not yet iterated beyond
+ # the iteration we are currently in, so
+ # interface.gravity_neumann_flux[neighbour] holds the "previous neumann gravity flux"
+ # we need to read these values from the interface.
+ gravity_flux_read_dict = interface.gravity_neumann_flux[neighbour][phase]
+
+ else:
+ # iteration == neighbour_iter_num:
+ # this is the only possible other case. In this case the
+ # neighbour has iterated once beyond the iteration we are
+ # currently in, so interface.gravity_neumann_flux_prev holds the "previous neumann gravity flux"
+ # we need to read from the interface.
+ gravity_flux_read_dict = interface.gravity_neumann_flux_prev[neighbour][phase]
+
+ flux_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
+ local2global_facet = interface.local_to_global_facet_map[subdomain]
+ for local_facet_index, flux_facet_dofs2coordinates in flux_interface_dofs.items():
+ global_facet_index = local2global_facet[local_facet_index]
+ # read gravity_neumann_flux_prev term from the neighbour.
+ neighbour_flux_on_facet = gravity_flux_read_dict[global_facet_index]
+
+ for flux_dof_index, flux_dof_coord in flux_facet_dofs2coordinates.items():
+ # same with the gravity_neumann_flux_previous_dof
+ gravity_neumann_flux_neighbour_dof = None
+ for flux_neighbour_coord_tupel, flux_neighbour_dof in neighbour_flux_on_facet.items():
+ flux_neighbour_coord = np.array([flux_neighbour_coord_tupel[0], flux_neighbour_coord_tupel[1]])
+ # due to the warning in the documentation of np.allclose
+ # that np.allclose is not symetric, we test if both
+ # np.allclose(a,b) and np.allclose(b,a) are true.
+ a_close_b = np.allclose(flux_dof_coord, flux_neighbour_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(flux_neighbour_coord, flux_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
+ gravity_neumann_flux_neighbour_dof = flux_neighbour_dof
+ break
+
+ if gravity_neumann_flux_neighbour_dof is None:
+ raise RuntimeError(f"didn't find gravity_neumann_flux_neighbour_dof corresponding to dict key ({flux_dof_coord})",
+ "something is wrong")
+
+ gravity_neumann_flux.vector()[flux_dof_index] = gravity_neumann_flux_neighbour_dof
+
+ # In the following case, the gravity_neumann_flux of the neighbour is saved to
+ # gravity_neumann_flux currently and will be overwritten in the next step,
+ # if we don't save it to gravity_neumann_flux_prev of the neighbour.
+ if iteration == neighbour_iter_num:
+ # print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
+ # f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
+ interface.gravity_neumann_flux_prev[neighbour][phase].update(#
+ {global_facet_index: interface.gravity_neumann_flux[neighbour][phase][global_facet_index].copy()}
+ )
+
+ return gravity_neumann_flux
+ ### END read_gravity_neumann_flux
+
+
def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
"""calculate the gl0 terms for all interfaces of the subdomain.
@@ -449,10 +669,10 @@ class DomainPatch(df.SubDomain):
S = self.saturation
phases_for_gl0 = self.has_phases
- if self.isRichards and self.has_TP_neighbour:
- interface_has_TP_neighbour = not interface.neighbour_isRichards[self.subdomain_index]
- if interface_has_TP_neighbour and self.include_gravity:
- phases_for_gl0 = ['wetting', 'nonwetting']
+ # if self.isRichards and self.TPR_occures:
+ # interface_has_TP_neighbour = not interface.neighbour_isRichards[self.subdomain_index]
+ # if interface_has_TP_neighbour and self.include_gravity:
+ # phases_for_gl0 = ['wetting', 'nonwetting']
for phase in phases_for_gl0:
@@ -461,10 +681,10 @@ class DomainPatch(df.SubDomain):
# relative permeability
ka = self.relative_permeability[phase]
# previous timestep pressure
- if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
- pa_prev = 0
- else:
- pa_prev = self.pressure_prev_timestep[phase]
+ # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+ # pa_prev = 0
+ # else:
+ pa_prev = self.pressure_prev_timestep[phase]
if self.include_gravity:
z_a = self.gravity_term[phase]
@@ -479,17 +699,17 @@ class DomainPatch(df.SubDomain):
flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
else:
# phase == 'nonwetting'
- if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
- if self.include_gravity:
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(-z_a)
- else:
- flux = 0
+ # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+ # if self.include_gravity:
+ # flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(-z_a)
+ # else:
+ # flux = 0
+ # else:
+ # we need anothre flux in this case
+ if self.include_gravity:
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
else:
- # we need anothre flux in this case
- if self.include_gravity:
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
- else:
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
flux_function = df.project(flux, self.function_space["flux"][phase])
flux_x, flux_y = flux_function.split()
@@ -518,13 +738,13 @@ class DomainPatch(df.SubDomain):
flux_y_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_y"]
p_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["pressure"]
- if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
- gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
- + flux_y.vector()[flux_y_dof_index]*normal[1]
- else:
- gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
- + flux_y.vector()[flux_y_dof_index]*normal[1] \
- - Lambda[phase]*pa_prev.vector()[p_dof_index]
+ # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
+ # gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+ # + flux_y.vector()[flux_y_dof_index]*normal[1]
+ # else:
+ gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+ + flux_y.vector()[flux_y_dof_index]*normal[1] \
+ - Lambda[phase]*pa_prev.vector()[p_dof_index]
if debug:
print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_index]}")
@@ -727,6 +947,7 @@ class DomainPatch(df.SubDomain):
# print(f"np.allclose({p_coord}, {gli_dof_coord}, rtol=1e-10, atol=1e-14): {b_close_a}")
if a_close_b and b_close_a:
p_previous_dof = p_prev_dof
+ break
if p_previous_dof is None:
raise RuntimeError(f"didn't find p_previous_dof corresponding to dict key ({gli_dof_coord})",
@@ -742,6 +963,7 @@ class DomainPatch(df.SubDomain):
b_close_a = np.allclose(gli_prev_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
if a_close_b and b_close_a:
gli_previous_dof = gli_prev_dof
+ break
if gli_previous_dof is None:
raise RuntimeError(f"didn't find gli_previous_dof corresponding to dict key ({gli_dof_coord})",
"something is wrong")
@@ -809,7 +1031,7 @@ class DomainPatch(df.SubDomain):
gli_space = df.FunctionSpace(self.mesh, 'DG', degree)
flux_space = df.VectorFunctionSpace(self.mesh, 'DG', degree)
- if self.isRichards and self.has_TP_neighbour:
+ if self.isRichards and self.TPR_occures:
subdom_has_phases = ["wetting", "nonwetting"]
else:
subdom_has_phases = self.has_phases
@@ -1153,7 +1375,7 @@ class DomainPatch(df.SubDomain):
""" calculate the gravity term if it is included"""
g = df.Constant(self.gravity_acceleration) #
has_phases = self.has_phases
- if self.isRichards and self.has_TP_neighbour:
+ if self.isRichards and self.TPR_occures:
has_phases = ["wetting", "nonwetting"]
for phase in has_phases: