Skip to content
Snippets Groups Projects
Commit a6c6b8bf authored by David Seus's avatar David Seus
Browse files

adjust readme, bring TP-R-2-patch-test.py up to par.

parent 2a5e9934
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Usecases/README.md
+ 254
− 1
View file @ a6c6b8bf
...@@ -69,3 +69,256 @@ simulation data. Chose a descriptive name for the usecase you are simulating.! ...@@ -69,3 +69,256 @@ simulation data. Chose a descriptive name for the usecase you are simulating.!
output the content of the script to stdout and into the logfile if the script output the content of the script to stdout and into the logfile if the script
`run_simulation` is used to start the simulation. The latter should be `run_simulation` is used to start the simulation. The latter should be
contained in all folders. contained in all folders.
### SOLVER CONFIG
~~~python
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
max_iter_num = 500
FEM_Lagrange_degree = 1
~~~
This block sets the maximal number of L-iterations `max_iter_num` that the
LDD solver gets for the calculation of each time step.
`FEM_Lagrange_degree` sets the degree of FEM ansatz functions that are used.
While the communication of dofs over interfaces has been writen in a way to
allow for general anasatz functions, other parts of the code have not.
If higher order FEM methods are to be used, the code will have to be adjusted.
### GRID AND MESH STUDY CONFIG
~~~python
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = True
resolutions = {
# 1: 1e-6,
# 2: 1e-6,
# 4: 1e-6,
# 8: 1e-6,
# 16: 5e-6,
# 32: 5e-6,
64: 2e-6,
# 128: 1e-6,
# 256: 1e-6,
}
# 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.0}
timestep_size = 0.001
number_of_timesteps = 20
~~~
This block sets mesh parameters, number of timesteps to be simulated, timestep
size as well as whether or not a mesh study is being performed.
In detail:
- ~python mesh_study = True/False~ Toggle mesh study mode.
This changes the output slightly. Have a look at examples in the `mesh_study`
folders.
- resolutions = {
# 1: 1e-6,
# 2: 1e-6,
# 4: 1e-6,
# 8: 1e-6,
# 16: 5e-6,
# 32: 5e-6,
64: 2e-6,
# 128: 1e-6,
# 256: 1e-6,
}
- # 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.0}
- timestep_size = 0.001
- number_of_timesteps = 20
### LDD SCHEME PARAMATERS
~~~python
# LDD scheme parameters ######################################################
Lw1 = 0.25 #/timestep_size
Lnw1= 0.25
Lw2 = 0.5 #/timestep_size
Lnw2= 0.25
lambda_w = 40
lambda_nw = 40
include_gravity = True
debugflag = True
analyse_condition = False
~~~
### I/O CONFIG
~~~python
# 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
# 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:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
'L_iterations_per_timestep': False,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'subsequent_errors': True
}
~~~
### OUTPUT FILE STRING
~~~python
# OUTPUT FILE STRING #########################################################
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
~~~
### DOMAIN SUBSTRUCTURING
~~~python
# DOMAIN AND INTERFACE #######################################################
substructuring = dss.twoSoilLayers()
interface_def_points = substructuring.interface_def_points
adjacent_subdomains = substructuring.adjacent_subdomains
subdomain_def_points = substructuring.subdomain_def_points
outer_boundary_def_points = substructuring.outer_boundary_def_points
~~
### MODEL CONFIGURATION
~~~python
# MODEL CONFIGURATION #########################################################
isRichards = {
1: True, #
2: False
}
viscosity = {#
# subdom_num : viscosity
1: {'wetting' :1,
'nonwetting': 1/50}, #
2: {'wetting' :1,
'nonwetting': 1/50}
}
porosity = {#
# subdom_num : porosity
1: 0.22,#
2: 0.22
}
# Dict of the form: { subdom_num : density }
densities = {
1: {'wetting': 997,
'nonwetting': 1.225},
2: {'wetting': 997,
'nonwetting': 1.225}
}
gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
1 : {'wetting' :Lw1,
'nonwetting': Lnw1},#
2 : {'wetting' :Lw2,
'nonwetting': Lnw2}
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
0 : {'wetting' :lambda_w,
'nonwetting': lambda_nw},#
}
intrinsic_permeability = {
1: 0.01,
2: 0.01,
}
# relative permeabilties
rel_perm_definition = {
1: {"wetting": "Spow2",
"nonwetting": "oneMinusSpow2"},
2: {"wetting": "Spow3",
"nonwetting": "oneMinusSpow3"},
}
rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
relative_permeability = rel_perm_dict["ka"]
ka_prime = rel_perm_dict["ka_prime"]
# S-pc relation
Spc_on_subdomains = {
1: {"testSpc": {"index": 1}},
2: {"testSpc": {"index": 2}},
}
Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
S_pc_sym = Spc["symbolic"]
S_pc_sym_prime = Spc["prime_symbolic"]
sat_pressure_relationship = Spc["dolfin"]
~~~
### MANUFACTURED SOLUTION/SOURCE EXPRESSION
~~~python
###############################################################################
# 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))}, #*(1-x)**2*(1+x)**2*(1-y)**2},
2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
'nonwetting': (-2-t*(1.1+y + x**2))*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)
~~~
### BOUNDARY CONDITIONS
~~~python
# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions. If an exact solution case is
# used, use the hlp.generate_exact_DirichletBC() method to generate the
# Dirichlet Boundary conditions from the exact solution.
dirichletBC = hlp.generate_exact_DirichletBC(
isRichards=isRichards,
outer_boundary_def_points=outer_boundary_def_points,
exact_solution=exact_solution
)
# If no exact solution is provided you need to provide a dictionary of boundary
# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see
# the structure.
~~~
Usecases/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+ 17
− 166
View file @ a6c6b8bf
...@@ -201,177 +201,28 @@ intrinsic_permeability = { ...@@ -201,177 +201,28 @@ intrinsic_permeability = {
2: 0.01, 2: 0.01,
} }
## relative permeabilty functions on subdomain 1 # relative permeabilties
def rel_perm1w(s): rel_perm_definition = {
# relative permeabilty wetting on subdomain1 1: {"wetting": "Spow2",
return intrinsic_permeability[1]*s**2 "nonwetting": "oneMinusSpow2"},
2: {"wetting": "Spow3",
def rel_perm1nw(s): "nonwetting": "oneMinusSpow3"},
# relative permeabilty nonwetting on subdomain1
return intrinsic_permeability[1]*(1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
subdomain1_rel_perm = {
'wetting': _rel_perm1w,#
'nonwetting': _rel_perm1nw
}
## relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
return intrinsic_permeability[2]*s**3
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdomain2
return intrinsic_permeability[2]*(1-s)**3
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)
subdomain2_rel_perm = {
'wetting': _rel_perm2w,#
'nonwetting': _rel_perm2nw
}
## dictionary of relative permeabilties on all domains.
relative_permeability = {#
1: subdomain1_rel_perm,
2: subdomain2_rel_perm
} }
rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
relative_permeability = rel_perm_dict["ka"]
ka_prime = rel_perm_dict["ka_prime"]
# definition of the derivatives of the relative permeabilities # S-pc relation
# relative permeabilty functions on subdomain 1 Spc_on_subdomains = {
def rel_perm1w_prime(s): 1: {"testSpc": {"index": 1}},
# relative permeabilty on subdomain1 2: {"testSpc": {"index": 2}},
return intrinsic_permeability[1]*2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
return -1*intrinsic_permeability[1]*2*(1-s)
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
def rel_perm2w_prime(s):
# relative permeabilty on subdomain2
return intrinsic_permeability[2]*3*s**2
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain2
return -3*intrinsic_permeability[2]*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
subdomain1_rel_perm_prime = {
'wetting': _rel_perm1w_prime,
'nonwetting': _rel_perm1nw_prime
}
subdomain2_rel_perm_prime = {
'wetting': _rel_perm2w_prime,
'nonwetting': _rel_perm2nw_prime
}
# dictionary of relative permeabilties on all domains.
ka_prime = {
1: subdomain1_rel_perm_prime,
2: subdomain2_rel_perm_prime,
}
# def saturation1(pc, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
#
# def saturation2(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
#
# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
# def saturation1_sym(pc, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return 1/((1 + pc)**(1/(subdomain_index + 1)))
#
#
# def saturation2_sym(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# #df.conditional(pc > 0,
# return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
#
#
# # derivative of S-pc relationship with respect to pc. This is needed for the
# # construction of a analytic solution.
# def saturation1_sym_prime(pc, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
#
#
# def saturation2_sym_prime(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
#
# # note that the conditional definition of S-pc in the nonsymbolic part will be
# # incorporated in the construction of the exact solution below.
# S_pc_sym = {
# 1: ft.partial(saturation1_sym, subdomain_index = 1),
# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
# }
#
# S_pc_sym_prime = {
# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
# }
#
# sat_pressure_relationship = {
# 1: ft.partial(saturation1, subdomain_index = 1),#,
# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
# }
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)))
# 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 -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
# note that the conditional definition of S-pc in the nonsymbolic part will be
# incorporated in the construction of the exact solution below.
S_pc_sym = {
1: ft.partial(saturation_sym, index=1),
2: ft.partial(saturation_sym, index=2),
# 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),
# 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),
# 3: ft.partial(saturation, index=2),
# 4: ft.partial(saturation, index=1)
} }
Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
S_pc_sym = Spc["symbolic"]
S_pc_sym_prime = Spc["prime_symbolic"]
sat_pressure_relationship = Spc["dolfin"]
############################################################################### ###############################################################################
# Manufacture source expressions with sympy # # Manufacture source expressions with sympy #
... ...
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment