From ed4e464b40840d4536d9abb1c8be95f2096e9103 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Thu, 27 Jun 2019 18:57:27 +0200
Subject: [PATCH] play_with_parameters
---
LDDsimulation/LDDsimulation.py | 50 +++++++++-----
.../TP-TP-2-patch-pure-dd.py | 66 +++++++++++++------
...ed_soil_with_inner_patch_const_solution.py | 14 ++--
.../TP-TP-layered_soil_with_inner_patch.py | 6 +-
4 files changed, 88 insertions(+), 48 deletions(-)
diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index a076bc7..71b503d 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -103,7 +103,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 = 100
# 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
@@ -148,20 +148,28 @@ class LDDsimulation(object):
# print("\nLinear Algebra Backends:")
# df.list_linear_algebra_backends()
# # print("\nLinear Solver Methods:")
- # df.list_linear_solver_methods()
+ df.list_linear_solver_methods()
+ print("\n")
+ df.list_krylov_solver_methods()
+ print("\n")
# # print("\nPeconditioners for Krylov Solvers:")
- # df.list_krylov_solver_preconditioners()
+ df.list_krylov_solver_preconditioners()
+ # df.LinearSolver_default_parameters()
+
+ df.info(df.parameters, True)
+
+ self.solver_type_is_Kryov = False
### Define the linear solver to be used.
- self.solver = 'bicgstab' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
- self.preconditioner = 'jacobi'#'hypre_amg' #'ilu'#'hypre_amg' # incomplete LU factorization
+ self.solver = 'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
+ self.preconditioner = 'default'#jacobi#'hypre_amg' #'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.
self.solver_parameters = {
'nonzero_initial_guess': True,
- 'absolute_tolerance': 1E-12,
- 'relative_tolerance': 1E-10,
+ 'absolute_tolerance': 1E-14,
+ 'relative_tolerance': 1E-12,
'maximum_iterations': 1000,
'report': False
}
@@ -937,18 +945,24 @@ class LDDsimulation(object):
}
)
# setup the linear solvers and set the solver parameters.
- self.subdomain[subdom_num].linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
- # we use the previous iteration as an initial guess for the linear solver.
- solver_param = self.subdomain[subdom_num].linear_solver.parameters
- solver_param['nonzero_initial_guess'] = self.solver_parameters['nonzero_initial_guess']
- solver_param['absolute_tolerance'] = self.solver_parameters['absolute_tolerance']
- solver_param['relative_tolerance'] = self.solver_parameters['relative_tolerance']
- solver_param['maximum_iterations'] = self.solver_parameters['maximum_iterations']
- solver_param['report'] = self.solver_parameters['report']
+ # df.LinearSolver(self.solver)
+ if self.solver_type_is_Kryov:
+ self.subdomain[subdom_num].linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
+ # we use the previous iteration as an initial guess for the linear solver.
+ solver_param = self.subdomain[subdom_num].linear_solver.parameters
+ df.info(solver_param, True)
+ solver_param['nonzero_initial_guess'] = self.solver_parameters['nonzero_initial_guess']
+ solver_param['absolute_tolerance'] = self.solver_parameters['absolute_tolerance']
+ solver_param['relative_tolerance'] = self.solver_parameters['relative_tolerance']
+ solver_param['maximum_iterations'] = self.solver_parameters['maximum_iterations']
+ solver_param['report'] = self.solver_parameters['report']
+ else:
+ self.subdomain[subdom_num].linear_solver = df.LUSolver()
+
# ## print out set solver parameters
- # for parameter, value in self.linear_solver.parameters.items():
- # print(f"parameter: {parameter} = {value}")
- # df.info(solver_param, True)
+ for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
+ print(f"parameter: {parameter} = {value}")
+
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+f"interface_marker_subdomain{subdom_num}.pvd"
df.File(filepath) << self.subdomain[subdom_num].interface_marker
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 cea1cda..a12fded 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
@@ -88,10 +88,12 @@ isRichards = {
}
+solver_tol = 1E-8
+
############ GRID #######################ΓΌ
-mesh_resolution = 51
-timestep_size = 0.01
-number_of_timesteps = 50
+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
@@ -121,15 +123,17 @@ densities = {
gravity_acceleration = 9.81
+Lw = 10/timestep_size
+Lnw = Lw
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
+l_param = 10
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :l_param,
@@ -157,10 +161,10 @@ subdomain1_rel_perm = {
## relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
- return s**3
+ return s**2
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
- return (1-s)**3
+ return (1-s)**2
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)
@@ -191,11 +195,11 @@ def rel_perm1nw_prime(s):
# # relative permeabilty functions on subdomain 1
def rel_perm2w_prime(s):
# relative permeabilty on subdomain1
- return 3*s**2
+ return 2*s
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return 3*(1-s)**2
+ return 2*(1-s)
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -226,6 +230,16 @@ def saturation(pc, index):
return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
+
+def pc_sat_rel_sym(S, index):
+ # capillary pressure-saturation-relationship
+ return 1/S**(index+1) -1
+
+pc_saturation_sym = {
+ 1: ft.partial(pc_sat_rel_sym, index=1),
+ 2: ft.partial(pc_sat_rel_sym, index=1),
+}
+
def saturation_sym(pc, index):
# inverse capillary pressure-saturation-relationship
return 1/((1 + pc)**(1/(index + 1)))
@@ -242,21 +256,21 @@ 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),
+ 2: ft.partial(saturation_sym, index=1),
# 3: ft.partial(saturation_sym, index=2),
# 4: 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),
+ 2: ft.partial(saturation_sym_prime, index=1),
# 3: ft.partial(saturation_sym_prime, index=2),
# 4: ft.partial(saturation_sym_prime, index=1)
}
sat_pressure_relationship = {
1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
+ 2: ft.partial(saturation, index=1),
# 3: ft.partial(saturation, index=2),
# 4: ft.partial(saturation, index=1)
}
@@ -317,11 +331,23 @@ sat_pressure_relationship = {
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
+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)
+ }
+
+Spc = {
+ 1: sym.Piecewise((pc_saturation_sym[1](sat_sym[1]), sat_sym[1] > 0), (pc_saturation_sym[1](sat_sym[1]), 1>=sat_sym[1]), (0, True)),
+ 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)
+p2w = p1w
p_e_sym = {
- 1: {'wetting': 1 - (1+t*t)*(1 + x*x + y*y),
- 'nonwetting': -t*(1-y - x**2)**2 - sym.sqrt(2+t**2)*(1-y)},
- 2: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x),
- 'nonwetting': -t*(1- x**2)**2 - sym.sqrt(2+t**2)*(1-y)},
+ 1: {'wetting': p1w,
+ 'nonwetting': p1w + Spc[1]},
+ 2: {'wetting': 1 - (1+t*t)*(1 + x*x + y*y),
+ 'nonwetting': p2w + Spc[2]},
}
pc_e_sym = {
@@ -448,7 +474,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 = False)
simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
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 8200e0e..d666eac 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,22 +23,22 @@ sym.init_printing()
# ----------------------------------------------------------------------------#
# ------------------- MESH ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-mesh_resolution = 51
+mesh_resolution = 14
# ----------------------------------------:-----------------------------------#
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
-timestep_size = 0.01
+timestep_size = 0.0001
number_of_timesteps = 10
# 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 = 0.25
-Lnw = 0.25
+Lw = 100/timestep_size
+Lnw = Lw
-l_param_w = 200
-l_param_nw = 200
+l_param_w = 40
+l_param_nw = 40
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
@@ -778,7 +778,7 @@ write_to_file = {
}
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=True, LDDsolver_tol=5E-5)
+simulation = ldd.LDDsimulation(tol=1E-14, debug=True, LDDsolver_tol=1E-7)
simulation.set_parameters(output_dir="./output/",
subdomain_def_points=subdomain_def_points,
isRichards=isRichards,
diff --git a/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py b/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
index abc1414..868fd4f 100755
--- a/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
+++ b/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
@@ -28,14 +28,14 @@ mesh_resolution = 50
# ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------#
timestep_size = 0.0001
-number_of_timesteps = 10
+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 = 4
starttime = 0
-Lw = 0.5
-Lnw = 0.4
+Lw = 10000
+Lnw = 10000
l_param_w = 30
l_param_nw = 40
--
GitLab