diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 3d1716be7acb17fa0b930e6e3133e0d5e2e4505b..c16233ac73295022f532c3748cac9ab55984579d 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -681,6 +681,8 @@ class LDDsimulation(object):
tmp_exact_pressure = df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])
tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
file.write(tmp_exact_pressure, t)
+
+ file.close()
t += self.timestep_size
def write_subsequent_errors_to_csv(self,#
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 46af28ac033f73133cbb9195ec9a63cff7fe3abc..57e5a220d4d962f17ace7e41dad9dac042bf6209 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -331,7 +331,7 @@ class Interface(BoundaryPart):
def dofs_and_coordinates(self, function_space: df.FunctionSpace,#
interface_marker: df.MeshFunction, #
interface_marker_value: int,
- debug: bool = True) -> tp.Dict[int, tp.Dict[int, np.ndarray]]:
+ debug: bool = False) -> tp.Dict[int, tp.Dict[int, np.ndarray]]:
""" for a given function space function_space, calculate the dof indices
and coordinates along the interface.
@@ -363,7 +363,7 @@ class Interface(BoundaryPart):
# this is the facet index relative to the Mesh
# that interface_marker is defined on.
facet_index = mesh_entity.index()
-
+
dofs_and_coords.update({facet_index: dict()})
# because mesh_entity doesn't have a method to access coordinates we
# need to cast it into an acutal facet object.
@@ -402,7 +402,7 @@ class Interface(BoundaryPart):
interface_marker: df.MeshFunction,#
interface_marker_value: int,
subdomain_index: int,
- debug: bool = True) -> None:
+ debug: bool = False) -> None:
""" create dictionary for each index subdomain_index containing the
index of each facet belonging to the interface (with respect to the mesh of
the subdomain subdomain_index) as key and stores a
@@ -467,7 +467,7 @@ class Interface(BoundaryPart):
print(f"outer normals per facet are: ", self.outer_normals[subdomain_index].items())
- def init_facet_maps(self, subdomain_index: int, debug: bool = True) -> None:
+ def init_facet_maps(self, subdomain_index: int, debug: bool = False) -> None:
""" calculate the mapping that maps the local numbering of interface
facets (w.r.t. the mesh of subdomain subdomain_index) to the global
index with respect to the global domain (subdomain 0)"""
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index 3a1af669afbced9982d570056c8d6fb7bd1a884d..439953b48bbe7c72eadb07e62b7b11c5a63e0266 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -86,10 +86,10 @@ isRichards = {
}
##################### MESH #####################################################
-mesh_resolution = 4
+mesh_resolution = 20
##################### TIME #####################################################
timestep_size = 0.01
-number_of_timesteps = 20
+number_of_timesteps = 200
# 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 = 0
diff --git a/RR-multi-patch-plus-gravity-const-solution/RR-multi-patch-with-gravity-constant-solution.py b/RR-multi-patch-plus-gravity-const-solution/RR-multi-patch-with-gravity-constant-solution.py
index 4ef083716191de8c2a4413d67b00e201b27613a5..e662d5bf3a96de94ba2732fa5d46ab4e28bb6cff 100755
--- a/RR-multi-patch-plus-gravity-const-solution/RR-multi-patch-with-gravity-constant-solution.py
+++ b/RR-multi-patch-plus-gravity-const-solution/RR-multi-patch-with-gravity-constant-solution.py
@@ -15,15 +15,15 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 29
+mesh_resolution = 40
# ----------------------------------------:-------------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.01
-number_of_timesteps = 1
+number_of_timesteps = 100
# 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 = 1
+number_of_timesteps_to_analyse = 11
starttime = 0
@@ -169,7 +169,7 @@ L = {
4: {'wetting': 0.25}
}
-lamdal_w = 30
+lamdal_w = 4
# subdom_num : lambda parameter for the L-scheme
lambda_param = {
1: {'wetting': lamdal_w},
diff --git a/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py b/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py
index b6e7bcb69f803c8882de204f053d57142f87ed1a..78d7cfcc272d1c3eb5b0f6af180cf101d9de67d4 100755
--- a/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py
+++ b/RR-multi-patch-plus-gravity/RR-multi-patch-with-gravity.py
@@ -20,7 +20,7 @@ mesh_resolution = 50
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.01
-number_of_timesteps = 500
+number_of_timesteps = 1000
# 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 = 11
@@ -163,10 +163,10 @@ porosity = {
# subdom_num : subdomain L for L-scheme
L = {
- 1: {'wetting': 0.5},
- 2: {'wetting': 0.5},
- 3: {'wetting': 0.5},
- 4: {'wetting': 0.5}
+ 1: {'wetting': 0.25},
+ 2: {'wetting': 0.25},
+ 3: {'wetting': 0.25},
+ 4: {'wetting': 0.25}
}
lamdal_w = 30
diff --git a/RR-multi-patch-with-inner-patch-const-solution/RR-multi-patch-with-inner-patch-constant-solution.py b/RR-multi-patch-with-inner-patch-const-solution/RR-multi-patch-with-inner-patch-constant-solution.py
index bacc24dc79fd40936ab4cf91dd3a255f5f912d43..ad5332aabf735695d93c5f5e9bb5f7aa805effaa 100755
--- a/RR-multi-patch-with-inner-patch-const-solution/RR-multi-patch-with-inner-patch-constant-solution.py
+++ b/RR-multi-patch-with-inner-patch-const-solution/RR-multi-patch-with-inner-patch-constant-solution.py
@@ -15,12 +15,12 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 40
+mesh_resolution = 20
# ----------------------------------------:-------------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.01
-number_of_timesteps = 30
+number_of_timesteps = 40
# 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
@@ -201,7 +201,7 @@ L = {
5: {'wetting': 0.25}
}
-lamdal_w = 32
+lamdal_w = 40
# subdom_num : lambda parameter for the L-scheme
lambda_param = {
@@ -330,25 +330,25 @@ def saturation_sym_prime(pc, index):
# incorporated in the construction of the exact solution below.
S_pc_sym = {
1: ft.partial(saturation_sym, index=1),
- 2: ft.partial(saturation_sym, index=2),
- 3: ft.partial(saturation_sym, index=2),
- 4: ft.partial(saturation_sym, index=2),
+ 2: ft.partial(saturation_sym, index=1),
+ 3: ft.partial(saturation_sym, index=1),
+ 4: ft.partial(saturation_sym, index=1),
5: ft.partial(saturation_sym, index=1)
}
S_pc_sym_prime = {
1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- 3: ft.partial(saturation_sym_prime, index=2),
- 4: ft.partial(saturation_sym_prime, index=2),
+ 2: ft.partial(saturation_sym_prime, index=1),
+ 3: ft.partial(saturation_sym_prime, index=1),
+ 4: ft.partial(saturation_sym_prime, index=1),
5: ft.partial(saturation_sym_prime, index=1)
}
sat_pressure_relationship = {
1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- 3: ft.partial(saturation, index=2),
- 4: ft.partial(saturation, index=2),
+ 2: ft.partial(saturation, index=1),
+ 3: ft.partial(saturation, index=1),
+ 4: ft.partial(saturation, index=1),
5: ft.partial(saturation, index=1)
}
diff --git a/TP-TP-layered-soil-case-const-solution/TP-TP-layered_soil-const-solution.py b/TP-TP-layered-soil-case-const-solution/TP-TP-layered_soil-const-solution.py
index 914239ae50b21012cc08632bb1496e2c5d35dea1..b17d1b80edc3f967bfee2248b304f2a27d9eaf69 100755
--- a/TP-TP-layered-soil-case-const-solution/TP-TP-layered_soil-const-solution.py
+++ b/TP-TP-layered-soil-case-const-solution/TP-TP-layered_soil-const-solution.py
@@ -27,7 +27,7 @@ mesh_resolution = 40
# ----------------------------------------:-------------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-timestep_size = 0.001
+timestep_size = 0.01
number_of_timesteps = 100
# decide how many timesteps you want analysed. Analysed means, that we write
# out subsequent errors of the L-iteration within the timestep.
diff --git a/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py b/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py
index 03c28f6bd464f155d9609f9c7f07725264a5f647..4f469ec029df53d19985c6ae0ff63952a0c52f8b 100755
--- a/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py
+++ b/TP-TP-layered-soil-case-with-inner-patch-constant-solution/TP-TP-layered_soil_with_inner_patch_const_solution.py
@@ -27,7 +27,7 @@ mesh_resolution = 41
# ----------------------------------------:-----------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-timestep_size = 0.001
+timestep_size = 0.00001
number_of_timesteps = 100
# decide how many timesteps you want analysed. Analysed means, that we write
# out subsequent errors of the L-iteration within the timestep.
@@ -37,8 +37,8 @@ starttime = 0
Lw = 0.25
Lnw = 0.25
-l_param_w = 40
-l_param_nw = 40
+l_param_w = 100
+l_param_nw = 100
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
@@ -519,98 +519,94 @@ ka_prime = {
# since we unify the treatment in the code for Richards and two-phase, we need
# the same requierment
# for both cases, two-phase and Richards.
-# def saturation(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
+def saturation(pc, n_index, alpha):
+ # inverse capillary pressure-saturation-relationship
+ return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
#
-# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
-# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
-# def saturation_sym(pc, n_index, alpha):
+# S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
+# we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
+def saturation_sym(pc, n_index, alpha):
+ # inverse capillary pressure-saturation-relationship
+ #df.conditional(pc > 0,
+ return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
+
+
+# derivative of S-pc relationship with respect to pc. This is needed for the
+# construction of a analytic solution.
+def saturation_sym_prime(pc, n_index, alpha):
+ # inverse capillary pressure-saturation-relationship
+ return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
+
+# note that the conditional definition of S-pc in the nonsymbolic part will be
+# incorporated in the construction of the exact solution below.
+S_pc_sym = {
+ 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+ 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+ 3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+ 4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+ 5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
+ 6: ft.partial(saturation_sym, n_index=3, alpha=0.001)
+}
+
+S_pc_sym_prime = {
+ 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+ 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+ 3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+ 4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+ 5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
+ 6: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
+}
+
+sat_pressure_relationship = {
+ 1: ft.partial(saturation, n_index=3, alpha=0.001),
+ 2: ft.partial(saturation, n_index=3, alpha=0.001),
+ 3: ft.partial(saturation, n_index=3, alpha=0.001),
+ 4: ft.partial(saturation, n_index=3, alpha=0.001),
+ 5: ft.partial(saturation, n_index=3, alpha=0.001),
+ 6: ft.partial(saturation, n_index=3, alpha=0.001)
+}
+
+# def saturation(pc, n_index):
# # inverse capillary pressure-saturation-relationship
-# #df.conditional(pc > 0,
-# return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
+# return df.conditional(pc > 0, 1/((1 + pc)**(1/(n_index + 1))), 1)
#
#
-# # derivative of S-pc relationship with respect to pc. This is needed for the
-# # construction of a analytic solution.
-# def saturation_sym_prime(pc, n_index, alpha):
+# def saturation_sym(pc, n_index):
# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
+# return 1/((1 + pc)**(1/(n_index + 1)))
+#
+# def saturation_sym_prime(pc, n_index):
+# # inverse capillary pressure-saturation-relationship
+# return -1/((n_index+1)*(1 + pc)**((n_index+2)/(n_index+1)))
#
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-
#
-# # note that the conditional definition of S-pc in the nonsymbolic part will be
-# # incorporated in the construction of the exact solution below.
# S_pc_sym = {
-# 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-# 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-# 3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-# 4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-# 5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
-# 6: ft.partial(saturation_sym, n_index=3, alpha=0.001)
+# 1: ft.partial(saturation_sym, n_index=1),
+# 2: ft.partial(saturation_sym, n_index=1),
+# 3: ft.partial(saturation_sym, n_index=1),
+# 4: ft.partial(saturation_sym, n_index=1),
+# 5: ft.partial(saturation_sym, n_index=1),
+# 6: ft.partial(saturation_sym, n_index=1)
# }
#
# S_pc_sym_prime = {
-# 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-# 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-# 3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-# 4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-# 5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
-# 6: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
+# 1: ft.partial(saturation_sym_prime, n_index=1),
+# 2: ft.partial(saturation_sym_prime, n_index=1),
+# 3: ft.partial(saturation_sym_prime, n_index=1),
+# 4: ft.partial(saturation_sym_prime, n_index=1),
+# 5: ft.partial(saturation_sym_prime, n_index=1),
+# 6: ft.partial(saturation_sym_prime, n_index=1)
# }
#
# sat_pressure_relationship = {
-# 1: ft.partial(saturation, n_index=3, alpha=0.001),
-# 2: ft.partial(saturation, n_index=3, alpha=0.001),
-# 3: ft.partial(saturation, n_index=3, alpha=0.001),
-# 4: ft.partial(saturation, n_index=3, alpha=0.001),
-# 5: ft.partial(saturation, n_index=3, alpha=0.001),
-# 6: ft.partial(saturation, n_index=3, alpha=0.001)
+# 1: ft.partial(saturation, n_index=1),
+# 2: ft.partial(saturation, n_index=1),
+# 3: ft.partial(saturation, n_index=1),
+# 4: ft.partial(saturation, n_index=1),
+# 5: ft.partial(saturation, n_index=1),
+# 6: ft.partial(saturation, n_index=1)
# }
-def saturation(pc, n_index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(n_index + 1))), 1)
-
-
-def saturation_sym(pc, n_index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(n_index + 1)))
-
-def saturation_sym_prime(pc, n_index):
- # inverse capillary pressure-saturation-relationship
- return -1/((n_index+1)*(1 + pc)**((n_index+2)/(n_index+1)))
-
-
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=1),
- 2: ft.partial(saturation_sym, n_index=1),
- 3: ft.partial(saturation_sym, n_index=1),
- 4: ft.partial(saturation_sym, n_index=1),
- 5: ft.partial(saturation_sym, n_index=1),
- 6: ft.partial(saturation_sym, n_index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=1),
- 2: ft.partial(saturation_sym_prime, n_index=1),
- 3: ft.partial(saturation_sym_prime, n_index=1),
- 4: ft.partial(saturation_sym_prime, n_index=1),
- 5: ft.partial(saturation_sym_prime, n_index=1),
- 6: ft.partial(saturation_sym_prime, n_index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=1),
- 2: ft.partial(saturation, n_index=1),
- 3: ft.partial(saturation, n_index=1),
- 4: ft.partial(saturation, n_index=1),
- 5: ft.partial(saturation, n_index=1),
- 6: ft.partial(saturation, n_index=1)
-}
-
#############################################
# Manufacture source expressions with sympy #
@@ -618,8 +614,8 @@ sat_pressure_relationship = {
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
-p_dict = {'wetting': -3.0 + 0.0*t,
- 'nonwetting': -1.0 + 0.0*t}
+p_dict = {'wetting': -3 + 0*t,
+ 'nonwetting': -1 + 0*t}
p_e_sym = dict()
pc_e_sym = dict()
for subdomain, isR in isRichards.items():