diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index d7bf0ec4f787af51e3798dec713b0c813d9b4d5f..061f36e4d4408effd17566a9ead25fb5e66b0e4d 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -171,6 +171,16 @@ class DomainPatch(df.SubDomain):
# 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']
# dictionary holding the dof indices corresponding to an interface of
@@ -437,15 +447,25 @@ class DomainPatch(df.SubDomain):
# n = df.FacetNormal(self.mesh)
S = self.saturation
- for phase in self.has_phases:
+
+ 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']
+
+ for phase in phases_for_gl0:
+
# viscosity
mu_a = self.viscosity[phase]
# relative permeability
ka = self.relative_permeability[phase]
# previous timestep pressure
- pa_prev = self.pressure_prev_timestep[phase]
- # testfunction
- # phi_a = self.testfunction[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]
if phase == 'wetting':
@@ -459,11 +479,17 @@ class DomainPatch(df.SubDomain):
flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
else:
# phase == 'nonwetting'
- # 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)
+ 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:
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+ # 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_function = df.project(flux, self.function_space["flux"][phase])
flux_x, flux_y = flux_function.split()
@@ -492,9 +518,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"]
- 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]}")
@@ -779,18 +809,8 @@ class DomainPatch(df.SubDomain):
gli_space = df.FunctionSpace(self.mesh, 'DG', degree)
flux_space = df.VectorFunctionSpace(self.mesh, 'DG', degree)
- if self.isRichards:
- # if our domain assumes Richards but the neighbour of the interface
- # assumes TP, we need to initialise interface values for the nonwetting
- # phase, too, for communication, see the article. Therefor we need to
- # initialise the nonwetting function space too.
- for interf in self.has_interface:
- neighbour_assumes_R = self.interface[interf].neighbour_isRichards[self.subdomain_index]
- if not neighbour_assumes_R:
- subdom_has_phases = ["wetting", "nonwetting"]
- # if only one neighbour assumes TP we ste up the nonwetting
- # space.
- break
+ if self.isRichards and self.has_TP_neighbour:
+ subdom_has_phases = ["wetting", "nonwetting"]
else:
subdom_has_phases = self.has_phases
@@ -1132,7 +1152,11 @@ class DomainPatch(df.SubDomain):
def _calc_gravity_expressions(self):
""" calculate the gravity term if it is included"""
g = df.Constant(self.gravity_acceleration) #
- for phase in self.has_phases:
+ has_phases = self.has_phases
+ if self.isRichards and self.has_TP_neighbour:
+ has_phases = ["wetting", "nonwetting"]
+
+ for phase in has_phases:
rho = df.Constant(self.densities[phase])
self.gravity_term.update(
{phase: df.Expression('-g*rho*x[1]',
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index cd8315d0ebd815df1a69cd2a8a487c09554e480d..6d3495afc0d06006d432a9665818fdb6df2628da 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -38,21 +38,21 @@ resolutions = {
############ GRID #######################
# mesh_resolution = 20
-timestep_size = 0.001
-number_of_timesteps = 20
-plot_timestep_every = 1
+timestep_size = 0.0001
+number_of_timesteps = 1200
+plot_timestep_every = 5
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
-number_of_timesteps_to_analyse = 10
+number_of_timesteps_to_analyse = 6
starttimes = [0.0]
-Lw = 0.025 #/timestep_size
+Lw = 10 #/timestep_size
Lnw=Lw
-lambda_w = 40
-lambda_nw = 40
+lambda_w = 4
+lambda_nw = 4
-include_gravity = False
+include_gravity = True
debugflag = True
analyse_condition = True
@@ -151,8 +151,8 @@ isRichards = {
viscosity = {#
# subdom_num : viscosity
- 1 : {'wetting' :1},
- #'nonwetting': 1}, #
+ 1 : {'wetting' :1,
+ 'nonwetting': 1/50}, #
2 : {'wetting' :1,
'nonwetting': 1/50}
}
@@ -165,7 +165,8 @@ porosity = {#
# Dict of the form: { subdom_num : density }
densities = {
- 1: {'wetting': 997},
+ 1: {'wetting': 997,
+ 'nonwetting': 1.225},
2: {'wetting': 997,
'nonwetting': 1.225},
}
@@ -194,15 +195,15 @@ def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
return s**2
-# def rel_perm1nw(s):
-# # relative permeabilty nonwetting on subdomain1
-# return (1-s)**2
+def rel_perm1nw(s):
+ # relative permeabilty nonwetting on subdomain1
+ return (1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
-# _rel_perm1nw = ft.partial(rel_perm1nw)
+_rel_perm1nw = ft.partial(rel_perm1nw)
subdomain1_rel_perm = {
'wetting': _rel_perm1w,#
- # 'nonwetting': _rel_perm1nw
+ 'nonwetting': _rel_perm1nw
}
## relative permeabilty functions on subdomain 2
def rel_perm2w(s):
@@ -233,9 +234,9 @@ def rel_perm1w_prime(s):
# relative permeabilty on subdomain1
return 2*s
-# def rel_perm1nw_prime(s):
-# # relative permeabilty on subdomain1
-# return 2*(1-s)
+def rel_perm1nw_prime(s):
+ # relative permeabilty on subdomain1
+ return -2*(1-s)
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
@@ -248,13 +249,13 @@ def rel_perm2nw_prime(s):
return -3*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-# _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
+_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime
- # 'nonwetting': _rel_perm1nw_prime
+ 'wetting': _rel_perm1w_prime,
+ 'nonwetting': _rel_perm1nw_prime
}
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
index bc102a36cd207b46c15f9d552bbfebfbe18c15ea..f3b556c5c1d7257eddd6db84a344e0e4c085c513 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
@@ -19,7 +19,7 @@ datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
-use_case = "TP-R-2-patch-realistic"
+use_case = "TP-R-2-patch-realistic-new-nw-gli"
# solver_tol = 6E-7
max_iter_num = 1000
FEM_Lagrange_degree = 1
@@ -32,12 +32,13 @@ resolutions = { 1: 7e-7,
32: 7e-7,
64: 7e-7,
128: 7e-7,
- 256: 7e-7}
+ #256: 7e-7,
+ }
############ GRID #######################
# mesh_resolution = 20
timestep_size = 0.001
-number_of_timesteps = 600
+number_of_timesteps = 800
plot_timestep_every = 2
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
@@ -182,8 +183,8 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :lambda_w},
- # 'nonwetting': lambda},#
+ 1 : {'wetting' :lambda_w,
+ 'nonwetting': lambda_nw},#
2 : {'wetting' :lambda_w,
'nonwetting': lambda_nw}
}