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

move some of the files

parent 7a45fbb8
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py
+ 8
− 1
View file @ fd028787
...@@ -562,4 +562,11 @@ for mesh_resolution in resolutions: ...@@ -562,4 +562,11 @@ for mesh_resolution in resolutions:
simulation.initialise() simulation.initialise()
# simulation.write_exact_solution_to_xdmf() # simulation.write_exact_solution_to_xdmf()
simulation.run(analyse_condition=analyse_condition) errornorms = simulation.run(analyse_condition=analyse_condition)
for subdomain_index in errornorms.keys():
for phase, different_errornorm in errornorms[subdomain_index].items():
for errortype, errornorm in errornorms[subdomain_index][phase].items():
eoc_filename = "{}_error".format(output_string)
eocfile = open("eoc", "a")
eocfile.write( str(dx) + " " + str(err) + "\n" )
eocfile.close()
TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py deleted 100755 → 0
+ 0
− 550
View file @ 7a45fbb8
#!/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 ufl as ufl
# init sympy session
sym.init_printing()
use_case = "TP-TP-2-patch-pure-dd"
solver_tol = 6E-6
max_iter_num = 1000
############ GRID #######################
mesh_resolution = 20
timestep_size = 0.0005
number_of_timesteps = 600
# 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
starttime = 0
Lw = 1 #/timestep_size
Lnw=Lw
lambda_w = 4
lambda_nw = 4
include_gravity = False
debugflag = False
analyse_condition = True
output_string = "./output/2019-08-23-nondirichlet_number_of_timesteps{}_".format(number_of_timesteps)
##### 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]]
}
# subdomain2_outer_boundary_verts = {
# 0: [interface12_vertices[0], df.Point(0.0,0.0)],#
# 1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], #
# 2: [df.Point(1.0,0.0), interface12_vertices[1]]
# }
# subdomain2_outer_boundary_verts = {
# 0: None
# }
# 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}, #
2 : {'wetting' :1,
'nonwetting': 1}
}
porosity = {#
# subdom_num : porosity
1 : 1,#
2 : 1
}
# Dict of the form: { subdom_num : density }
densities = {
1: {'wetting': 1, #997,
'nonwetting': 1}, #1225},
2: {'wetting': 1, #997,
'nonwetting': 1}, #1225},
}
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 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
}
## relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
return s**2
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
return (1-s)**2
_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 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 2*s
def rel_perm2nw_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)
_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 saturation(pc, index):
# inverse capillary pressure-saturation-relationship
return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
def pc_sat_rel_sym(S, index):
# capillary pressure-saturation-relationship
return 1/S**(index+1) -1
pc_saturation_sym = {
1: ft.partial(pc_sat_rel_sym, index=1),
2: ft.partial(pc_sat_rel_sym, index=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=1),
# 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=1),
# 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=1),
# 3: ft.partial(saturation, index=2),
# 4: ft.partial(saturation, index=1)
}
#
# def saturation(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 saturation_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 saturation_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(saturation_sym, n_index=3, alpha=0.001),
# 2: ft.partial(saturation_sym, n_index=6, alpha=0.001),
# # 3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
# # 4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
# # 5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
# # 6: ft.partial(saturation_sym, n_index=3, alpha=0.001)
# }
#
# S_pc_sym_prime = {
# 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# 2: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
# # 3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# # 4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# # 5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# # 6: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
# }
#
# sat_pressure_relationship = {
# 1: ft.partial(saturation, n_index=3, alpha=0.001),
# 2: ft.partial(saturation, n_index=6, alpha=0.001),p1w + Spc[1]
# # 3: ft.partial(saturation, n_index=3, alpha=0.001),
# # 4: ft.partial(saturation, n_index=3, alpha=0.001),
# # 5: ft.partial(saturation, n_index=3, alpha=0.001),
# # 6: ft.partial(saturation, n_index=3, alpha=0.001)
# }
#
#############################################
# Manufacture source expressions with sympy #
#############################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
symbols = { "x": x,
"y": y,
"t": t}
# 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)
#
#
# sat_sym = {
# 1: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t),
# 2: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t)
# }
#
# Spc = {
# 1: sym.Piecewise((pc_saturation_sym[1](sat_sym[1]), sat_sym[1] > 0), (pc_saturation_sym[1](sat_sym[1]), 1>=sat_sym[1]), (0, True)),
# 2: sym.Piecewise((pc_saturation_sym[2](sat_sym[2]), sat_sym[2] > 0), (pc_saturation_sym[2](sat_sym[2]), 2>=sat_sym[2]), (0, True))
# }
#
# p1w = (-1 - (1+t*t)*(1 + x*x + y*y))#*cutoff
# p2w = p1w
# p_e_sym = {
# 1: {'wetting': p1w,
# 'nonwetting': (p1w + Spc[1])}, #*cutoff},
# 2: {'wetting': p2w,
# 'nonwetting': (p2w + Spc[2])}, #*cutoff},
# }
p_e_sym = {
1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
'nonwetting': (-1 -t*(1.1+ y + x**2))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
2: {'wetting': (-6.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': (-1 -t*(1.1 + x**2) - sym.sqrt(2+t**2)*(1.1+y)**2*y**2)}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
# 1: {'wetting': (-5 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
# 'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff},
# 2: {'wetting': (-5 - (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']})
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
write_to_file = {
'meshes_and_markers': True,
'L_iterations': True
}
# initialise LDD simulation class
simulation = ldd.LDDsimulation(
tol=1E-14,
LDDsolver_tol=solver_tol,
debug=debugflag,
max_iter_num=max_iter_num
)
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,
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()
# simulation.write_exact_solution_to_xdmf()
simulation.run(analyse_condition=analyse_condition)
TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch-realistic.py
+ 6
− 6
View file @ fd028787
...@@ -22,18 +22,18 @@ import helpers as hlp ...@@ -22,18 +22,18 @@ import helpers as hlp
sym.init_printing() sym.init_printing()
use_case = "TP-TP-layered-soil-with-inner-patch-realistic-with-gravity" use_case = "TP-TP-layered-soil-with-inner-patch-realistic-with-gravity"
solver_tol = 5E-6 solver_tol = 2E-6
############ GRID #######################ü ############ GRID #######################ü
mesh_resolution = 60 mesh_resolution = 20
timestep_size = 0.0001 timestep_size = 0.00005
number_of_timesteps = 30 number_of_timesteps = 1000
# decide how many timesteps you want analysed. Analysed means, that we write out # decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep. # subsequent errors of the L-iteration within the timestep.
number_of_timesteps_to_analyse = 10 number_of_timesteps_to_analyse = 10
starttime = 0 starttime = 0
Lw = 10 #/timestep_size Lw = 1 #/timestep_size
Lnw=Lw Lnw=Lw
lambda_w = 4 lambda_w = 4
...@@ -43,7 +43,7 @@ include_gravity = True ...@@ -43,7 +43,7 @@ include_gravity = True
debugflag = True debugflag = True
analyse_condition = False analyse_condition = False
output_string = "./output/realistic-with_gravity-post-bugfix-nondirichlet_number_of_timesteps{}_".format(number_of_timesteps) output_string = "./output/2019-08-30-{}_timesteps{}_".format(use_case, number_of_timesteps)
# global domain # global domain
subdomain0_vertices = [df.Point(-1.0,-1.0), # subdomain0_vertices = [df.Point(-1.0,-1.0), #
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment