diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py
index 5b066bfe166b15aab1681e67efaa7af4d1a6175f..0737454e106a9e38a8c24c6fe193871ea834db79 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py
@@ -13,7 +13,13 @@ import datetime
import os
import pandas as pd
-# check if output directory exists
+
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
@@ -21,57 +27,57 @@ else:
print("Directory ", './output', " already exists. Will use as output \
directory")
-
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
-# init sympy session
-sym.init_printing()
-# solver_tol = 6E-7
-use_case = "TP-R-layered-soil-with-inner-patch-all-params-one"
-# name of this very file. Needed for log output.
-thisfile = "TP-R-layered_soil_with_inner_patch-all-params-one.py"
+# Name of the usecase that will be printed during simulation.
+use_case = "TP-R-layered-soil-with-inner-patch-realistic"
+# The name of this very file. Needed for creating log output.
+thisfile = "TP-R-layered_soil_with_inner_patch-realistic.py"
-max_iter_num = 150
+# GENERAL SOLVER CONFIG ######################################################
+# maximal iteration per timestep
+max_iter_num = 700
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False
resolutions = {
# 1: 2e-6, # h=2
# 2: 2e-6, # h=1.1180
# 4: 2e-6, # h=0.5590
# 8: 2e-6, # h=0.2814
- # 16: 7e-6, # h=0.1412
- # 32: 5e-6,
- 64: 2e-6,
+ 16: 8e-6, # h=0.1412
+ # 32: 2e-6,
+ # 64: 2e-6,
# 128: 2e-6
}
-# GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.00025
-number_of_timesteps = 4000
-plot_timestep_every = 4
-# 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 = 5
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
starttimes = [0.0]
+timestep_size = 0.01
+number_of_timesteps = 5
-Lw1 = 0.25 # /timestep_size
+# LDD scheme parameters ######################################################
+
+Lw1 = 0.025 # /timestep_size
Lnw1 = Lw1
-Lw2 = 0.25 # /timestep_size
+Lw2 = 0.025 # /timestep_size
Lnw2 = Lw2
-Lw3 = 0.25 # /timestep_size
+Lw3 = 0.025 # /timestep_size
Lnw3 = Lw3
-Lw4 = 0.25 # /timestep_size
+Lw4 = 0.025 # /timestep_size
Lnw4 = Lw4
-Lw5 = 0.25 # /timestep_size
+Lw5 = 0.025 # /timestep_size
Lnw5 = Lw5
-Lw6 = 0.25 # /timestep_size
+Lw6 = 0.025 # /timestep_size
Lnw6 = Lw6
lambda12_w = 4
@@ -101,45 +107,62 @@ lambda46_nw = 4
lambda56_w = 4
lambda56_nw = 4
-include_gravity = False
-debugflag = False
+include_gravity = True
+debugflag = True
analyse_condition = False
-# if mesh_study:
-output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-# else:
-# # for tol in resolutions.values():
-# # solver_tol = tol
-# output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
-# datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
-# )
-
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
'space_errornorms': True,
+ # save the mesh and marker functions to disk
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
+ 'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True
}
else:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ 'L_iterations_per_timestep': False,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'subsequent_errors': True
}
+# OUTPUT FILE STRING #########################################################
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+
+# DOMAIN AND INTERFACES #######################################################
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(1.0,-1.0),#
@@ -189,14 +212,18 @@ interface56_vertices = [interface46_vertices[1],
interface34_vertices = [interface36_vertices[1],
interface23_vertices[2]]
-# interface36
-interface45_vertices = [interface56_vertices[0],
- df.Point(0.7, -0.2),
+# Interface 45 needs to be split, because of the shape. There can be triangles
+# with two facets on the interface and this creates a rogue dof type error when
+# integrating over that particular interface. Accordingly, the lambda_param
+# dictionary has two entries for that interface.
+interface45_vertices_a = [interface56_vertices[0],
+ df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
+ ]
+interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
interface25_vertices[0]
]
-
# interface_vertices introduces a global numbering of interfaces.
interface_def_points = [interface12_vertices,
interface23_vertices,
@@ -204,11 +231,11 @@ interface_def_points = [interface12_vertices,
interface25_vertices,
interface34_vertices,
interface36_vertices,
- interface45_vertices,
+ interface45_vertices_a,
+ interface45_vertices_b,
interface46_vertices,
interface56_vertices,
]
-
adjacent_subdomains = [[1,2],
[2,3],
[2,4],
@@ -216,6 +243,7 @@ adjacent_subdomains = [[1,2],
[3,4],
[3,6],
[4,5],
+ [4,5],
[4,6],
[5,6]
]
@@ -225,18 +253,17 @@ subdomain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
interface12_vertices[2],
interface12_vertices[3],
- interface12_vertices[4], # southern boundary, 12 interface
- subdomain0_vertices[2], # eastern boundary, outer boundary
- subdomain0_vertices[3]] # northern boundary, outer on_boundary
+ interface12_vertices[4], # southern boundary, 12 interface
+ subdomain0_vertices[2], # eastern boundary, outer boundary
+ subdomain0_vertices[3]] # northern boundary, outer on_boundary
# vertex coordinates of the outer boundaries. If it can not be specified as a
-# polygon, use an entry per boundary polygon.
-# This information is used for defining
+# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
- 0: [subdomain1_vertices[4],
- subdomain1_vertices[5], # eastern boundary, outer boundary
+ 0: [subdomain1_vertices[4], #
+ subdomain1_vertices[5], # eastern boundary, outer boundary
subdomain1_vertices[6],
subdomain1_vertices[0]]
}
@@ -246,12 +273,12 @@ subdomain2_vertices = [interface23_vertices[0],
interface23_vertices[1],
interface23_vertices[2],
interface24_vertices[1],
- interface25_vertices[1], # southern boundary, 23 interface
- subdomain1_vertices[4], # eastern boundary, outer boundary
+ interface25_vertices[1], # southern boundary, 23 interface
+ subdomain1_vertices[4], # eastern boundary, outer boundary
subdomain1_vertices[3],
subdomain1_vertices[2],
subdomain1_vertices[1],
- subdomain1_vertices[0] ] # northern boundary, 12 interface
+ subdomain1_vertices[0] ] # northern boundary, 12 interface
subdomain2_outer_boundary_verts = {
0: [subdomain2_vertices[9],
@@ -279,7 +306,8 @@ subdomain3_outer_boundary_verts = {
# subdomain3
subdomain4_vertices = [interface46_vertices[0],
interface46_vertices[1],
- interface45_vertices[1],
+ # interface45_vertices[1],
+ interface45_vertices_a[1],
interface24_vertices[1],
interface24_vertices[0],
interface34_vertices[1]
@@ -292,10 +320,11 @@ subdomain5_vertices = [interface56_vertices[0],
interface56_vertices[2],
interface25_vertices[1],
interface25_vertices[0],
- interface45_vertices[1],
- interface45_vertices[0]
+ interface45_vertices_b[1],
+ interface45_vertices_b[0]
]
+
subdomain5_outer_boundary_verts = {
0: [subdomain5_vertices[2],
subdomain5_vertices[3]]
@@ -343,7 +372,7 @@ outer_boundary_def_points = {
6: subdomain6_outer_boundary_verts
}
-
+# MODEL CONFIGURATION #########################################################
isRichards = {
1: True,
2: True,
@@ -435,32 +464,36 @@ L = {
# interface25_vertices,
# interface34_vertices,
# interface36_vertices,
-# interface45_vertices,
+# interface45_vertices_a,
+# interface45_vertices_b,
# interface46_vertices,
# interface56_vertices,
# ]
lambda_param = {
0: {'wetting': lambda12_w,
- 'nonwetting': lambda12_nw},#
+ 'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
- 'nonwetting': lambda23_nw},#
+ 'nonwetting': lambda23_nw},
2: {'wetting': lambda24_w,
- 'nonwetting': lambda24_nw},#
+ 'nonwetting': lambda24_nw},
3: {'wetting': lambda25_w,
- 'nonwetting': lambda25_nw},#
+ 'nonwetting': lambda25_nw},
4: {'wetting': lambda34_w,
- 'nonwetting': lambda34_nw},#
+ 'nonwetting': lambda34_nw},
5: {'wetting': lambda36_w,
- 'nonwetting': lambda36_nw},#
+ 'nonwetting': lambda36_nw},
6: {'wetting': lambda45_w,
- 'nonwetting': lambda45_nw},#
- 7: {'wetting': lambda46_w,
- 'nonwetting': lambda46_nw},#
- 8: {'wetting': lambda56_w,
- 'nonwetting': lambda56_nw},#
+ 'nonwetting': lambda45_nw},
+ 7: {'wetting': lambda45_w,
+ 'nonwetting': lambda45_nw},
+ 8: {'wetting': lambda46_w,
+ 'nonwetting': lambda46_nw},
+ 9: {'wetting': lambda56_w,
+ 'nonwetting': lambda56_nw},
}
+# after Lewis, see pdf file
# after Lewis, see pdf file
intrinsic_permeability = {
1: 1, # sand
@@ -824,37 +857,36 @@ sat_pressure_relationship = {
6: ft.partial(saturation, n_index=2)
}
-#############################################
+###############################################################################
# Manufacture source expressions with sympy #
-#############################################
+###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym = {
- 1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
- 'nonwetting': 0.0*t },
- 2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
- 'nonwetting': 0.0*t },
+ 1: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 'nonwetting': 0*t },
+ 2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 'nonwetting': 0*t },
3: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
4: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
5: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
6: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
}
+
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
- pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- pc_e_sym.update(
- {subdomain: p_e_sym[subdomain]['nonwetting']
- - p_e_sym[subdomain]['wetting']}
- )
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
symbols = {"x": x,
@@ -881,6 +913,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
@@ -897,8 +930,7 @@ dirichletBC = dict()
# subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys():
- # if subdomain has no outer boundary, outer_boundary_def_points[subdomain]
- # is None
+ # subdomain can have no outer boundary
if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None})
else:
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
index d81858530db4e5bedfa155611eedc707a5585338..8d330bf9bb1cd9fab6afcb08ecacce858304a99c 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -13,7 +13,13 @@ import datetime
import os
import pandas as pd
-# check if output directory exists
+
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
@@ -24,15 +30,17 @@ else:
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
-
-# init sympy session
-sym.init_printing()
-# solver_tol = 6E-7
+# Name of the usecase that will be printed during simulation.
use_case = "TP-R-layered-soil-with-inner-patch-realistic"
-# name of this very file. Needed for log output.
+# The name of this very file. Needed for creating log output.
thisfile = "TP-R-layered_soil_with_inner_patch-realistic.py"
+
+# GENERAL SOLVER CONFIG ######################################################
+# maximal iteration per timestep
max_iter_num = 700
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False
resolutions = {
# 1: 2e-6, # h=2
@@ -45,15 +53,14 @@ resolutions = {
# 128: 2e-6
}
-# GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.001
-number_of_timesteps = 5
-plot_timestep_every = 1
-# 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 = 5
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
starttimes = [0.0]
+timestep_size = 0.01
+number_of_timesteps = 5
+
+# LDD scheme parameters ######################################################
Lw1 = 0.025 # /timestep_size
Lnw1 = Lw1
@@ -73,58 +80,69 @@ Lnw5 = Lw5
Lw6 = 0.025 # /timestep_size
Lnw6 = Lw6
-lambda12_w = 40
-lambda12_nw = 40
+lambda12_w = 4
+lambda12_nw = 4
-lambda23_w = 40
-lambda23_nw = 40
+lambda23_w = 4
+lambda23_nw = 4
-lambda24_w = 40
-lambda24_nw= 40
+lambda24_w = 4
+lambda24_nw= 4
-lambda25_w= 40
-lambda25_nw= 40
+lambda25_w= 4
+lambda25_nw= 4
-lambda34_w = 40
-lambda34_nw = 40
+lambda34_w = 4
+lambda34_nw = 4
-lambda36_w = 40
-lambda36_nw = 40
+lambda36_w = 4
+lambda36_nw = 4
-lambda45_w = 40
-lambda45_nw = 40
+lambda45_w = 4
+lambda45_nw = 4
-lambda46_w = 40
-lambda46_nw = 40
+lambda46_w = 4
+lambda46_nw = 4
-lambda56_w = 40
-lambda56_nw = 40
+lambda56_w = 4
+lambda56_nw = 4
-include_gravity = False
+include_gravity = True
debugflag = True
-analyse_condition = True
-
-if mesh_study:
- output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-else:
- for tol in resolutions.values():
- solver_tol = tol
- output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
- )
-
+analyse_condition = False
+
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
'space_errornorms': True,
+ # save the mesh and marker functions to disk
'meshes_and_markers': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True
}
else:
@@ -138,7 +156,13 @@ else:
'subsequent_errors': True
}
+# OUTPUT FILE STRING #########################################################
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+
+# DOMAIN AND INTERFACES #######################################################
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(1.0,-1.0),#
@@ -188,14 +212,18 @@ interface56_vertices = [interface46_vertices[1],
interface34_vertices = [interface36_vertices[1],
interface23_vertices[2]]
-# interface36
-interface45_vertices = [interface56_vertices[0],
- df.Point(0.7, -0.2),
+# Interface 45 needs to be split, because of the shape. There can be triangles
+# with two facets on the interface and this creates a rogue dof type error when
+# integrating over that particular interface. Accordingly, the lambda_param
+# dictionary has two entries for that interface.
+interface45_vertices_a = [interface56_vertices[0],
+ df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
+ ]
+interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
interface25_vertices[0]
]
-
# interface_vertices introduces a global numbering of interfaces.
interface_def_points = [interface12_vertices,
interface23_vertices,
@@ -203,11 +231,11 @@ interface_def_points = [interface12_vertices,
interface25_vertices,
interface34_vertices,
interface36_vertices,
- interface45_vertices,
+ interface45_vertices_a,
+ interface45_vertices_b,
interface46_vertices,
interface56_vertices,
]
-
adjacent_subdomains = [[1,2],
[2,3],
[2,4],
@@ -215,6 +243,7 @@ adjacent_subdomains = [[1,2],
[3,4],
[3,6],
[4,5],
+ [4,5],
[4,6],
[5,6]
]
@@ -224,18 +253,17 @@ subdomain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
interface12_vertices[2],
interface12_vertices[3],
- interface12_vertices[4], # southern boundary, 12 interface
- subdomain0_vertices[2], # eastern boundary, outer boundary
- subdomain0_vertices[3]] # northern boundary, outer on_boundary
+ interface12_vertices[4], # southern boundary, 12 interface
+ subdomain0_vertices[2], # eastern boundary, outer boundary
+ subdomain0_vertices[3]] # northern boundary, outer on_boundary
# vertex coordinates of the outer boundaries. If it can not be specified as a
-# polygon, use an entry per boundary polygon.
-# This information is used for defining
+# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
- 0: [subdomain1_vertices[4],
- subdomain1_vertices[5], # eastern boundary, outer boundary
+ 0: [subdomain1_vertices[4], #
+ subdomain1_vertices[5], # eastern boundary, outer boundary
subdomain1_vertices[6],
subdomain1_vertices[0]]
}
@@ -245,12 +273,12 @@ subdomain2_vertices = [interface23_vertices[0],
interface23_vertices[1],
interface23_vertices[2],
interface24_vertices[1],
- interface25_vertices[1], # southern boundary, 23 interface
- subdomain1_vertices[4], # eastern boundary, outer boundary
+ interface25_vertices[1], # southern boundary, 23 interface
+ subdomain1_vertices[4], # eastern boundary, outer boundary
subdomain1_vertices[3],
subdomain1_vertices[2],
subdomain1_vertices[1],
- subdomain1_vertices[0] ] # northern boundary, 12 interface
+ subdomain1_vertices[0] ] # northern boundary, 12 interface
subdomain2_outer_boundary_verts = {
0: [subdomain2_vertices[9],
@@ -278,7 +306,8 @@ subdomain3_outer_boundary_verts = {
# subdomain3
subdomain4_vertices = [interface46_vertices[0],
interface46_vertices[1],
- interface45_vertices[1],
+ # interface45_vertices[1],
+ interface45_vertices_a[1],
interface24_vertices[1],
interface24_vertices[0],
interface34_vertices[1]
@@ -291,10 +320,11 @@ subdomain5_vertices = [interface56_vertices[0],
interface56_vertices[2],
interface25_vertices[1],
interface25_vertices[0],
- interface45_vertices[1],
- interface45_vertices[0]
+ interface45_vertices_b[1],
+ interface45_vertices_b[0]
]
+
subdomain5_outer_boundary_verts = {
0: [subdomain5_vertices[2],
subdomain5_vertices[3]]
@@ -342,7 +372,7 @@ outer_boundary_def_points = {
6: subdomain6_outer_boundary_verts
}
-
+# MODEL CONFIGURATION #########################################################
isRichards = {
1: True,
2: True,
@@ -379,18 +409,18 @@ viscosity = {
# Dict of the form: { subdom_num : density }
densities = {
- 1: {'wetting': 1, #997
- 'nonwetting': 1}, #1}, #1.225},
- 2: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225},
- 3: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225},
- 4: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225}
- 5: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225},
- 6: {'wetting': 1, #997
- 'nonwetting': 1} #1.225}
+ 1: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1}, #1.225},
+ 2: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225},
+ 3: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225},
+ 4: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225}
+ 5: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225},
+ 6: {'wetting': 997.0, #997
+ 'nonwetting': 1.225} #1.225}
}
gravity_acceleration = 9.81
@@ -398,12 +428,12 @@ gravity_acceleration = 9.81
# https://www.geotechdata.info/parameter/soil-porosity.html
# Dict of the form: { subdom_num : porosity }
porosity = {
- 1: 1, #0.2, # Clayey gravels, clayey sandy gravels
- 2: 1, #0.22, # Silty gravels, silty sandy gravels
- 3: 1, #0.37, # Clayey sands
- 4: 1, #0.2 # Silty or sandy clay
- 5: 1, #
- 6: 1, #
+ 1: 0.2, #0.2, # Clayey gravels, clayey sandy gravels
+ 2: 0.2, #0.22, # Silty gravels, silty sandy gravels
+ 3: 0.2, #0.37, # Clayey sands
+ 4: 0.2, #0.2 # Silty or sandy clay
+ 5: 0.2, #
+ 6: 0.2, #
}
# subdom_num : subdomain L for L-scheme
@@ -434,40 +464,43 @@ L = {
# interface25_vertices,
# interface34_vertices,
# interface36_vertices,
-# interface45_vertices,
+# interface45_vertices_a,
+# interface45_vertices_b,
# interface46_vertices,
# interface56_vertices,
# ]
lambda_param = {
0: {'wetting': lambda12_w,
- 'nonwetting': lambda12_nw},#
+ 'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
- 'nonwetting': lambda23_nw},#
+ 'nonwetting': lambda23_nw},
2: {'wetting': lambda24_w,
- 'nonwetting': lambda24_nw},#
+ 'nonwetting': lambda24_nw},
3: {'wetting': lambda25_w,
- 'nonwetting': lambda25_nw},#
+ 'nonwetting': lambda25_nw},
4: {'wetting': lambda34_w,
- 'nonwetting': lambda34_nw},#
+ 'nonwetting': lambda34_nw},
5: {'wetting': lambda36_w,
- 'nonwetting': lambda36_nw},#
+ 'nonwetting': lambda36_nw},
6: {'wetting': lambda45_w,
- 'nonwetting': lambda45_nw},#
- 7: {'wetting': lambda46_w,
- 'nonwetting': lambda46_nw},#
- 8: {'wetting': lambda56_w,
- 'nonwetting': lambda56_nw},#
+ 'nonwetting': lambda45_nw},
+ 7: {'wetting': lambda45_w,
+ 'nonwetting': lambda45_nw},
+ 8: {'wetting': lambda46_w,
+ 'nonwetting': lambda46_nw},
+ 9: {'wetting': lambda56_w,
+ 'nonwetting': lambda56_nw},
}
# after Lewis, see pdf file
intrinsic_permeability = {
- 1: 0.1, # sand
- 2: 0.1, # sand, there is a range
- 3: 0.001, #10e-2, # clay has a range
- 4: 0.001, #10e-3
- 5: 0.001, #10e-2, # clay has a range
- 6: 0.001, #10e-3
+ 1: 0.01, # sand
+ 2: 0.01, # sand, there is a range
+ 3: 0.01, #10e-2, # clay has a range
+ 4: 0.01, #10e-3
+ 5: 0.01, #10e-2, # clay has a range
+ 6: 0.01, #10e-3
}
@@ -823,9 +856,9 @@ sat_pressure_relationship = {
6: ft.partial(saturation, n_index=2)
}
-#############################################
+###############################################################################
# Manufacture source expressions with sympy #
-#############################################
+###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
@@ -845,15 +878,14 @@ p_e_sym = {
'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
}
+
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
- pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- pc_e_sym.update(
- {subdomain: p_e_sym[subdomain]['nonwetting']
- - p_e_sym[subdomain]['wetting']}
- )
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
symbols = {"x": x,
@@ -880,6 +912,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
@@ -896,8 +929,7 @@ dirichletBC = dict()
# subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys():
- # if subdomain has no outer boundary, outer_boundary_def_points[subdomain]
- # is None
+ # subdomain can have no outer boundary
if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None})
else:
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py
index 89e336628299608941e9eec1acdf9e2dd4471543..0737454e106a9e38a8c24c6fe193871ea834db79 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py
@@ -13,7 +13,13 @@ import datetime
import os
import pandas as pd
-# check if output directory exists
+
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
@@ -21,58 +27,57 @@ else:
print("Directory ", './output', " already exists. Will use as output \
directory")
-
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
-# init sympy session
-sym.init_printing()
-# solver_tol = 6E-7
-use_case = "TP-R-layered-soil-with-inner-patch-all-params-one"
-# name of this very file. Needed for log output.
-thisfile = "TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py"
+# Name of the usecase that will be printed during simulation.
+use_case = "TP-R-layered-soil-with-inner-patch-realistic"
+# The name of this very file. Needed for creating log output.
+thisfile = "TP-R-layered_soil_with_inner_patch-realistic.py"
-max_iter_num = 150
+# GENERAL SOLVER CONFIG ######################################################
+# maximal iteration per timestep
+max_iter_num = 700
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False
resolutions = {
# 1: 2e-6, # h=2
# 2: 2e-6, # h=1.1180
# 4: 2e-6, # h=0.5590
# 8: 2e-6, # h=0.2814
- # 16: 7e-6, # h=0.1412
- 32: 5e-6,
- 64: 2e-6,
- 128: 2e-6,
- 256: 2e-6
+ 16: 8e-6, # h=0.1412
+ # 32: 2e-6,
+ # 64: 2e-6,
+ # 128: 2e-6
}
-# GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.00025
-number_of_timesteps = 4000
-plot_timestep_every = 4
-# 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 = 5
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
starttimes = [0.0]
+timestep_size = 0.01
+number_of_timesteps = 5
-Lw1 = 0.25 # /timestep_size
+# LDD scheme parameters ######################################################
+
+Lw1 = 0.025 # /timestep_size
Lnw1 = Lw1
-Lw2 = 0.25 # /timestep_size
+Lw2 = 0.025 # /timestep_size
Lnw2 = Lw2
-Lw3 = 0.25 # /timestep_size
+Lw3 = 0.025 # /timestep_size
Lnw3 = Lw3
-Lw4 = 0.25 # /timestep_size
+Lw4 = 0.025 # /timestep_size
Lnw4 = Lw4
-Lw5 = 0.25 # /timestep_size
+Lw5 = 0.025 # /timestep_size
Lnw5 = Lw5
-Lw6 = 0.25 # /timestep_size
+Lw6 = 0.025 # /timestep_size
Lnw6 = Lw6
lambda12_w = 4
@@ -102,45 +107,62 @@ lambda46_nw = 4
lambda56_w = 4
lambda56_nw = 4
-include_gravity = False
-debugflag = False
+include_gravity = True
+debugflag = True
analyse_condition = False
-# if mesh_study:
-output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-# else:
-# # for tol in resolutions.values():
-# # solver_tol = tol
-# output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
-# datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
-# )
-
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
'space_errornorms': True,
+ # save the mesh and marker functions to disk
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
+ 'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True
}
else:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ 'L_iterations_per_timestep': False,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'subsequent_errors': True
}
+# OUTPUT FILE STRING #########################################################
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+
+# DOMAIN AND INTERFACES #######################################################
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(1.0,-1.0),#
@@ -190,14 +212,18 @@ interface56_vertices = [interface46_vertices[1],
interface34_vertices = [interface36_vertices[1],
interface23_vertices[2]]
-# interface36
-interface45_vertices = [interface56_vertices[0],
- df.Point(0.7, -0.2),
+# Interface 45 needs to be split, because of the shape. There can be triangles
+# with two facets on the interface and this creates a rogue dof type error when
+# integrating over that particular interface. Accordingly, the lambda_param
+# dictionary has two entries for that interface.
+interface45_vertices_a = [interface56_vertices[0],
+ df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
+ ]
+interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
interface25_vertices[0]
]
-
# interface_vertices introduces a global numbering of interfaces.
interface_def_points = [interface12_vertices,
interface23_vertices,
@@ -205,11 +231,11 @@ interface_def_points = [interface12_vertices,
interface25_vertices,
interface34_vertices,
interface36_vertices,
- interface45_vertices,
+ interface45_vertices_a,
+ interface45_vertices_b,
interface46_vertices,
interface56_vertices,
]
-
adjacent_subdomains = [[1,2],
[2,3],
[2,4],
@@ -217,6 +243,7 @@ adjacent_subdomains = [[1,2],
[3,4],
[3,6],
[4,5],
+ [4,5],
[4,6],
[5,6]
]
@@ -226,18 +253,17 @@ subdomain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
interface12_vertices[2],
interface12_vertices[3],
- interface12_vertices[4], # southern boundary, 12 interface
- subdomain0_vertices[2], # eastern boundary, outer boundary
- subdomain0_vertices[3]] # northern boundary, outer on_boundary
+ interface12_vertices[4], # southern boundary, 12 interface
+ subdomain0_vertices[2], # eastern boundary, outer boundary
+ subdomain0_vertices[3]] # northern boundary, outer on_boundary
# vertex coordinates of the outer boundaries. If it can not be specified as a
-# polygon, use an entry per boundary polygon.
-# This information is used for defining
+# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
- 0: [subdomain1_vertices[4],
- subdomain1_vertices[5], # eastern boundary, outer boundary
+ 0: [subdomain1_vertices[4], #
+ subdomain1_vertices[5], # eastern boundary, outer boundary
subdomain1_vertices[6],
subdomain1_vertices[0]]
}
@@ -247,12 +273,12 @@ subdomain2_vertices = [interface23_vertices[0],
interface23_vertices[1],
interface23_vertices[2],
interface24_vertices[1],
- interface25_vertices[1], # southern boundary, 23 interface
- subdomain1_vertices[4], # eastern boundary, outer boundary
+ interface25_vertices[1], # southern boundary, 23 interface
+ subdomain1_vertices[4], # eastern boundary, outer boundary
subdomain1_vertices[3],
subdomain1_vertices[2],
subdomain1_vertices[1],
- subdomain1_vertices[0] ] # northern boundary, 12 interface
+ subdomain1_vertices[0] ] # northern boundary, 12 interface
subdomain2_outer_boundary_verts = {
0: [subdomain2_vertices[9],
@@ -280,7 +306,8 @@ subdomain3_outer_boundary_verts = {
# subdomain3
subdomain4_vertices = [interface46_vertices[0],
interface46_vertices[1],
- interface45_vertices[1],
+ # interface45_vertices[1],
+ interface45_vertices_a[1],
interface24_vertices[1],
interface24_vertices[0],
interface34_vertices[1]
@@ -293,10 +320,11 @@ subdomain5_vertices = [interface56_vertices[0],
interface56_vertices[2],
interface25_vertices[1],
interface25_vertices[0],
- interface45_vertices[1],
- interface45_vertices[0]
+ interface45_vertices_b[1],
+ interface45_vertices_b[0]
]
+
subdomain5_outer_boundary_verts = {
0: [subdomain5_vertices[2],
subdomain5_vertices[3]]
@@ -344,7 +372,7 @@ outer_boundary_def_points = {
6: subdomain6_outer_boundary_verts
}
-
+# MODEL CONFIGURATION #########################################################
isRichards = {
1: True,
2: True,
@@ -436,32 +464,36 @@ L = {
# interface25_vertices,
# interface34_vertices,
# interface36_vertices,
-# interface45_vertices,
+# interface45_vertices_a,
+# interface45_vertices_b,
# interface46_vertices,
# interface56_vertices,
# ]
lambda_param = {
0: {'wetting': lambda12_w,
- 'nonwetting': lambda12_nw},#
+ 'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
- 'nonwetting': lambda23_nw},#
+ 'nonwetting': lambda23_nw},
2: {'wetting': lambda24_w,
- 'nonwetting': lambda24_nw},#
+ 'nonwetting': lambda24_nw},
3: {'wetting': lambda25_w,
- 'nonwetting': lambda25_nw},#
+ 'nonwetting': lambda25_nw},
4: {'wetting': lambda34_w,
- 'nonwetting': lambda34_nw},#
+ 'nonwetting': lambda34_nw},
5: {'wetting': lambda36_w,
- 'nonwetting': lambda36_nw},#
+ 'nonwetting': lambda36_nw},
6: {'wetting': lambda45_w,
- 'nonwetting': lambda45_nw},#
- 7: {'wetting': lambda46_w,
- 'nonwetting': lambda46_nw},#
- 8: {'wetting': lambda56_w,
- 'nonwetting': lambda56_nw},#
+ 'nonwetting': lambda45_nw},
+ 7: {'wetting': lambda45_w,
+ 'nonwetting': lambda45_nw},
+ 8: {'wetting': lambda46_w,
+ 'nonwetting': lambda46_nw},
+ 9: {'wetting': lambda56_w,
+ 'nonwetting': lambda56_nw},
}
+# after Lewis, see pdf file
# after Lewis, see pdf file
intrinsic_permeability = {
1: 1, # sand
@@ -825,37 +857,36 @@ sat_pressure_relationship = {
6: ft.partial(saturation, n_index=2)
}
-#############################################
+###############################################################################
# Manufacture source expressions with sympy #
-#############################################
+###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym = {
- 1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
- 'nonwetting': 0.0*t },
- 2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
- 'nonwetting': 0.0*t },
+ 1: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 'nonwetting': 0*t },
+ 2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 'nonwetting': 0*t },
3: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
4: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
5: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
6: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
}
+
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
- pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- pc_e_sym.update(
- {subdomain: p_e_sym[subdomain]['nonwetting']
- - p_e_sym[subdomain]['wetting']}
- )
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
symbols = {"x": x,
@@ -882,6 +913,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
@@ -898,8 +930,7 @@ dirichletBC = dict()
# subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys():
- # if subdomain has no outer boundary, outer_boundary_def_points[subdomain]
- # is None
+ # subdomain can have no outer boundary
if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None})
else:
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py
index 9b11c1aff69216de780b2b5f0818f1bb01b18b44..62b02c79ed8d73487db951758b813a72ab1f5c92 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py
@@ -13,7 +13,13 @@ import datetime
import os
import pandas as pd
-# check if output directory exists
+
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
@@ -21,19 +27,19 @@ else:
print("Directory ", './output', " already exists. Will use as output \
directory")
-
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
-# init sympy session
-sym.init_printing()
-# solver_tol = 6E-7
+# Name of the usecase that will be printed during simulation.
use_case = "TP-R-layered-soil-with-inner-patch-pure-dd"
-# name of this very file. Needed for log output.
+# The name of this very file. Needed for creating log output.
thisfile = "TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py"
+# GENERAL SOLVER CONFIG ######################################################
max_iter_num = 300
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = True
resolutions = {
# 1: 5e-6, # h=2
@@ -46,16 +52,14 @@ resolutions = {
128: 3e-6
}
-# GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.001
-number_of_timesteps = 800
-plot_timestep_every = 4
-# 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 = 5
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
starttimes = [0.0]
+timestep_size = 0.01
+number_of_timesteps = 5
+# LDD scheme parameters ######################################################
Lw1 = 0.25 # /timestep_size
Lnw1 = Lw1
@@ -101,45 +105,62 @@ lambda46_nw = 4
lambda56_w = 4
lambda56_nw = 4
-include_gravity = False
-debugflag = False
+include_gravity = True
+debugflag = True
analyse_condition = False
-if mesh_study:
- output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-else:
- for tol in resolutions.values():
- solver_tol = tol
- output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
- )
-
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
'space_errornorms': True,
+ # save the mesh and marker functions to disk
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
+ 'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True
}
else:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ 'L_iterations_per_timestep': False,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'subsequent_errors': True
}
+# OUTPUT FILE STRING #########################################################
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+
+# DOMAIN AND INTERFACES #######################################################
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(1.0,-1.0),#
@@ -189,14 +210,18 @@ interface56_vertices = [interface46_vertices[1],
interface34_vertices = [interface36_vertices[1],
interface23_vertices[2]]
-# interface36
-interface45_vertices = [interface56_vertices[0],
- df.Point(0.7, -0.2),
+# Interface 45 needs to be split, because of the shape. There can be triangles
+# with two facets on the interface and this creates a rogue dof type error when
+# integrating over that particular interface. Accordingly, the lambda_param
+# dictionary has two entries for that interface.
+interface45_vertices_a = [interface56_vertices[0],
+ df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
+ ]
+interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
interface25_vertices[0]
]
-
# interface_vertices introduces a global numbering of interfaces.
interface_def_points = [interface12_vertices,
interface23_vertices,
@@ -204,11 +229,11 @@ interface_def_points = [interface12_vertices,
interface25_vertices,
interface34_vertices,
interface36_vertices,
- interface45_vertices,
+ interface45_vertices_a,
+ interface45_vertices_b,
interface46_vertices,
interface56_vertices,
]
-
adjacent_subdomains = [[1,2],
[2,3],
[2,4],
@@ -216,6 +241,7 @@ adjacent_subdomains = [[1,2],
[3,4],
[3,6],
[4,5],
+ [4,5],
[4,6],
[5,6]
]
@@ -225,18 +251,17 @@ subdomain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
interface12_vertices[2],
interface12_vertices[3],
- interface12_vertices[4], # southern boundary, 12 interface
- subdomain0_vertices[2], # eastern boundary, outer boundary
- subdomain0_vertices[3]] # northern boundary, outer on_boundary
+ interface12_vertices[4], # southern boundary, 12 interface
+ subdomain0_vertices[2], # eastern boundary, outer boundary
+ subdomain0_vertices[3]] # northern boundary, outer on_boundary
# vertex coordinates of the outer boundaries. If it can not be specified as a
-# polygon, use an entry per boundary polygon.
-# This information is used for defining
+# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
- 0: [subdomain1_vertices[4],
- subdomain1_vertices[5], # eastern boundary, outer boundary
+ 0: [subdomain1_vertices[4], #
+ subdomain1_vertices[5], # eastern boundary, outer boundary
subdomain1_vertices[6],
subdomain1_vertices[0]]
}
@@ -246,12 +271,12 @@ subdomain2_vertices = [interface23_vertices[0],
interface23_vertices[1],
interface23_vertices[2],
interface24_vertices[1],
- interface25_vertices[1], # southern boundary, 23 interface
- subdomain1_vertices[4], # eastern boundary, outer boundary
+ interface25_vertices[1], # southern boundary, 23 interface
+ subdomain1_vertices[4], # eastern boundary, outer boundary
subdomain1_vertices[3],
subdomain1_vertices[2],
subdomain1_vertices[1],
- subdomain1_vertices[0] ] # northern boundary, 12 interface
+ subdomain1_vertices[0] ] # northern boundary, 12 interface
subdomain2_outer_boundary_verts = {
0: [subdomain2_vertices[9],
@@ -279,7 +304,8 @@ subdomain3_outer_boundary_verts = {
# subdomain3
subdomain4_vertices = [interface46_vertices[0],
interface46_vertices[1],
- interface45_vertices[1],
+ # interface45_vertices[1],
+ interface45_vertices_a[1],
interface24_vertices[1],
interface24_vertices[0],
interface34_vertices[1]
@@ -292,10 +318,11 @@ subdomain5_vertices = [interface56_vertices[0],
interface56_vertices[2],
interface25_vertices[1],
interface25_vertices[0],
- interface45_vertices[1],
- interface45_vertices[0]
+ interface45_vertices_b[1],
+ interface45_vertices_b[0]
]
+
subdomain5_outer_boundary_verts = {
0: [subdomain5_vertices[2],
subdomain5_vertices[3]]
@@ -343,7 +370,7 @@ outer_boundary_def_points = {
6: subdomain6_outer_boundary_verts
}
-
+# MODEL CONFIGURATION #########################################################
isRichards = {
1: True,
2: True,
@@ -435,29 +462,32 @@ L = {
# interface25_vertices,
# interface34_vertices,
# interface36_vertices,
-# interface45_vertices,
+# interface45_vertices_a,
+# interface45_vertices_b,
# interface46_vertices,
# interface56_vertices,
# ]
lambda_param = {
0: {'wetting': lambda12_w,
- 'nonwetting': lambda12_nw},#
+ 'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
- 'nonwetting': lambda23_nw},#
+ 'nonwetting': lambda23_nw},
2: {'wetting': lambda24_w,
- 'nonwetting': lambda24_nw},#
+ 'nonwetting': lambda24_nw},
3: {'wetting': lambda25_w,
- 'nonwetting': lambda25_nw},#
+ 'nonwetting': lambda25_nw},
4: {'wetting': lambda34_w,
- 'nonwetting': lambda34_nw},#
+ 'nonwetting': lambda34_nw},
5: {'wetting': lambda36_w,
- 'nonwetting': lambda36_nw},#
+ 'nonwetting': lambda36_nw},
6: {'wetting': lambda45_w,
- 'nonwetting': lambda45_nw},#
- 7: {'wetting': lambda46_w,
- 'nonwetting': lambda46_nw},#
- 8: {'wetting': lambda56_w,
- 'nonwetting': lambda56_nw},#
+ 'nonwetting': lambda45_nw},
+ 7: {'wetting': lambda45_w,
+ 'nonwetting': lambda45_nw},
+ 8: {'wetting': lambda46_w,
+ 'nonwetting': lambda46_nw},
+ 9: {'wetting': lambda56_w,
+ 'nonwetting': lambda56_nw},
}
@@ -824,9 +854,9 @@ sat_pressure_relationship = {
6: ft.partial(saturation, n_index=2)
}
-#############################################
+###############################################################################
# Manufacture source expressions with sympy #
-#############################################
+###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
@@ -849,12 +879,10 @@ p_e_sym = {
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
- pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- pc_e_sym.update(
- {subdomain: p_e_sym[subdomain]['nonwetting']
- - p_e_sym[subdomain]['wetting']}
- )
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
symbols = {"x": x,
@@ -881,6 +909,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
@@ -897,8 +926,7 @@ dirichletBC = dict()
# subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys():
- # if subdomain has no outer boundary, outer_boundary_def_points[subdomain]
- # is None
+ # subdomain can have no outer boundary
if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None})
else:
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one-finer-dt.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-realistic.py
similarity index 79%
rename from Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one-finer-dt.py
rename to Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-realistic.py
index ddd85b3b4aa53a5a517af046631c1fbe95671323..8d330bf9bb1cd9fab6afcb08ecacce858304a99c 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one-finer-dt.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -1,16 +1,3 @@
-Directory ./output created
- type of dtpc is <class 'sympy.core.mul.Mul'>
- type of dS is Piecewise
- type of dtpc is <class 'sympy.core.mul.Mul'>
- type of dS is Piecewise
- type of dtpc is <class 'sympy.core.add.Add'>
- type of dS is Piecewise
- type of dtpc is <class 'sympy.core.add.Add'>
- type of dS is Piecewise
- type of dtpc is <class 'sympy.core.add.Add'>
- type of dS is Piecewise
- type of dtpc is <class 'sympy.core.add.Add'>
- type of dS is Piecewise
#!/usr/bin/python3
"""Layered soil simulation with inner patch.
@@ -26,7 +13,13 @@ import datetime
import os
import pandas as pd
-# check if output directory exists
+
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
@@ -34,57 +27,57 @@ else:
print("Directory ", './output', " already exists. Will use as output \
directory")
-
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
-# init sympy session
-sym.init_printing()
-# solver_tol = 6E-7
-use_case = "TP-R-layered-soil-with-inner-patch-all-params-one"
-# name of this very file. Needed for log output.
-thisfile = "TP-R-layered_soil_with_inner_patch-all-params-one-finer-dt.py"
+# Name of the usecase that will be printed during simulation.
+use_case = "TP-R-layered-soil-with-inner-patch-realistic"
+# The name of this very file. Needed for creating log output.
+thisfile = "TP-R-layered_soil_with_inner_patch-realistic.py"
-max_iter_num = 150
+# GENERAL SOLVER CONFIG ######################################################
+# maximal iteration per timestep
+max_iter_num = 700
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False
resolutions = {
# 1: 2e-6, # h=2
# 2: 2e-6, # h=1.1180
# 4: 2e-6, # h=0.5590
# 8: 2e-6, # h=0.2814
- # 16: 7e-6, # h=0.1412
- # 32: 5e-6,
- 64: 2e-6,
+ 16: 8e-6, # h=0.1412
+ # 32: 2e-6,
+ # 64: 2e-6,
# 128: 2e-6
}
-# GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.00025
-number_of_timesteps = 4000
-plot_timestep_every = 4
-# 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 = 5
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
starttimes = [0.0]
+timestep_size = 0.01
+number_of_timesteps = 5
-Lw1 = 0.25 # /timestep_size
+# LDD scheme parameters ######################################################
+
+Lw1 = 0.025 # /timestep_size
Lnw1 = Lw1
-Lw2 = 0.25 # /timestep_size
+Lw2 = 0.025 # /timestep_size
Lnw2 = Lw2
-Lw3 = 0.25 # /timestep_size
+Lw3 = 0.025 # /timestep_size
Lnw3 = Lw3
-Lw4 = 0.25 # /timestep_size
+Lw4 = 0.025 # /timestep_size
Lnw4 = Lw4
-Lw5 = 0.25 # /timestep_size
+Lw5 = 0.025 # /timestep_size
Lnw5 = Lw5
-Lw6 = 0.25 # /timestep_size
+Lw6 = 0.025 # /timestep_size
Lnw6 = Lw6
lambda12_w = 4
@@ -114,45 +107,62 @@ lambda46_nw = 4
lambda56_w = 4
lambda56_nw = 4
-include_gravity = False
-debugflag = False
+include_gravity = True
+debugflag = True
analyse_condition = False
-# if mesh_study:
-output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-# else:
-# # for tol in resolutions.values():
-# # solver_tol = tol
-# output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
-# datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
-# )
-
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
'space_errornorms': True,
+ # save the mesh and marker functions to disk
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
+ 'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True
}
else:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
+ 'L_iterations_per_timestep': False,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'subsequent_errors': True
}
+# OUTPUT FILE STRING #########################################################
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+
+# DOMAIN AND INTERFACES #######################################################
# global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), #
df.Point(1.0,-1.0),#
@@ -202,14 +212,18 @@ interface56_vertices = [interface46_vertices[1],
interface34_vertices = [interface36_vertices[1],
interface23_vertices[2]]
-# interface36
-interface45_vertices = [interface56_vertices[0],
- df.Point(0.7, -0.2),
+# Interface 45 needs to be split, because of the shape. There can be triangles
+# with two facets on the interface and this creates a rogue dof type error when
+# integrating over that particular interface. Accordingly, the lambda_param
+# dictionary has two entries for that interface.
+interface45_vertices_a = [interface56_vertices[0],
+ df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
+ ]
+interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
interface25_vertices[0]
]
-
# interface_vertices introduces a global numbering of interfaces.
interface_def_points = [interface12_vertices,
interface23_vertices,
@@ -217,11 +231,11 @@ interface_def_points = [interface12_vertices,
interface25_vertices,
interface34_vertices,
interface36_vertices,
- interface45_vertices,
+ interface45_vertices_a,
+ interface45_vertices_b,
interface46_vertices,
interface56_vertices,
]
-
adjacent_subdomains = [[1,2],
[2,3],
[2,4],
@@ -229,6 +243,7 @@ adjacent_subdomains = [[1,2],
[3,4],
[3,6],
[4,5],
+ [4,5],
[4,6],
[5,6]
]
@@ -238,18 +253,17 @@ subdomain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
interface12_vertices[2],
interface12_vertices[3],
- interface12_vertices[4], # southern boundary, 12 interface
- subdomain0_vertices[2], # eastern boundary, outer boundary
- subdomain0_vertices[3]] # northern boundary, outer on_boundary
+ interface12_vertices[4], # southern boundary, 12 interface
+ subdomain0_vertices[2], # eastern boundary, outer boundary
+ subdomain0_vertices[3]] # northern boundary, outer on_boundary
# vertex coordinates of the outer boundaries. If it can not be specified as a
-# polygon, use an entry per boundary polygon.
-# This information is used for defining
+# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
- 0: [subdomain1_vertices[4],
- subdomain1_vertices[5], # eastern boundary, outer boundary
+ 0: [subdomain1_vertices[4], #
+ subdomain1_vertices[5], # eastern boundary, outer boundary
subdomain1_vertices[6],
subdomain1_vertices[0]]
}
@@ -259,12 +273,12 @@ subdomain2_vertices = [interface23_vertices[0],
interface23_vertices[1],
interface23_vertices[2],
interface24_vertices[1],
- interface25_vertices[1], # southern boundary, 23 interface
- subdomain1_vertices[4], # eastern boundary, outer boundary
+ interface25_vertices[1], # southern boundary, 23 interface
+ subdomain1_vertices[4], # eastern boundary, outer boundary
subdomain1_vertices[3],
subdomain1_vertices[2],
subdomain1_vertices[1],
- subdomain1_vertices[0] ] # northern boundary, 12 interface
+ subdomain1_vertices[0] ] # northern boundary, 12 interface
subdomain2_outer_boundary_verts = {
0: [subdomain2_vertices[9],
@@ -292,7 +306,8 @@ subdomain3_outer_boundary_verts = {
# subdomain3
subdomain4_vertices = [interface46_vertices[0],
interface46_vertices[1],
- interface45_vertices[1],
+ # interface45_vertices[1],
+ interface45_vertices_a[1],
interface24_vertices[1],
interface24_vertices[0],
interface34_vertices[1]
@@ -305,10 +320,11 @@ subdomain5_vertices = [interface56_vertices[0],
interface56_vertices[2],
interface25_vertices[1],
interface25_vertices[0],
- interface45_vertices[1],
- interface45_vertices[0]
+ interface45_vertices_b[1],
+ interface45_vertices_b[0]
]
+
subdomain5_outer_boundary_verts = {
0: [subdomain5_vertices[2],
subdomain5_vertices[3]]
@@ -356,7 +372,7 @@ outer_boundary_def_points = {
6: subdomain6_outer_boundary_verts
}
-
+# MODEL CONFIGURATION #########################################################
isRichards = {
1: True,
2: True,
@@ -378,33 +394,33 @@ isRichards = {
# Dict of the form: { subdom_num : viscosity }
viscosity = {
1: {'wetting' :1,
- 'nonwetting': 1},
+ 'nonwetting': 1/50},
2: {'wetting' :1,
- 'nonwetting': 1},
+ 'nonwetting': 1/50},
3: {'wetting' :1,
- 'nonwetting': 1},
+ 'nonwetting': 1/50},
4: {'wetting' :1,
- 'nonwetting': 1},
+ 'nonwetting': 1/50},
5: {'wetting' :1,
- 'nonwetting': 1},
+ 'nonwetting': 1/50},
6: {'wetting' :1,
- 'nonwetting': 1},
+ 'nonwetting': 1/50},
}
# Dict of the form: { subdom_num : density }
densities = {
- 1: {'wetting': 1, #997
- 'nonwetting': 1}, #1}, #1.225},
- 2: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225},
- 3: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225},
- 4: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225}
- 5: {'wetting': 1, #997
- 'nonwetting': 1}, #1.225},
- 6: {'wetting': 1, #997
- 'nonwetting': 1} #1.225}
+ 1: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1}, #1.225},
+ 2: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225},
+ 3: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225},
+ 4: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225}
+ 5: {'wetting': 997.0, #997
+ 'nonwetting': 1.225}, #1.225},
+ 6: {'wetting': 997.0, #997
+ 'nonwetting': 1.225} #1.225}
}
gravity_acceleration = 9.81
@@ -412,12 +428,12 @@ gravity_acceleration = 9.81
# https://www.geotechdata.info/parameter/soil-porosity.html
# Dict of the form: { subdom_num : porosity }
porosity = {
- 1: 1, #0.2, # Clayey gravels, clayey sandy gravels
- 2: 1, #0.22, # Silty gravels, silty sandy gravels
- 3: 1, #0.37, # Clayey sands
- 4: 1, #0.2 # Silty or sandy clay
- 5: 1, #
- 6: 1, #
+ 1: 0.2, #0.2, # Clayey gravels, clayey sandy gravels
+ 2: 0.2, #0.22, # Silty gravels, silty sandy gravels
+ 3: 0.2, #0.37, # Clayey sands
+ 4: 0.2, #0.2 # Silty or sandy clay
+ 5: 0.2, #
+ 6: 0.2, #
}
# subdom_num : subdomain L for L-scheme
@@ -448,40 +464,43 @@ L = {
# interface25_vertices,
# interface34_vertices,
# interface36_vertices,
-# interface45_vertices,
+# interface45_vertices_a,
+# interface45_vertices_b,
# interface46_vertices,
# interface56_vertices,
# ]
lambda_param = {
0: {'wetting': lambda12_w,
- 'nonwetting': lambda12_nw},#
+ 'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
- 'nonwetting': lambda23_nw},#
+ 'nonwetting': lambda23_nw},
2: {'wetting': lambda24_w,
- 'nonwetting': lambda24_nw},#
+ 'nonwetting': lambda24_nw},
3: {'wetting': lambda25_w,
- 'nonwetting': lambda25_nw},#
+ 'nonwetting': lambda25_nw},
4: {'wetting': lambda34_w,
- 'nonwetting': lambda34_nw},#
+ 'nonwetting': lambda34_nw},
5: {'wetting': lambda36_w,
- 'nonwetting': lambda36_nw},#
+ 'nonwetting': lambda36_nw},
6: {'wetting': lambda45_w,
- 'nonwetting': lambda45_nw},#
- 7: {'wetting': lambda46_w,
- 'nonwetting': lambda46_nw},#
- 8: {'wetting': lambda56_w,
- 'nonwetting': lambda56_nw},#
+ 'nonwetting': lambda45_nw},
+ 7: {'wetting': lambda45_w,
+ 'nonwetting': lambda45_nw},
+ 8: {'wetting': lambda46_w,
+ 'nonwetting': lambda46_nw},
+ 9: {'wetting': lambda56_w,
+ 'nonwetting': lambda56_nw},
}
# after Lewis, see pdf file
intrinsic_permeability = {
- 1: 1, # sand
- 2: 1, # sand, there is a range
- 3: 1, #10e-2, # clay has a range
- 4: 1, #10e-3
- 5: 1, #10e-2, # clay has a range
- 6: 1, #10e-3
+ 1: 0.01, # sand
+ 2: 0.01, # sand, there is a range
+ 3: 0.01, #10e-2, # clay has a range
+ 4: 0.01, #10e-3
+ 5: 0.01, #10e-2, # clay has a range
+ 6: 0.01, #10e-3
}
@@ -837,37 +856,36 @@ sat_pressure_relationship = {
6: ft.partial(saturation, n_index=2)
}
-#############################################
+###############################################################################
# Manufacture source expressions with sympy #
-#############################################
+###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym = {
- 1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
- 'nonwetting': 0.0*t },
- 2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
- 'nonwetting': 0.0*t },
+ 1: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 'nonwetting': 0*t },
+ 2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 'nonwetting': 0*t },
3: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
4: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
5: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
6: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)),
- 'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2)*y**2 },
+ 'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
}
+
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
- pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
- pc_e_sym.update(
- {subdomain: p_e_sym[subdomain]['nonwetting']
- - p_e_sym[subdomain]['wetting']}
- )
+ pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
+ - p_e_sym[subdomain]['wetting'].copy()})
symbols = {"x": x,
@@ -894,6 +912,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
@@ -910,8 +929,7 @@ dirichletBC = dict()
# subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys():
- # if subdomain has no outer boundary, outer_boundary_def_points[subdomain]
- # is None
+ # subdomain can have no outer boundary
if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None})
else:
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 5c439052ea5a8235fd6ed909ae5f74426c57c088..a3e6306eaa8d6fed2c84d2344d74ea4074fe7fb5 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -120,16 +120,9 @@ else:
}
# OUTPUT FILE STRING #########################################################
-if mesh_study:
- output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-else:
- for tol in resolutions.values():
- solver_tol = tol
- output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
- )
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
# DOMAIN AND INTERFACE #######################################################
diff --git a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
index d1fcfbc0da9eef1a902af8b6e0170540b746ae4f..7f57d615342ea219370704f223c714c9fe7df600 100755
--- a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
+++ b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
@@ -1,5 +1,5 @@
#!/usr/bin/python3
-"""TP-TP 2 patch soil simulation.
+"""TP-TP layered soil with inner patch simulation.
This program sets up an LDD simulation
"""
@@ -47,7 +47,7 @@ resolutions = {
# 4: 1e-6,
# 8: 1e-6,
# 16: 5e-6,
- 32: 3e-6,
+ 32: 8e-6,
# 64: 2e-6,
# 128: 1e-6,
# 256: 1e-6,
@@ -56,65 +56,65 @@ resolutions = {
# starttimes gives a list of starttimes to run the simulation from.
# The list is looped over and a simulation is run with t_0 as initial time
# for each element t_0 in starttimes.
-starttimes = [0.0]
+starttimes = [0.3]
timestep_size = 0.001
-number_of_timesteps = 1000
+number_of_timesteps = 5
# LDD scheme parameters ######################################################
-Lw1 = 0.25 # /timestep_size
-Lnw1 = Lw1
+Lw1 = 0.5 # /timestep_size
+Lnw1 = 0.025
-Lw2 = 0.25 # /timestep_size
-Lnw2 = Lw2
+Lw2 = 0.5 # /timestep_size
+Lnw2 = 0.025
-Lw3 = 0.05 # /timestep_size
-Lnw3 = Lw3
+Lw3 = 0.5 # /timestep_size
+Lnw3 = 0.025
-Lw4 = 0.05 # /timestep_size
-Lnw4 = Lw4
+Lw4 = 0.5 # /timestep_size
+Lnw4 = 0.025
-Lw5 = 0.05 # /timestep_size
-Lnw5 = Lw5
+Lw5 = 0.5 # /timestep_size
+Lnw5 = 0.025
-Lw6 = 0.05 # /timestep_size
-Lnw6 = Lw6
+Lw6 = 0.5 # /timestep_size
+Lnw6 = 0.025
-lambda12_w = 40
-lambda12_nw = 40
+lambda12_w = 4
+lambda12_nw = 4
-lambda23_w = 40
-lambda23_nw = 40
+lambda23_w = 4
+lambda23_nw = 4
-lambda24_w = 40
-lambda24_nw= 40
+lambda24_w = 4
+lambda24_nw= 4
-lambda25_w= 40
-lambda25_nw= 40
+lambda25_w= 4
+lambda25_nw= 4
-lambda34_w = 40
-lambda34_nw = 40
+lambda34_w = 4
+lambda34_nw = 4
-lambda36_w = 40
-lambda36_nw = 40
+lambda36_w = 4
+lambda36_nw = 4
-lambda45_w = 40
-lambda45_nw = 40
+lambda45_w = 4
+lambda45_nw = 4
-lambda46_w = 40
-lambda46_nw = 40
+lambda46_w = 4
+lambda46_nw = 4
-lambda56_w = 40
-lambda56_nw = 40
+lambda56_w = 4
+lambda56_nw = 4
-include_gravity = False
-debugflag = False
+include_gravity = True
+debugflag = True
analyse_condition = False
# I/O CONFIG #################################################################
# when number_of_timesteps is high, it might take a long time to write all
# timesteps to disk. Therefore, you can choose to only write data of every
# plot_timestep_every timestep to disk.
-plot_timestep_every = 4
+plot_timestep_every = 1
# Decide how many timesteps you want analysed. Analysed means, that
# subsequent errors of the L-iteration within the timestep are written out.
number_of_timesteps_to_analyse = 5
@@ -156,16 +156,9 @@ else:
}
# OUTPUT FILE STRING #########################################################
-if mesh_study:
- output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-else:
- for tol in resolutions.values():
- solver_tol = tol
- output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
- )
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
# DOMAIN AND INTERFACE #######################################################