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

add TPR inner patch mesh study

parent 567ae212
No related branches found
No related tags found
No related merge requests found
#!/usr/bin/python3
"""Layered soil simulation with inner patch.
This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
import functions as fts
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
import multiprocessing as mp
import domainSubstructuring as dss
# 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 ")
else:
print("Directory ", './output', " already exists. Will use as output \
directory")
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
# Name of the usecase that will be printed during simulation.
use_case = "TP-R-layered-soil-with-inner-patch-realistic-same-intrinsic"
# The name of this very file. Needed for creating log output.
thisfile = "TP-R-layered_soil_with_inner_patch-realistic-same-intrinsic.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
max_iter_num = 1000
FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = True
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: 1e-6, # h=0.1412
32: 1e-6,
64: 5e-7,
128: 5e-7
}
# 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}
# starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
timestep_size = 0.001
number_of_timesteps = 1000
# LDD scheme parameters ######################################################
Lw1 = 0.01 # /timestep_size
Lnw1 = Lw1
Lw2 = 0.01 # /timestep_size
Lnw2 = Lw2
Lw3 = 0.01 # /timestep_size
Lnw3 = 0.003
Lw4 = 0.01 # /timestep_size
Lnw4 = 0.003
Lw5 = 0.01 # /timestep_size
Lnw5 = 0.003
Lw6 = 0.01 # /timestep_size
Lnw6 = 0.003
lambda12_w = 0.5
lambda12_nw = 0.5
lambda23_w = 0.5
lambda23_nw = 0.5
lambda24_w = 0.5
lambda24_nw= 0.5
lambda25_w= 0.5
lambda25_nw= 0.5
lambda34_w = 0.5
lambda34_nw = 0.5
lambda36_w = 0.5
lambda36_nw = 0.5
lambda45_w = 0.5
lambda45_nw = 0.5
lambda46_w = 0.5
lambda46_nw = 0.5
lambda56_w = 0.5
lambda56_nw = 0.5
include_gravity = False
debugflag = False
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 = 3
# 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 #########################################################
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
# DOMAIN AND INTERFACES #######################################################
substructuring = dss.layeredSoilInnerPatch()
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 #########################################################
isRichards = {
1: True,
2: True,
3: False,
4: False,
5: False,
6: False
}
# isRichards = {
# 1: True,
# 2: True,
# 3: True,
# 4: True,
# 5: True,
# 6: True
# }
# Dict of the form: { subdom_num : viscosity }
viscosity = {
1: {'wetting' :1,
'nonwetting': 1/50},
2: {'wetting' :1,
'nonwetting': 1/50},
3: {'wetting' :1,
'nonwetting': 1/50},
4: {'wetting' :1,
'nonwetting': 1/50},
5: {'wetting' :1,
'nonwetting': 1/50},
6: {'wetting' :1,
'nonwetting': 1/50},
}
# Dict of the form: { subdom_num : density }
densities = {
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
# porosities taken from
# https://www.geotechdata.info/parameter/soil-porosity.html
# Dict of the form: { subdom_num : porosity }
porosity = {
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
L = {
1: {'wetting' :Lw1,
'nonwetting': Lnw1},
2: {'wetting' :Lw2,
'nonwetting': Lnw2},
3: {'wetting' :Lw3,
'nonwetting': Lnw3},
4: {'wetting' :Lw4,
'nonwetting': Lnw4},
5: {'wetting' :Lw5,
'nonwetting': Lnw5},
6: {'wetting' :Lw6,
'nonwetting': Lnw6}
}
# interface_num : lambda parameter for the L-scheme on that interface.
# Note that interfaces are numbered starting from 0, because
# adjacent_subdomains is a list and not a dict. Historic fuckup, I know
# We have defined above as interfaces
# # interface_vertices introduces a global numbering of interfaces.
# interface_def_points = [interface12_vertices,
# interface23_vertices,
# interface24_vertices,
# interface25_vertices,
# interface34_vertices,
# interface36_vertices,
# interface45_vertices_a,
# interface45_vertices_b,
# interface46_vertices,
# interface56_vertices,
# ]
lambda_param = {
0: {'wetting': lambda12_w,
'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
'nonwetting': lambda23_nw},
2: {'wetting': lambda24_w,
'nonwetting': lambda24_nw},
3: {'wetting': lambda25_w,
'nonwetting': lambda25_nw},
4: {'wetting': lambda34_w,
'nonwetting': lambda34_nw},
5: {'wetting': lambda36_w,
'nonwetting': lambda36_nw},
6: {'wetting': lambda45_w,
'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.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
}
# relative permeabilties
rel_perm_definition = {
1: {"wetting": "Spow2",
"nonwetting": "oneMinusSpow2"},
2: {"wetting": "Spow2",
"nonwetting": "oneMinusSpow2"},
3: {"wetting": "Spow3",
"nonwetting": "oneMinusSpow3"},
4: {"wetting": "Spow3",
"nonwetting": "oneMinusSpow3"},
5: {"wetting": "Spow3",
"nonwetting": "oneMinusSpow3"},
6: {"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": 1}},
3: {"testSpc": {"index": 2}},
4: {"testSpc": {"index": 2}},
5: {"testSpc": {"index": 2}},
6: {"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"]
###############################################################################
# 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*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.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.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.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.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**2 },
}
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
)
symbols = {"x": x,
"y": y,
"t": t}
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
exact_solution_example = hlp.generate_exact_solution_expressions(
symbols=symbols,
isRichards=isRichards,
symbolic_pressure=p_e_sym,
symbolic_capillary_pressure=pc_e_sym,
saturation_pressure_relationship=S_pc_sym,
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
intrinsic_permeability=intrinsic_permeability,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
gravity_acceleration=gravity_acceleration,
include_gravity=include_gravity,
)
source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
# 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.
# LOG FILE OUTPUT #############################################################
# read this file and print it to std out. This way the simulation can produce a
# log file with ./TP-R-layered_soil.py | tee simulation.log
f = open(thisfile, 'r')
print(f.read())
f.close()
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
# mesh_resolution and starttime are excluded, as they get passed explicitly
# to achieve parallelisation in these parameters in these parameters for
# mesh studies etc.
simulation_parameter = {
"tol": 1E-14,
"debugflag": debugflag,
"max_iter_num": max_iter_num,
"FEM_Lagrange_degree": FEM_Lagrange_degree,
"mesh_study": mesh_study,
"use_case": use_case,
"output_string": output_string,
"subdomain_def_points": subdomain_def_points,
"isRichards": isRichards,
"interface_def_points": interface_def_points,
"outer_boundary_def_points": outer_boundary_def_points,
"adjacent_subdomains": adjacent_subdomains,
# "mesh_resolution": mesh_resolution,
"viscosity": viscosity,
"porosity": porosity,
"L": L,
"lambda_param": lambda_param,
"relative_permeability": relative_permeability,
"intrinsic_permeability": intrinsic_permeability,
"sat_pressure_relationship": sat_pressure_relationship,
# "starttime": starttime,
"number_of_timesteps": number_of_timesteps,
"number_of_timesteps_to_analyse": number_of_timesteps_to_analyse,
"plot_timestep_every": plot_timestep_every,
"timestep_size": timestep_size,
"source_expression": source_expression,
"initial_condition": initial_condition,
"dirichletBC": dirichletBC,
"exact_solution": exact_solution,
"densities": densities,
"include_gravity": include_gravity,
"gravity_acceleration": gravity_acceleration,
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
for number_shift, starttime in starttimes.items():
simulation_parameter.update(
{"starttime_timestep_number_shift": number_shift}
)
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
processQueue = mp.Queue()
LDDsim = mp.Process(
target=hlp.run_simulation,
args=(
simulation_parameter,
processQueue,
starttime,
mesh_resolution
)
)
LDDsim.start()
# LDDsim.join()
# hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
# )
LDDsim.join()
if mesh_study:
simulation_output_dir = processQueue.get()
hlp.merge_spacetime_errornorms(isRichards=isRichards,
resolutions=resolutions,
output_dir=simulation_output_dir)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment