Select Git revision
TP-one-patch.py
TP-one-patch.py 12.27 KiB
#!/usr/bin/python3
"""TP one patch soil simulation, Copyright 2020, David Seus
This program runs an L simulation on one domain without domain decomposition.
# LICENCE #####################################################################
Copyright 2020, David Seus
david.seus[at]ians.uni-stuttgart.de
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
###############################################################################
"""
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-one-patch-realistic-same-intrinsic"
# The name of this very file. Needed for creating log output.
thisfile = "TP-one-patch.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-6,
# 16: 5e-6,
32: 1e-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}
# starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
timestep_size = 0.001
number_of_timesteps = 2000
# LDD scheme parameters ######################################################
Lw = 0.01 #/timestep_size
Lnw= 0.01
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 = 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 = 10
# 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': True,
'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.globalDomain()
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 = {
0: False #
}
viscosity = {#
# subdom_num : viscosity
0: {'wetting' :1.0,
'nonwetting': 1/50} #
}
porosity = {#
# subdom_num : porosity
0: 0.2
}
# Dict of the form: { subdom_num : density }
densities = {
0: {'wetting': 997.0,
'nonwetting': 1.225}
}
gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
0 : {'wetting' :Lw,
'nonwetting': Lnw}
}
lambda_param = None
intrinsic_permeability = {
0: 0.01
}
# RELATIVE PEMRMEABILITIES
rel_perm_definition = {
0: {"wetting": "Spow2",
"nonwetting": "oneMinusSpow2"}
}
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 = {
0: {"testSpc": {"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 #
###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym = {
0: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
'nonwetting': (-1 -t*(1.0 + x**2) - sym.sin(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)