Select Git revision
domainPatch.py
R-one-patch-mesh-study.py 17.04 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 = "R-one-patch-mesh-study"
# solver_tol = 5E-9
max_iter_num = 1000
FEM_Lagrange_degree = 1
mesh_study = True
# resolutions = {128: 1e-7} #[1,2,3,4,5,10,20,40,75,100]
resolutions = {
# 1: 1e-7,
2: 1e-7,
4: 1e-7,
8: 1e-7,
16: 1e-7,
32: 1e-7,
64: 1e-7,
128: 1e-7,
256: 1e-7,
512: 1e-7,
}
############ GRID #######################
# mesh_resolution = 20
timestep_size = 0.01
number_of_timesteps = 70
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 = 5
starttimes = [0.0,0.25,0.5]
# starttimes = [0.0, 0.05]
# starttimes = {
# 1: 0.0
# 2: 0.05
# 4: 0.1
# 8: 0.2
# 16: 0.4
# 32: 0.7
# 64: 1.0
# 128: 1.3
# }
Lw = 0.025 #/timestep_size
Lnw=Lw
lambda_w = 0
lambda_nw = 0
include_gravity = True
debugflag = False
analyse_condition = False
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)]
subdomain0_outer_boundary_verts = {
0: [sub_domain0_vertices[0],
sub_domain0_vertices[1],
sub_domain0_vertices[2],
sub_domain0_vertices[3],
sub_domain0_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]
# in the below list, index 0 corresponds to the 12 interface which has index 1
interface_def_points = None
# 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
0 : subdomain0_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 = None
isRichards = {
0: True, #
}
viscosity = {#
# subdom_num : viscosity
0 : {'wetting' :1,
'nonwetting': 1}, #
}
porosity = {#
# subdom_num : porosity
0: 1,#
}
# Dict of the form: { subdom_num : density }
densities = {
0: {'wetting': 1, #997,
'nonwetting': 1}, #1225}
}
gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
0: {'wetting' :Lw,
'nonwetting': Lnw},#
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
0: {'wetting' :lambda_w,
'nonwetting': lambda_nw},#
}
## relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
return s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
return (1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
subdomain1_rel_perm = {
'wetting': _rel_perm1w,#
'nonwetting': _rel_perm1nw
}
## dictionary of relative permeabilties on all domains.
relative_permeability = {#
0: subdomain1_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)
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
subdomain1_rel_perm_prime = {
'wetting': _rel_perm1w_prime,
'nonwetting': _rel_perm1nw_prime
}
# dictionary of relative permeabilties on all domains.
ka_prime = {
0: subdomain1_rel_perm_prime,
}
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)))
# def saturation(pc, index):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pc > 0, -index*pc, 1)
#
#
# def saturation_sym(pc, index):
# # inverse capillary pressure-saturation-relationship
# return -index*pc
#
#
# # 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 -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 = {
0: ft.partial(saturation_sym, index=1),
}
S_pc_sym_prime = {
0: ft.partial(saturation_sym_prime, index=1),
}
sat_pressure_relationship = {
0: 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)
epsilon_x_inner = 0.7
epsilon_x_outer = 0.99
epsilon_y_inner = epsilon_x_inner
epsilon_y_outer = epsilon_x_outer
def mollifier(x, epsilon):
""" one d mollifier """
out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1)
return out_expr
mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner)
pw_sym_x = sym.Piecewise(
(mollifier_handle(x), x**2 < epsilon_x_outer**2),
(0, True)
)
pw_sym_y = sym.Piecewise(
(mollifier_handle(y), y**2 < epsilon_y_outer**2),
(0, True)
)
def mollifier2d(x, y, epsilon):
""" one d mollifier """
out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1)
return out_expr
mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer)
pw_sym2d_x = sym.Piecewise(
(mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2),
(0, True)
)
zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise(
(mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))),
(0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))),
(mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))),
(1, True),
)
zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise(
(mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))),
(0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))),
(mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))),
(1, True),
)
zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise(
(1, y<=-2*epsilon_x_inner),
(mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))),
(0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))),
(mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))),
(1, True),
)
zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y
gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x
cutoff = gaussian/(gaussian + zero_on_shrinking)
# # construction of differentiable characteristic function.
# def smooth_characteristic_func_on_epsilon_shrinking_of_subdomain0(x, y, epsilon_x_inner, epsilon_y_inner, epsilon_x_outer, epsilon_y_outer):
# dist_to_complement_x = ft.partial(mollifier, epsilon=epsilon_x_inner)
# dist_to_complement_y = ft.partial(mollifier, epsilon=epsilon_y_inner)
# dist_to_complement = dist_to_complement_y(y)*dist_to_complement_x(x)
# dist_to_eps_shrinking_x = ft.partial(zero_outside_epsilon_thickening_of_subdomain, epsilon=epsilon_x_outer)
# dist_to_eps_shrinking_y = ft.partial(zero_outside_epsilon_thickening_of_subdomain, epsilon=epsilon_y_outer)
# dist_to_eps_shrinking = dist_to_eps_shrinking_y(y)*dist_to_eps_shrinking_x(x)
# return dist_to_complement/(dist_to_eps_shrinking + dist_to_complement)
#
# def dist_to_epsilon_thickening_of_subdomain0_complement(x, y, epsilon):
# """ calculates the (euklidian distance)^2 of a point x,y to the epsilon
# thickening of the complement of the domain.
# """
# is_inside = ((1-sym.Abs(x) > epsilon) & (1-sym.Abs(y) > epsilon))
# sym.Piecewise((0, is_inside))
p_e_sym = {
0: {'wetting': (-7 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff},
}
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]}
)
# def saturation(pressure, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
#
# sa
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)