From b4538afe189c296f4e6546b95049227614190d2a Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Tue, 23 Jul 2019 17:21:57 +0200
Subject: [PATCH] vor dem urlaub
---
LDDsimulation/LDDsimulation.py | 7 +-
LDDsimulation/solutionFile.py | 2 +-
.../RR-multi-patch-with-inner-patch.py | 104 ++++++-
.../TP-R-2-patch-test-constant-solution.py | 150 ++++-----
TP-R-two-patch-test-case/TP-R-2-patch-test.py | 155 ++++------
.../TP-TP-2-patch-constant-solution.py | 285 ++++--------------
.../TP-TP-2-patch-pure-dd.py | 213 +++++++------
...ed_soil_with_inner_patch_const_solution.py | 4 +-
...h-with-gravity-same-wetting-phase-as-RR.py | 131 +++-----
TP-one-patch/TP-one-patch.py | 268 ++++++++++------
10 files changed, 611 insertions(+), 708 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 647574d..1c6d16e 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -54,6 +54,7 @@ class LDDsimulation(object):
# # To be able to run DG in parallel
# df.parameters["ghost_mode"] = "shared_facet"
# df.parameters["ghost_mode"] = "none"
+ df.parameters["mesh_partitioner"] = "ParMETIS"
# Form compiler options
df.parameters["form_compiler"]["optimize"] = True
df.parameters["form_compiler"]["cpp_optimize"] = True
@@ -103,7 +104,7 @@ class LDDsimulation(object):
## Private variables
# maximal number of L-iterations that the LDD solver uses.
- self._max_iter_num = 2000
+ self._max_iter_num = 500
# TODO rewrite this with regard to the mesh sizes
# self.calc_tol = self.tol
# list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps
@@ -161,8 +162,8 @@ class LDDsimulation(object):
self.solver_type_is_Kryov = True
### Define the linear solver to be used.
- self.solver = 'bicgstab' #'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
- self.preconditioner = 'jacobi' #'default'#jacobi#'hypre_amg' #'ilu'#'hypre_amg' # incomplete LU factorization
+ self.solver = 'bicgstab'#'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
+ self.preconditioner = 'jacobi' #'hypre_amg' 'default' #'ilu'#'hypre_amg' # incomplete LU factorization
# dictionary of solver parametrs. This is passed to self._init_subdomains,
# where for each subdomain a sovler object of type self.solver is created
# with these parameters.
diff --git a/LDDsimulation/solutionFile.py b/LDDsimulation/solutionFile.py
index 7ae0a95..e71450f 100644
--- a/LDDsimulation/solutionFile.py
+++ b/LDDsimulation/solutionFile.py
@@ -9,5 +9,5 @@ class SolutionFile(df.XDMFFile):
df.XDMFFile.__init__(self, mpicomm, filepath)
self.parameters["functions_share_mesh"] = True
- self.parameters["flush_output"] = False
+ self.parameters["flush_output"] = True
self.path = filepath # Mimic the file path attribute from a `file` returned by `open`
diff --git a/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py b/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py
index f9d20f8..d83d619 100755
--- a/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py
+++ b/RR-multi-patch-with-inner-patch/RR-multi-patch-with-inner-patch.py
@@ -20,12 +20,15 @@ mesh_resolution = 51
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.01
-number_of_timesteps = 160
+number_of_timesteps = 150
# 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
+number_of_timesteps_to_analyse = 10
starttime = 0
+# Lw = 0.25/timestep_size
+Lw = 0.25
+
# ----------------------------------------------------------------------------#
# ------------------- Domain and Interface -----------------------------------#
@@ -192,13 +195,14 @@ porosity = {
5: 1
}
+
# subdom_num : subdomain L for L-scheme
L = {
- 1: {'wetting': 0.25},
- 2: {'wetting': 0.25},
- 3: {'wetting': 0.25},
- 4: {'wetting': 0.25},
- 5: {'wetting': 0.25}
+ 1: {'wetting': Lw},
+ 2: {'wetting': Lw},
+ 3: {'wetting': Lw},
+ 4: {'wetting': Lw},
+ 5: {'wetting': Lw}
}
lamdal_w = 32
@@ -374,14 +378,88 @@ sat_pressure_relationship = {
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
+# cutoff function
+
+epsilon_x_inner = 0.7
+epsilon_x_outer = 0.99
+epsilon_y_inner = epsilon_x_inner
+epsilon_y_outer = epsilon_x_outer
+
+def mollifier(x, epsilon):
+ """ one d mollifier """
+ out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1)
+ return out_expr
+
+mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner)
+
+pw_sym_x = sym.Piecewise(
+ (mollifier_handle(x), x**2 < epsilon_x_outer**2),
+ (0, True)
+)
+pw_sym_y = sym.Piecewise(
+ (mollifier_handle(y), y**2 < epsilon_y_outer**2),
+ (0, True)
+)
+
+def mollifier2d(x, y, epsilon):
+ """ one d mollifier """
+ out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1)
+ return out_expr
+
+mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer)
+
+pw_sym2d_x = sym.Piecewise(
+ (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2),
+ (0, True)
+)
+
+zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise(
+ (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))),
+ (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise(
+ (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))),
+ (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise(
+ (1, y<=-2*epsilon_x_inner),
+ (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))),
+ (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y
+gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x
+cutoff = gaussian/(gaussian + zero_on_shrinking)
+
+
+epsilon = 0.5
+
+# *(1-(1-y)/epsilon)**2*(1-(1-x)/epsilon)**2*(1-(-1 - y)/epsilon)**2*(1-(-1 - x)/epsilon)**2
+# p_e_sym = {
+# 1: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
+# 2: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
+# 3: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
+# 4: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
+# 5: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2}
+# }
+
p_e_sym = {
- 1: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)},
- 2: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x)},
- 3: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x)},
- 4: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x)},
- 5: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)}
+ 1: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*cutoff},
+ 2: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*cutoff},
+ 3: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*cutoff},
+ 4: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*cutoff},
+ 5: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*cutoff}
}
+
pc_e_sym = {
1: -1*p_e_sym[1]['wetting'],
2: -1*p_e_sym[2]['wetting'],
@@ -570,7 +648,7 @@ write_to_file = {
# initialise LDD simulation class
simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-7)
-simulation.set_parameters(output_dir="./output/",
+simulation.set_parameters(output_dir="./output/with_cutoff_function",
subdomain_def_points=subdomain_def_points,
isRichards=isRichards,
interface_def_points=interface_def_points,
diff --git a/TP-R-two-patch-test-case-constant-solution/TP-R-2-patch-test-constant-solution.py b/TP-R-two-patch-test-case-constant-solution/TP-R-2-patch-test-constant-solution.py
index aab3e17..11dbcf5 100755
--- a/TP-R-two-patch-test-case-constant-solution/TP-R-2-patch-test-constant-solution.py
+++ b/TP-R-two-patch-test-case-constant-solution/TP-R-2-patch-test-constant-solution.py
@@ -7,11 +7,31 @@ import typing as tp
import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
+import helpers as hlp
#import ufl as ufl
# init sympy session
sym.init_printing()
+solver_tol = 5e-7
+############ GRID #######################ü
+mesh_resolution = 30
+timestep_size = 0.0002
+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 = 10
+starttime = 0
+
+Lw = 100/timestep_size
+Lnw=Lw
+
+l_param_w = 40
+l_param_nw = l_param_w
+
+include_gravity = False
+
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0, -1.0),
@@ -88,15 +108,6 @@ isRichards = {
}
-############ GRID #######################ü
-mesh_resolution = 20
-timestep_size = 0.001
-number_of_timesteps = 50
-# 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
-starttime = 0
-
viscosity = {#
# subdom_num : viscosity
1 : {'wetting' :1},
@@ -122,19 +133,19 @@ gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
- 1 : {'wetting' :0.25},
+ 1 : {'wetting' :Lw},
# 'nonwetting': 0.25},#
- 2 : {'wetting' :0.25,
- 'nonwetting': 0.25}
+ 2 : {'wetting' :Lw,
+ 'nonwetting': Lnw}
}
-l_param = 40
+
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :l_param},
+ 1 : {'wetting' :l_param_w},
# 'nonwetting': l_param},#
- 2 : {'wetting' :l_param,
- 'nonwetting': l_param}
+ 2 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw}
}
## relative permeabilty functions on subdomain 1
@@ -193,7 +204,7 @@ def rel_perm2w_prime(s):
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
# _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -281,81 +292,38 @@ p_e_sym = {
} #-(y-0.5)*(y-0.5)*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
-pc_e_sym = {
- 1: -1*p_e_sym[1]['wetting'],
- 2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting']
-}
-
-# turn above symbolic code into exact solution for dolphin and
-# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
+pc_e_sym = dict()
for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
if isR:
- subdomain_has_phases = ["wetting"]
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
+ - p_e_sym[subdomain]['wetting']})
+
+
+symbols = {"x": x,
+ "y": y,
+ "t": t}
+# turn above symbolic code into exact solution for dolphin and
+# construct the rhs that matches the above exact solution.
+exact_solution_example = hlp.generate_exact_solution_expressions(
+ symbols=symbols,
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym,
+ symbolic_capillary_pressure=pc_e_sym,
+ saturation_pressure_relationship=S_pc_sym,
+ saturation_pressure_relationship_prime=S_pc_sym_prime,
+ viscosity=viscosity,
+ porosity=porosity,
+ relative_permeability=relative_permeability,
+ relative_permeability_prime=ka_prime,
+ densities=densities,
+ gravity_acceleration=gravity_acceleration,
+ include_gravity=include_gravity,
+ )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
@@ -399,7 +367,7 @@ write_to_file = {
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol = 5E-4, debug = True)
+simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=False)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
@@ -422,7 +390,7 @@ simulation.set_parameters(output_dir = "./output/",#
dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
densities=densities,
- include_gravity=True,
+ include_gravity=include_gravity,
write2file = write_to_file,#
)
diff --git a/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 459f2dd..ef696e3 100755
--- a/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -7,11 +7,30 @@ import typing as tp
import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
+import helpers as hlp
#import ufl as ufl
# init sympy session
sym.init_printing()
+solver_tol = 5e-7
+############ GRID #######################ü
+mesh_resolution = 30
+timestep_size = 0.0001
+number_of_timesteps = 50
+# 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
+starttime = 0
+
+Lw = 1/timestep_size
+Lnw=Lw
+
+l_param_w = 40
+l_param_nw = l_param_w
+
+include_gravity = True
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0, -1.0),
@@ -80,15 +99,6 @@ isRichards = {
}
-############ GRID #######################ü
-mesh_resolution = 50
-timestep_size = 0.01
-number_of_timesteps = 160
-# 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
-starttime = 0
-
viscosity = {#
# subdom_num : viscosity
1 : {'wetting' :1},
@@ -114,19 +124,19 @@ gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
- 1 : {'wetting' :0.25},
+ 1 : {'wetting' :Lw},
# 'nonwetting': 0.25},#
- 2 : {'wetting' :0.25,
- 'nonwetting': 0.25}
+ 2 : {'wetting' :Lw,
+ 'nonwetting': Lnw}
}
-l_param = 40
+
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :l_param},
+ 1 : {'wetting' :l_param_w},
# 'nonwetting': l_param},#
- 2 : {'wetting' :l_param,
- 'nonwetting': l_param}
+ 2 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw}
}
## relative permeabilty functions on subdomain 1
@@ -185,7 +195,7 @@ def rel_perm2w_prime(s):
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return 3*(1-s)**2
+ return -3*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
# _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -307,87 +317,44 @@ x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym = {
- 1: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)},
- 2: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x),
- 'nonwetting': (-t*(1-y - x**2)**2 - sym.sqrt(2+t**2))*y},
+ 1: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*(1-x)**2*(1+x)**2*(1-y)**2},
+ 2: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*(1-x)**2*(1+x)**2*(1+y)**2,
+ 'nonwetting': (-t**2*(1+y + x**2)**2 - sym.sqrt(2+t**4))*y**2*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
-pc_e_sym = {
- 1: -1*p_e_sym[1]['wetting'],
- 2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting']
-}
-
-# turn above symbolic code into exact solution for dolphin and
-# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
+pc_e_sym = dict()
for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
if isR:
- subdomain_has_phases = ["wetting"]
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
+
+
+symbols = {"x": x,
+ "y": y,
+ "t": t}
+# turn above symbolic code into exact solution for dolphin and
+# construct the rhs that matches the above exact solution.
+exact_solution_example = hlp.generate_exact_solution_expressions(
+ symbols=symbols,
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym,
+ symbolic_capillary_pressure=pc_e_sym,
+ saturation_pressure_relationship=S_pc_sym,
+ saturation_pressure_relationship_prime=S_pc_sym_prime,
+ viscosity=viscosity,
+ porosity=porosity,
+ relative_permeability=relative_permeability,
+ relative_permeability_prime=ka_prime,
+ densities=densities,
+ gravity_acceleration=gravity_acceleration,
+ include_gravity=include_gravity,
+ )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
@@ -431,7 +398,7 @@ write_to_file = {
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol = 1E-7, debug = False)
+simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=True)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
@@ -454,7 +421,7 @@ simulation.set_parameters(output_dir = "./output/",#
dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
densities=densities,
- include_gravity=True,
+ include_gravity=include_gravity,
write2file = write_to_file,#
)
diff --git a/TP-TP-2-patch-constant-solution/TP-TP-2-patch-constant-solution.py b/TP-TP-2-patch-constant-solution/TP-TP-2-patch-constant-solution.py
index b60c4a7..0dfa344 100755
--- a/TP-TP-2-patch-constant-solution/TP-TP-2-patch-constant-solution.py
+++ b/TP-TP-2-patch-constant-solution/TP-TP-2-patch-constant-solution.py
@@ -7,11 +7,32 @@ import typing as tp
import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
+import helpers as hlp
#import ufl as ufl
# init sympy session
sym.init_printing()
+solver_tol = 5E-6
+
+############ GRID #######################ü
+mesh_resolution = 20
+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.
+number_of_timesteps_to_analyse = 10
+starttime = 0
+
+Lw = 1/timestep_size
+Lnw=Lw
+
+l_param_w = 40
+l_param_nw = 40
+
+include_gravity = True
+
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0,-1.0), #
@@ -80,15 +101,6 @@ isRichards = {
}
-############ GRID #######################ü
-mesh_resolution = 41
-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.
-number_of_timesteps_to_analyse = 11
-starttime = 0
-
viscosity = {#
# subdom_num : viscosity
1 : {'wetting' :1,
@@ -116,19 +128,19 @@ porosity = {#
L = {#
# subdom_num : subdomain L for L-scheme
- 1 : {'wetting' :0.25,
- 'nonwetting': 0.25},#
- 2 : {'wetting' :0.25,
- 'nonwetting': 0.25}
+ 1 : {'wetting' :Lw,
+ 'nonwetting': Lnw},#
+ 2 : {'wetting' :Lw,
+ 'nonwetting': Lnw}
}
-l_param = 40
+
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :l_param,
- 'nonwetting': l_param},#
- 2 : {'wetting' :l_param,
- 'nonwetting': l_param}
+ 1 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw},#
+ 2 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw}
}
## relative permeabilty functions on subdomain 1
@@ -177,7 +189,7 @@ def rel_perm1w_prime(s):
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
# # definition of the derivatives of the relative permeabilities
# # relative permeabilty functions on subdomain 1
@@ -187,7 +199,7 @@ def rel_perm1nw_prime(s):
#
# def rel_perm2nw_prime(s):
# # relative permeabilty on subdomain1
-# return 2*(l_param_w1-s)
+# return -2*(l_param_w1-s)
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -364,211 +376,36 @@ p_e_sym = {
# 5: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)}
}
-# pc_e_sym = {
-# 1: -1*p_e_sym[1]['wetting'],
-# 2: -1*p_e_sym[2]['wetting'],
-# # 3: -1*p_e_sym[3]['wetting'],
-# # 4: -1*p_e_sym[4]['wetting'],
-# # 5: -1*p_e_sym[5]['wetting']
-# }
-
-pc_e_sym = {
- 1: p_e_sym[1]['nonwetting'] - p_e_sym[1]['wetting'],
- 2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting'],
- # 3: -1*p_e_sym[3]['wetting'],
- # 4: -1*p_e_sym[4]['wetting'],
- # 5: -1*p_e_sym[5]['wetting']
-}
-
-
-# #### Manufacture source expressions with sympy
-# ###############################################################################
-# ## subdomain1
-# x, y = sym.symbols('x[0], x[1]') # needed by UFL
-# t = sym.symbols('t', positive=True)
-# #f = -sym.diff(u, x, 2) - sym.diff(u, y, 2) # -Laplace(u)
-# #f = sym.simplify(f) # simplify f
-# p1_w = 1 - (1+t**2)*(1 + x**2 + (y-0.5)**2)
-# p1_nw = t*(1-(y-0.5) - x**2)**2 - sym.sqrt(2+t**2)*(1-(y-0.5))
-#
-# #dtS1_w = sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)
-# #dtS1_nw = -sym.diff(S_pc_rel_sym[1](p1_nw - p1_w), t, 1)
-# dtS1_w = porosity[1]*sym.diff(sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ), t, 1)
-# dtS1_nw = -porosity[1]*sym.diff(sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ), t, 1)
-# print("dtS1_w = ", dtS1_w, "\n")
-# print("dtS1_nw = ", dtS1_nw, "\n")
-#
-# #dxdxflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, x, 1), x, 1)
-# #dydyflux1_w = -sym.diff(relative_permeability[1]['wetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_w, y, 1), y, 1)
-# dxdxflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_w, x, 1), x, 1)
-# dydyflux1_w = -1/viscosity[1]['wetting']*sym.diff(relative_permeability[1]['wetting'](sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_w, y, 1), y, 1)
-#
-# rhs1_w = dtS1_w + dxdxflux1_w + dydyflux1_w
-# rhs1_w = sym.printing.ccode(rhs1_w)
-# print("rhs_w = ", rhs1_w, "\n")
-# #rhs_w = sym.expand(rhs_w)
-# #print("rhs_w", rhs_w, "\n")
-# #rhs_w = sym.collect(rhs_w, x)
-# #print("rhs_w", rhs_w, "\n")
-#
-# #dxdxflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, x, 1), x, 1)
-# #dydyflux1_nw = -sym.diff(relative_permeability[1]['nonwetting'](S_pc_rel_sym[1](p1_nw - p1_w))*sym.diff(p1_nw, y, 1), y, 1)
-# dxdxflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_nw, x, 1), x, 1)
-# dydyflux1_nw = -1/viscosity[1]['nonwetting']*sym.diff(relative_permeability[1]['nonwetting'](1-sym.Piecewise((S_pc_rel[1](p1_nw - p1_w), (p1_nw - p1_w) > 0), (1, True) ))*sym.diff(p1_nw, y, 1), y, 1)
-#
-# rhs1_nw = dtS1_nw + dxdxflux1_nw + dydyflux1_nw
-# rhs1_nw = sym.printing.ccode(rhs1_nw)
-# print("rhs_nw = ", rhs1_nw, "\n")
-#
-# ## subdomain2
-# p2_w = 1 - (1+t**2)*(1 + x**2)
-# p2_nw = t*(1- x**2)**2 - sym.sqrt(2+t**2)*(1-(y-0.5))
-#
-# #dtS2_w = sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)
-# #dtS2_nw = -sym.diff(S_pc_rel_sym[2](p2_nw - p2_w), t, 1)
-# dtS2_w = porosity[2]*sym.diff(sym.Piecewise((sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ), (p2_nw - p2_w) > 0), (1, True) ), t, 1)
-# dtS2_nw = -porosity[2]*sym.diff(sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ), t, 1)
-# print("dtS2_w = ", dtS2_w, "\n")
-# print("dtS2_nw = ", dtS2_nw, "\n")
-#
-# #dxdxflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, x, 1), x, 1)
-# #dydyflux2_w = -sym.diff(relative_permeability[2]['wetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_w, y, 1), y, 1)
-# dxdxflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_w, x, 1), x, 1)
-# dydyflux2_w = -1/viscosity[2]['wetting']*sym.diff(relative_permeability[2]['wetting'](sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_w, y, 1), y, 1)
-#
-# rhs2_w = dtS2_w + dxdxflux2_w + dydyflux2_w
-# rhs2_w = sym.printing.ccode(rhs2_w)
-# print("rhs2_w = ", rhs2_w, "\n")
-# #rhs_w = sym.expand(rhs_w)
-# #print("rhs_w", rhs_w, "\n")
-# #rhs_w = sym.collect(rhs_w, x)
-# #print("rhs_w", rhs_w, "\n")
-#
-# #dxdxflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, x, 1), x, 1)
-# #dydyflux2_nw = -sym.diff(relative_permeability[2]['nonwetting'](S_pc_rel_sym[2](p2_nw - p2_w))*sym.diff(p2_nw, y, 1), y, 1)
-# dxdxflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_nw, x, 1), x, 1)
-# dydyflux2_nw = -1/viscosity[2]['nonwetting']*sym.diff(relative_permeability[2]['nonwetting'](1-sym.Piecewise((S_pc_rel[2](p2_nw - p2_w), (p2_nw - p2_w) > 0), (1, True) ))*sym.diff(p2_nw, y, 1), y, 1)
-#
-# rhs2_nw = dtS2_nw + dxdxflux2_nw + dydyflux2_nw
-# rhs2_nw = sym.printing.ccode(rhs2_nw)
-# print("rhs2_nw = ", rhs2_nw, "\n")
-#
-#
-# ###############################################################################
-#
-# source_expression = {
-# 1: {'wetting': rhs1_w,
-# 'nonwetting': rhs1_nw},
-# 2: {'wetting': rhs2_w,
-# 'nonwetting': rhs2_nw}
-# }
-#
-# p1_w_00 = p1_w.subs(t, 0)
-# p1_nw_00 = p1_nw.subs(t, 0)
-# p2_w_00 = p2_w.subs(t, 0)
-# p2_nw_00 = p2_nw.subs(t, 0)
-# # p1_w_00 = sym.printing.ccode(p1_w_00)
-#
-# initial_condition = {
-# 1: {'wetting': sym.printing.ccode(p1_w_00),
-# 'nonwetting': sym.printing.ccode(p1_nw_00)},#
-# 2: {'wetting': sym.printing.ccode(p2_w_00),
-# 'nonwetting': sym.printing.ccode(p2_nw_00)}
-# }
-#
-# exact_solution = {
-# 1: {'wetting': sym.printing.ccode(p1_w),
-# 'nonwetting': sym.printing.ccode(p1_nw)},#
-# 2: {'wetting': sym.printing.ccode(p2_w),
-# 'nonwetting': sym.printing.ccode(p2_nw)}
-# }
-#
-# # similary to the outer boundary dictionary, if a patch has no outer boundary
-# # None should be written instead of an expression. This is a bit of a brainfuck:
-# # dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
-# # Since a domain patch can have several disjoint outer boundary parts, the expressions
-# # need to get an enumaration index which starts at 0. So dirichletBC[ind][j] is
-# # the dictionary of outer dirichlet conditions of subdomain ind and boundary part j.
-# # finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting'] return
-# # the actual expression needed for the dirichlet condition for both phases if present.
-# dirichletBC = {
-# #subdomain index: {outer boudary part index: {phase: expression}}
-# 1: { 0: {'wetting': sym.printing.ccode(p1_w),
-# 'nonwetting': sym.printing.ccode(p1_nw)}},
-# 2: { 0: {'wetting': sym.printing.ccode(p2_w),
-# 'nonwetting': sym.printing.ccode(p2_nw)}}
-# }
-
-# turn above symbolic code into exact solution for dolphin and
-# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
+pc_e_sym = dict()
for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
if isR:
- subdomain_has_phases = ["wetting"]
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'] - p_e_sym[subdomain]['wetting']})
+
+symbols = {"x": x,
+ "y": y,
+ "t": t}
+# turn above symbolic code into exact solution for dolphin and
+# construct the rhs that matches the above exact solution.
+exact_solution_example = hlp.generate_exact_solution_expressions(
+ symbols=symbols,
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym,
+ symbolic_capillary_pressure=pc_e_sym,
+ saturation_pressure_relationship=S_pc_sym,
+ saturation_pressure_relationship_prime=S_pc_sym_prime,
+ viscosity=viscosity,
+ porosity=porosity,
+ relative_permeability=relative_permeability,
+ relative_permeability_prime=ka_prime,
+ densities=densities,
+ gravity_acceleration=gravity_acceleration,
+ include_gravity=include_gravity,
+ )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
@@ -612,7 +449,7 @@ write_to_file = {
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol = 1E-6, debug = False)
+simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol = solver_tol, debug = True)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
@@ -635,7 +472,7 @@ simulation.set_parameters(output_dir = "./output/",#
dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
densities=densities,
- include_gravity=True,
+ include_gravity=include_gravity,
write2file = write_to_file,#
)
diff --git a/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py b/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py
index 1b31d37..a2f841a 100755
--- a/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py
+++ b/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py
@@ -7,11 +7,30 @@ import typing as tp
import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
+import helpers as hlp
#import ufl as ufl
# init sympy session
sym.init_printing()
+solver_tol = 1E-6
+
+############ GRID #######################ü
+mesh_resolution = 31
+timestep_size = 0.001
+number_of_timesteps = 15
+# 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
+starttime = 0
+
+include_gravity = True
+Lw = 10/timestep_size
+Lnw = Lw
+
+l_param_w = 50
+l_param_nw = l_param_w
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0,-1.0), #
@@ -88,17 +107,6 @@ isRichards = {
}
-solver_tol = 1E-8
-
-############ GRID #######################ü
-mesh_resolution = 31
-timestep_size = 0.0001
-number_of_timesteps = 1500
-# 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
-starttime = 0
-
viscosity = {#
# subdom_num : viscosity
1 : {'wetting' :1,
@@ -123,8 +131,7 @@ densities = {
gravity_acceleration = 9.81
-Lw = 10/timestep_size
-Lnw = Lw
+
L = {#
# subdom_num : subdomain L for L-scheme
1 : {'wetting' :Lw,
@@ -133,13 +140,13 @@ L = {#
'nonwetting': Lnw}
}
-l_param = 10
+
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :l_param,
- 'nonwetting': l_param},#
- 2 : {'wetting' :l_param,
- 'nonwetting': l_param}
+ 1 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw},#
+ 2 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw}
}
## relative permeabilty functions on subdomain 1
@@ -189,7 +196,7 @@ def rel_perm1w_prime(s):
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
# # definition of the derivatives of the relative permeabilities
# # relative permeabilty functions on subdomain 1
@@ -199,7 +206,7 @@ def rel_perm2w_prime(s):
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -331,6 +338,70 @@ sat_pressure_relationship = {
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
+symbols = { "x": x,
+ "y": y,
+ "t": t}
+
+epsilon_x_inner = 0.7
+epsilon_x_outer = 0.99
+epsilon_y_inner = epsilon_x_inner
+epsilon_y_outer = epsilon_x_outer
+
+def mollifier(x, epsilon):
+ """ one d mollifier """
+ out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1)
+ return out_expr
+
+mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner)
+
+pw_sym_x = sym.Piecewise(
+ (mollifier_handle(x), x**2 < epsilon_x_outer**2),
+ (0, True)
+)
+pw_sym_y = sym.Piecewise(
+ (mollifier_handle(y), y**2 < epsilon_y_outer**2),
+ (0, True)
+)
+
+def mollifier2d(x, y, epsilon):
+ """ one d mollifier """
+ out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1)
+ return out_expr
+
+mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer)
+
+pw_sym2d_x = sym.Piecewise(
+ (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2),
+ (0, True)
+)
+
+zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise(
+ (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))),
+ (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise(
+ (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))),
+ (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise(
+ (1, y<=-2*epsilon_x_inner),
+ (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))),
+ (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y
+gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x
+cutoff = gaussian/(gaussian + zero_on_shrinking)
+
+
sat_sym = {
1: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t),
2: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t)
@@ -341,13 +412,13 @@ Spc = {
2: sym.Piecewise((pc_saturation_sym[2](sat_sym[2]), sat_sym[2] > 0), (pc_saturation_sym[2](sat_sym[2]), 2>=sat_sym[2]), (0, True))
}
-p1w = 1 - (1+t*t)*(1 + x*x + y*y)
+p1w = (-1 - (1+t*t)*(1 + x*x + y*y))*cutoff
p2w = p1w
p_e_sym = {
1: {'wetting': p1w,
- 'nonwetting': p1w + Spc[1]},
- 2: {'wetting': 1 - (1+t*t)*(1 + x*x + y*y),
- 'nonwetting': p2w + Spc[2]},
+ 'nonwetting': (p1w + Spc[1])*cutoff},
+ 2: {'wetting': p2w,
+ 'nonwetting': (p2w + Spc[2])*cutoff},
}
pc_e_sym = {
@@ -361,76 +432,24 @@ pc_e_sym = {
# 2: -1*p_e_sym[2]['wetting'],
# }
-# turn above symbolic code into exact solution for dolphin and
-# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
-for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
- if isR:
- subdomain_has_phases = ["wetting"]
- else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+exact_solution_example = hlp.generate_exact_solution_expressions(
+ symbols=symbols,
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym,
+ symbolic_capillary_pressure=pc_e_sym,
+ saturation_pressure_relationship=S_pc_sym,
+ saturation_pressure_relationship_prime=S_pc_sym_prime,
+ viscosity=viscosity,#
+ porosity=porosity,
+ relative_permeability=relative_permeability,#
+ relative_permeability_prime=ka_prime,
+ densities=densities,#
+ gravity_acceleration=gravity_acceleration,
+ include_gravity=include_gravity,
+ )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
@@ -474,8 +493,8 @@ write_to_file = {
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug = False)
-simulation.set_parameters(output_dir = "./output/",#
+simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug = True)
+simulation.set_parameters(output_dir = "./output/with_dirichlet_zero",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
interface_def_points = interface_def_points,#
@@ -497,7 +516,7 @@ simulation.set_parameters(output_dir = "./output/",#
dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
densities=densities,
- include_gravity=True,
+ include_gravity=include_gravity,
write2file = write_to_file,#
)
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 d666eac..a12790d 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
@@ -23,12 +23,12 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 14
+mesh_resolution = 40
# ----------------------------------------:-----------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.0001
-number_of_timesteps = 10
+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 = 10
diff --git a/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py b/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py
index 74af206..3510d8a 100755
--- a/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py
+++ b/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py
@@ -7,6 +7,7 @@ import sympy as sym
# import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
+import helpers as hlp
# import ufl as ufl
# init sympy session
@@ -15,22 +16,25 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 51
+mesh_resolution = 50
# ----------------------------------------:-------------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-timestep_size = 0.005
-number_of_timesteps = 160
+timestep_size = 0.001
+number_of_timesteps = 1500
# 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
+number_of_timesteps_to_analyse = 10
starttime = 0
-Lw = 1000
+Lw = 1/timestep_size
Lnw = Lw
-l_param_w = 80
-l_param_nw = 80
+l_param_w = 40
+l_param_nw = 40
+
+solver_tol = 5e-8
+include_gravity = True
# ----------------------------------------------------------------------------#
# ------------------- Domain and Interface -----------------------------------#
@@ -258,7 +262,7 @@ def rel_perm1w_prime(s):
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
@@ -268,7 +272,7 @@ def rel_perm2w_prime(s):
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return 3*(1-s)**2
+ return -3*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -354,13 +358,13 @@ x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym = {
- 1: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 1: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
- 2: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x),
+ 2: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
- 3: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x),
+ 3: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
- 4: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 4: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t}
}
@@ -376,78 +380,33 @@ for subdomain, isR in isRichards.items():
if isR:
pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - p_e_sym[subdomain]['wetting'].copy()})
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
+symbols = {"x": x,
+ "y": y,
+ "t": t}
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
-for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
- if isR:
- subdomain_has_phases = ["wetting"]
- else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+exact_solution_example = hlp.generate_exact_solution_expressions(
+ symbols=symbols,
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym,
+ symbolic_capillary_pressure=pc_e_sym,
+ saturation_pressure_relationship=S_pc_sym,
+ saturation_pressure_relationship_prime=S_pc_sym_prime,
+ viscosity=viscosity,
+ porosity=porosity,
+ relative_permeability=relative_permeability,
+ relative_permeability_prime=ka_prime,
+ densities=densities,
+ gravity_acceleration=gravity_acceleration,
+ include_gravity=include_gravity,
+ )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
+
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
@@ -484,9 +443,11 @@ write_to_file = {
'L_iterations': True
}
+output_string = "./output/like_RR_number_of_timesteps{}_".format(number_of_timesteps)
+
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=True, LDDsolver_tol=1E-6)
-simulation.set_parameters(output_dir="./output/",
+simulation = ldd.LDDsimulation(tol=1E-14, LDDsolver_tol=solver_tol, debug=False)
+simulation.set_parameters(output_dir=output_string,
subdomain_def_points=subdomain_def_points,
isRichards=isRichards,
interface_def_points=interface_def_points,
@@ -508,7 +469,7 @@ simulation.set_parameters(output_dir="./output/",
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
- include_gravity=True,
+ include_gravity=include_gravity,
write2file=write_to_file,
)
diff --git a/TP-one-patch/TP-one-patch.py b/TP-one-patch/TP-one-patch.py
index c46a8e6..ad98e9a 100755
--- a/TP-one-patch/TP-one-patch.py
+++ b/TP-one-patch/TP-one-patch.py
@@ -7,11 +7,32 @@ import typing as tp
import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
+import helpers as hlp
#import ufl as ufl
# init sympy session
sym.init_printing()
+
+solver_tol = 5E-6
+
+############ GRID #######################ü
+mesh_resolution = 30
+timestep_size = 0.001
+number_of_timesteps = 600
+# 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
+starttime = 0
+
+Lw = 1/timestep_size
+Lnw=Lw
+
+l_param_w = 80
+l_param_nw = 80
+
+include_gravity = True
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0,-1.0), #
@@ -20,7 +41,11 @@ sub_domain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(-1.0,1.0)]
subdomain0_outer_boundary_verts = {
- 0: sub_domain0_vertices
+ 0: [sub_domain0_vertices[0],
+ sub_domain0_vertices[1],
+ sub_domain0_vertices[2],
+ sub_domain0_vertices[3],
+ sub_domain0_vertices[0]]
}
# list of subdomains given by the boundary polygon vertices.
@@ -48,21 +73,6 @@ isRichards = {
0: False, #
}
-
-solver_tol = 1E-6
-
-############ GRID #######################ü
-mesh_resolution = 30
-timestep_size = 0.0005
-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 = 11
-starttime = 0
-
-Lw = 0.01/timestep_size
-Lnw=Lw
-
viscosity = {#
# subdom_num : viscosity
0 : {'wetting' :1,
@@ -88,8 +98,6 @@ L = {#
'nonwetting': Lnw},#
}
-l_param_w = 100
-l_param_nw = 100
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
0: {'wetting' :l_param_w,
@@ -126,7 +134,7 @@ def rel_perm1w_prime(s):
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -147,7 +155,6 @@ def saturation(pc, index):
# inverse capillary pressure-saturation-relationship
return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
def saturation_sym(pc, index):
# inverse capillary pressure-saturation-relationship
return 1/((1 + pc)**(1/(index + 1)))
@@ -160,6 +167,23 @@ def saturation_sym_prime(pc, index):
return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
+# def saturation(pc, index):
+# # inverse capillary pressure-saturation-relationship
+# return df.conditional(pc > 0, -index*pc, 1)
+#
+#
+# def saturation_sym(pc, index):
+# # inverse capillary pressure-saturation-relationship
+# return -index*pc
+#
+#
+# # derivative of S-pc relationship with respect to pc. This is needed for the
+# # construction of a analytic solution.
+# def saturation_sym_prime(pc, index):
+# # inverse capillary pressure-saturation-relationship
+# return -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 = {
@@ -181,91 +205,138 @@ sat_pressure_relationship = {
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
+epsilon_x_inner = 0.7
+epsilon_x_outer = 0.99
+epsilon_y_inner = epsilon_x_inner
+epsilon_y_outer = epsilon_x_outer
+
+def mollifier(x, epsilon):
+ """ one d mollifier """
+ out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1)
+ return out_expr
+
+mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner)
+
+pw_sym_x = sym.Piecewise(
+ (mollifier_handle(x), x**2 < epsilon_x_outer**2),
+ (0, True)
+)
+pw_sym_y = sym.Piecewise(
+ (mollifier_handle(y), y**2 < epsilon_y_outer**2),
+ (0, True)
+)
+
+def mollifier2d(x, y, epsilon):
+ """ one d mollifier """
+ out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1)
+ return out_expr
+
+mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer)
+
+pw_sym2d_x = sym.Piecewise(
+ (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2),
+ (0, True)
+)
+
+zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise(
+ (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))),
+ (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise(
+ (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))),
+ (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise(
+ (1, y<=-2*epsilon_x_inner),
+ (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))),
+ (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))),
+ (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))),
+ (1, True),
+)
+
+zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y
+gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x
+cutoff = gaussian/(gaussian + zero_on_shrinking)
+
+# # construction of differentiable characteristic function.
+# def smooth_characteristic_func_on_epsilon_shrinking_of_subdomain0(x, y, epsilon_x_inner, epsilon_y_inner, epsilon_x_outer, epsilon_y_outer):
+# dist_to_complement_x = ft.partial(mollifier, epsilon=epsilon_x_inner)
+# dist_to_complement_y = ft.partial(mollifier, epsilon=epsilon_y_inner)
+# dist_to_complement = dist_to_complement_y(y)*dist_to_complement_x(x)
+# dist_to_eps_shrinking_x = ft.partial(zero_outside_epsilon_thickening_of_subdomain, epsilon=epsilon_x_outer)
+# dist_to_eps_shrinking_y = ft.partial(zero_outside_epsilon_thickening_of_subdomain, epsilon=epsilon_y_outer)
+# dist_to_eps_shrinking = dist_to_eps_shrinking_y(y)*dist_to_eps_shrinking_x(x)
+# return dist_to_complement/(dist_to_eps_shrinking + dist_to_complement)
+#
+
+# def dist_to_epsilon_thickening_of_subdomain0_complement(x, y, epsilon):
+# """ calculates the (euklidian distance)^2 of a point x,y to the epsilon
+# thickening of the complement of the domain.
+# """
+# is_inside = ((1-sym.Abs(x) > epsilon) & (1-sym.Abs(y) > epsilon))
+# sym.Piecewise((0, is_inside))
+
+# p_e_sym = {
+# 0: {'wetting': (-3 - (1+t*t)*(1 + x*x + y*y))*cutoff,
+# 'nonwetting': (-1 -t*(1+y + x**2)**2)*cutoff},
+# }
+
p_e_sym = {
- 0: {'wetting': -3 - (1+t*t)*(1 + x*x + y*y),
- 'nonwetting': -1 -t*(1+y + x**2)**2*(1-y)**2},
+ 0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2,
+ 'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
}
-pc_e_sym = {
- 0: p_e_sym[0]['nonwetting'] - p_e_sym[0]['wetting'],
-}
+print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
+# # pw_sym_x*pw_sym_y
+# p_e_sym = {
+# 0: {'wetting': -3*pw_sym2d_x + 0*t,
+# 'nonwetting': -1*pw_sym_x*pw_sym_y+ 0*t},
+# }
-# pc_e_sym = {
-# 0: -1*p_e_sym[1]['wetting'],
-# 2: -1*p_e_sym[2]['wetting'],
+# p_e_sym = {
+# 0: {'wetting': -3*cutoff + 0*t,
+# 'nonwetting': -1*zero_on_shrinking+ 0*t},
# }
-# turn above symbolic code into exact solution for dolphin and
-# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
+
+pc_e_sym = dict()
for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
if isR:
- subdomain_has_phases = ["wetting"]
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
+
+symbols = {"x": x,
+ "y": y,
+ "t": t}
+# turn above symbolic code into exact solution for dolphin and
+# construct the rhs that matches the above exact solution.
+exact_solution_example = hlp.generate_exact_solution_expressions(
+ symbols=symbols,
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym,
+ symbolic_capillary_pressure=pc_e_sym,
+ saturation_pressure_relationship=S_pc_sym,
+ saturation_pressure_relationship_prime=S_pc_sym_prime,
+ viscosity=viscosity,
+ porosity=porosity,
+ relative_permeability=relative_permeability,
+ relative_permeability_prime=ka_prime,
+ densities=densities,
+ gravity_acceleration=gravity_acceleration,
+ include_gravity=include_gravity,
+ )
+source_expression = exact_solution_example['source']
+exact_solution = exact_solution_example['exact_solution']
+initial_condition = exact_solution_example['initial_condition']
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
@@ -307,10 +378,11 @@ write_to_file = {
'L_iterations': True
}
+output_string = "./output/number_of_timesteps{}_".format(number_of_timesteps)
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol = solver_tol, debug = True)
-simulation.set_parameters(output_dir = "./output/",#
+simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug=False)
+simulation.set_parameters(output_dir = output_string,#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
interface_def_points = interface_def_points,#
@@ -332,7 +404,7 @@ simulation.set_parameters(output_dir = "./output/",#
dirichletBC_expression_strings = dirichletBC,#
exact_solution = exact_solution,#
densities=densities,
- include_gravity=True,
+ include_gravity=include_gravity,
write2file = write_to_file,#
)
--
GitLab