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

set-up layered soil different intrinsic

parent d958a53d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case/TP-TP-layered_soil-different-intrinsic.py 0 → 100755
+ 394
− 0
View file @ a41a6494
#!/usr/bin/python3
""" TP-TP Layered soil simulation.
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-TP-layered_soil-realistic-different-intrinsic"
# The name of this very file. Needed for creating log output.
thisfile = "TP-TP-layered_soil-different-intrinsic.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
max_iter_num = 1000
FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False
resolutions = {
# 1: 1e-6,
# 2: 1e-6,
# 4: 1e-6,
# 8: 1e-5,
# 16: 5e-6,
32: 2e-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, 1: 0.3, 2: 0.6, 3: 0.9}
# starttimes = {0: 0.0}
timestep_size = 0.001
number_of_timesteps = 1
# LDD scheme parameters ######################################################
Lw1 = 0.007 # /timestep_size
Lnw1 = 0.005
Lw2 = 0.007 # /timestep_size
Lnw2 = 0.005
Lw3 = 0.0007 # /timestep_size
Lnw3 = 0.0005
Lw4 = 0.0007 # /timestep_size
Lnw4 = 0.0005
lambda12_w = 0.5
lambda12_nw = 0.5
lambda23_w = 0.5
lambda23_nw = 0.5
lambda34_w = 0.5
lambda34_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 = 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 = 1
# 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 INTERFACE #######################################################
substructuring = dss.layeredSoil()
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: False,
2: False,
3: False,
4: False
}
# 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},
}
densities = {
1: {'wetting': 997, # 997
'nonwetting': 1.225}, # 1.225}},
2: {'wetting': 997, # 997
'nonwetting': 1.225}, # 1.225}},
3: {'wetting': 997, # 997
'nonwetting': 1.225}, # 1.225}},
4: {'wetting': 997, # 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.37, # 0.2, # Clayey gravels, clayey sandy gravels
2: 0.022, # 0.22, # Silty gravels, silty sandy gravels
3: 0.002, # 0.37, # Clayey sands
4: 0.00022, # 0.2 # Silty or sandy clay
}
# 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}
}
# 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
lambda_param = {
0: {'wetting': lambda12_w,
'nonwetting': lambda12_nw},
1: {'wetting': lambda23_w,
'nonwetting': lambda23_nw},
2: {'wetting': lambda34_w,
'nonwetting': lambda34_nw},
}
# after Lewis, see pdf file
intrinsic_permeability = {
1: 0.1, # sand
2: 0.01, # sand, there is a range
3: 0.001, #10e-2, # clay has a range
4: 0.0001, #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"},
}
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: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
2: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
3: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
4: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
}
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_2patch = {
1: {'wetting': -7 - (1+t*t)*(1 + x*x + y*y),
'nonwetting': -1-t*(1.1 + y + x**2)**2},
2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x),
'nonwetting': -1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2)*y**2},
}
p_e_sym = {
1: {'wetting': p_e_sym_2patch[1]['wetting'],
'nonwetting': p_e_sym_2patch[1]['nonwetting']},
2: {'wetting': p_e_sym_2patch[1]['wetting'],
'nonwetting': p_e_sym_2patch[1]['nonwetting']},
3: {'wetting': p_e_sym_2patch[2]['wetting'],
'nonwetting': p_e_sym_2patch[2]['nonwetting']},
4: {'wetting': p_e_sym_2patch[2]['wetting'],
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
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% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment