Select Git revision
run_experiments.jl
TP-TP-2-patch-mesh-study.py 18.57 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")
#import ufl as ufl
# init sympy session
sym.init_printing()
use_case = "TP-TP-2-patch-realistic"
# solver_tol = 6E-7
max_iter_num = 300
FEM_Lagrange_degree = 1
mesh_study = False
resolutions = {
# 1: 1e-6,
# 2: 1e-6,
# 4: 1e-6,
# 8: 1e-6,
# 16: 1e-6,
32: 1e-6,
# 64: 1e-6,
# 128: 1e-6,
# 256: 1e-6,
}
############ GRID #######################
# mesh_resolution = 20
timestep_size = 0.001
number_of_timesteps = 20
plot_timestep_every = 2
# 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 = 6
starttimes = [0.0]
Lw = 0.025 #/timestep_size
Lnw=Lw
lambda_w = 4
lambda_nw = 4
include_gravity = True
debugflag = True
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': True,
'solutions': True,
'absolute_differences': True,
'condition_numbers': analyse_condition,
'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
}
##### 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)]
# interface between subdomain1 and subdomain2
interface12_vertices = [df.Point(-1.0, 0.0),
df.Point(1.0, 0.0) ]
# subdomain1.
sub_domain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
sub_domain0_vertices[2],
sub_domain0_vertices[3]]
# 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 internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
0: [interface12_vertices[1], #
sub_domain0_vertices[2],
sub_domain0_vertices[3], #
interface12_vertices[0]]
}
# subdomain2
sub_domain2_vertices = [sub_domain0_vertices[0],
sub_domain0_vertices[1],
interface12_vertices[1],
interface12_vertices[0] ]
subdomain2_outer_boundary_verts = {
0: [interface12_vertices[0], #
sub_domain0_vertices[0],
sub_domain0_vertices[1],
interface12_vertices[1]]
}
# 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]
# in the below list, index 0 corresponds to the 12 interface which has index 1
interface_def_points = [interface12_vertices]
# 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
}
# 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]]
isRichards = {
1: False, #
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.0022
}
# Dict of the form: { subdom_num : density }
densities = {
1: {'wetting': 997,
'nonwetting': 1.225},
2: {'wetting': 997,
'nonwetting': 1.225},
}
intrinsic_permeability = {
1: {"wetting": 1,
"nonwetting": 1},
2: {"wetting": 1,
"nonwetting": 1},
}
gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
1 : {'wetting' :Lw,
'nonwetting': Lnw},#
2 : {'wetting' :Lw,
'nonwetting': Lnw}
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :lambda_w,
'nonwetting': lambda_nw},#
2 : {'wetting' :lambda_w,
'nonwetting': lambda_nw}
}
## relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
return intrinsic_permeability[1]["wetting"]*s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
return intrinsic_permeability[1]["nonwetting"]*(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]["wetting"]*s**3
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
return intrinsic_permeability[2]["nonwetting"]*(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
}
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
def rel_perm1w_prime(s):
# relative permeabilty on subdomain1
return intrinsic_permeability[1]["wetting"]*2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
return -intrinsic_permeability[1]["nonwetting"]*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 intrinsic_permeability[2]["wetting"]*3*s**2
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
return -intrinsic_permeability[2]["nonwetting"]*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
}
# 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)
}
#############################################
# 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': (-6 - (1+t*t)*(1 + x*x + y*y)),
'nonwetting': -1 -t*(1.1 + y*y)}, #*(1-x)**2*(1+x)**2*(1-y)**2},
2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**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)
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
else:
pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy()
- p_e_sym[subdomain]['wetting'].copy()})
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]}
)
# def saturation(pressure, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
#
# sa
f = open('TP-TP-2-patch-mesh-study.py', 'r')
print(f.read())
f.close()
for starttime in starttimes:
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)