Select Git revision
TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py 18.77 KiB
#!/usr/bin/python3
import dolfin as df
# import mshr
# import numpy as np
import sympy as sym
# import typing as tp
# import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
import helpers as hlp
import datetime
import os
import pandas as pd
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
# solver_tol = 6E-7
use_case = "TP-TP-4-patch-nonwetting-zero"
max_iter_num = 1000
FEM_Lagrange_degree = 1
mesh_study = False
resolutions = {
# 1: 1e-7, # h=2
# 2: 2e-5, # h=1.1180
# 4: 1e-6, # h=0.5590
# 8: 1e-6, # h=0.2814
# 16: 5e-7, # h=0.1412
32: 7e-7, # h=0.0706
# 64: 1e-6,
# 128: 5e-7
}
############ GRID #######################
# mesh_resolution = 20
timestep_size = 0.005
number_of_timesteps = 250
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 = 2
starttime = 0.0
Lw = 0.05 #/timestep_size
Lnw=Lw
lambda_w = 20
lambda_nw = 20
include_gravity = True
debugflag = False
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)
# toggle what should be written to files
if mesh_study:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
'L_iterations_per_timestep': False,
'solutions': False,
'absolute_differences': False,
'condition_numbers': analyse_condition,
'subsequent_errors': False
}
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
}
# ----------------------------------------------------------------------------#
# ------------------- Domain and Interface -----------------------------------#
# ----------------------------------------------------------------------------#
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0, -1.0),
df.Point(1.0, -1.0),
df.Point(1.0, 1.0),
df.Point(-1.0, 1.0)]
# interfaces
interface12_vertices = [df.Point(0.0, 0.0),
df.Point(1.0, 0.0)]
interface14_vertices = [df.Point(0.0, 0.0),
df.Point(0.0, 1.0)]
interface23_vertices = [df.Point(0.0, 0.0),
df.Point(0.0, -1.0)]
interface34_vertices = [df.Point(-1.0, 0.0),
df.Point(0.0, 0.0)]
# subdomain1.
sub_domain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
sub_domain0_vertices[2],
df.Point(0.0, 1.0)]
# 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 the Dirichlet boundary conditions. If a domain is completely inter-
# nal, the dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
0: [interface12_vertices[1],
sub_domain0_vertices[2],
df.Point(0.0, 1.0)]
}
# subdomain2
sub_domain2_vertices = [interface23_vertices[1],
sub_domain0_vertices[1],
interface12_vertices[1],
interface12_vertices[0]]
subdomain2_outer_boundary_verts = {
0: [df.Point(0.0, -1.0),
sub_domain0_vertices[1],
interface12_vertices[1]]
}
sub_domain3_vertices = [interface34_vertices[0],
sub_domain0_vertices[0],
interface23_vertices[1],
interface23_vertices[0]]
subdomain3_outer_boundary_verts = {
0: [interface34_vertices[0],
sub_domain0_vertices[0],
interface23_vertices[1]]
}
sub_domain4_vertices = [interface34_vertices[0],
interface34_vertices[1],
interface14_vertices[1],
sub_domain0_vertices[3]]
subdomain4_outer_boundary_verts = {
0: [interface14_vertices[1],
sub_domain0_vertices[3],
interface34_vertices[0]]
}
# list of subdomains given by the boundary polygon vertices.
# Subdomains are given as a list of dolfin points forming
# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
# to create the subdomain. subdomain_def_points[0] contains the
# vertices of the global simulation domain and subdomain_def_points[i] contains
# the vertices of the subdomain i.
subdomain_def_points = [sub_domain0_vertices,
sub_domain1_vertices,
sub_domain2_vertices,
sub_domain3_vertices,
sub_domain4_vertices]
# in the below list, index 0 corresponds to the 12 interface which has global
# marker value 1
interface_def_points = [interface12_vertices,
interface14_vertices,
interface23_vertices,
interface34_vertices]
# adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1, 2], [1, 4], [2, 3], [3, 4]]
# if a subdomain has no outer boundary write None instead, i.e.
# i: None
# if i is the index of the inner subdomain.
outer_boundary_def_points = {
# subdomain number
1: subdomain1_outer_boundary_verts,
2: subdomain2_outer_boundary_verts,
3: subdomain3_outer_boundary_verts,
4: subdomain4_outer_boundary_verts
}
isRichards = {
1: False,
2: False,
3: False,
4: False
}
viscosity = {
1: {'wetting' :1,
'nonwetting': 1},
2: {'wetting' :1,
'nonwetting': 1},
3: {'wetting' :1,
'nonwetting': 1},
4: {'wetting' :1,
'nonwetting': 1},
}
# Dict of the form: { subdom_num : density }
densities = {
1: {'wetting': 1, #997
'nonwetting':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}}
}
gravity_acceleration = 9.81
# porosities taken from
# 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
}
# subdom_num : subdomain L for L-scheme
L = {
1: {'wetting' :Lw,
'nonwetting': Lnw},
2: {'wetting' :Lw,
'nonwetting': Lnw},
3: {'wetting' :Lw,
'nonwetting': Lnw},
4: {'wetting' :Lw,
'nonwetting': Lnw}
}
# subdom_num : lambda parameter for the L-scheme
lambda_param = {
1: {'wetting': lambda_w,
'nonwetting': lambda_nw},#
2: {'wetting': lambda_w,
'nonwetting': lambda_nw},#
3: {'wetting': lambda_w,
'nonwetting': lambda_nw},#
4: {'wetting': lambda_w,
'nonwetting': lambda_nw},#
}
# relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty on subdomain1
return s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
return (1-s)**2
## relative permeabilty functions on subdomain 2
# relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty on subdomain2
return s**3
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
return (1-s)**3
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)
subdomain1_rel_perm = {
'wetting': _rel_perm1w,#
'nonwetting': _rel_perm1nw
}
subdomain2_rel_perm = {
'wetting': _rel_perm2w,#
'nonwetting': _rel_perm2nw
}
subdomain3_rel_perm = subdomain2_rel_perm.copy()
subdomain4_rel_perm = subdomain1_rel_perm.copy()
# dictionary of relative permeabilties on all domains.
relative_permeability = {
1: subdomain1_rel_perm,
2: subdomain2_rel_perm,
3: subdomain3_rel_perm,
4: subdomain4_rel_perm
}
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
def rel_perm1w_prime(s):
# relative permeabilty on subdomain1
return 2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
return -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 subdomain1
return 3*s**2
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
return -3*(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
}
# _rel_perm3_prime = ft.partial(rel_perm2_prime)
subdomain3_rel_perm_prime = subdomain2_rel_perm_prime.copy()
# _rel_perm4_prime = ft.partial(rel_perm1_prime)
subdomain4_rel_perm_prime = subdomain1_rel_perm_prime.copy()
# dictionary of relative permeabilties on all domains.
ka_prime = {
1: subdomain1_rel_perm_prime,
2: subdomain2_rel_perm_prime,
3: subdomain3_rel_perm_prime,
4: subdomain4_rel_perm_prime
}
# this function needs to be monotonically decreasing in the capillary_pressure.
# since in the richards case pc=-pw, this becomes as a function of pw a mono
# tonically INCREASING function like in our Richards-Richards paper. However
# since we unify the treatment in the code for Richards and two-phase, we need
# the same requierment
# for both cases, two-phase and Richards.
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)
}
#############################################
# 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': (-3.0 - (1.0 + t*t)*(1.0 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
2: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
3: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
4: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t}
}
# pc_e_sym = {
# 1: -1*p_e_sym[1]['wetting'],
# 2: -1*p_e_sym[2]['wetting'],
# 3: -1*p_e_sym[3]['wetting'],
# 4: -1*p_e_sym[4]['wetting']
# }
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
- p_e_sym[subdomain]['wetting']})
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,
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']
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
# None should be written instead of an expression.
# This is a bit of a brainfuck:
# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
# Since a domain patch can have several disjoint outer boundary parts, the
# expressions need to get an enumaration index which starts at 0.
# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of
# subdomain ind and boundary part j.
# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting']
# return the actual expression needed for the dirichlet condition for both
# phases if present.
# 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
if outer_boundary_def_points[subdomain] is None:
dirichletBC.update({subdomain: None})
else:
dirichletBC.update({subdomain: dict()})
# set the dirichlet conditions to be the same code as exact solution on
# the subdomain.
for outer_boundary_ind in outer_boundary_def_points[subdomain].keys():
dirichletBC[subdomain].update(
{outer_boundary_ind: exact_solution[subdomain]}
)
for mesh_resolution, solver_tol in resolutions.items():
# initialise LDD simulation class
simulation = ldd.LDDsimulation(
tol=1E-14,
LDDsolver_tol=solver_tol,
debug=debugflag,
max_iter_num=max_iter_num,
FEM_Lagrange_degree=FEM_Lagrange_degree,
mesh_study=mesh_study
)
simulation.set_parameters(use_case=use_case,
output_dir=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,
saturation=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,
sources=source_expression,
initial_conditions=initial_condition,
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
include_gravity=include_gravity,
write2file=write_to_file,
)
simulation.initialise()
output_dir = simulation.output_dir
# simulation.write_exact_solution_to_xdmf()
output = simulation.run(analyse_condition=analyse_condition)
for subdomain_index, subdomain_output in output.items():
mesh_h = subdomain_output['mesh_size']
for phase, different_errornorms in subdomain_output['errornorm'].items():
filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
# for errortype, errornorm in different_errornorms.items():
# eocfile = open("eoc_filename", "a")
# eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
# eocfile.close()
# if subdomain.isRichards:mesh_h
data_dict = {
'mesh_parameter': mesh_resolution,
'mesh_h': mesh_h,
}
for error_type, errornorms in different_errornorms.items():
data_dict.update(
{error_type: errornorms}
)
errors = pd.DataFrame(data_dict, index=[mesh_resolution])
# check if file exists
if os.path.isfile(filename) == True:
with open(filename, 'a') as f:
errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
else:
errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)