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

fix weird git fuckug

parent 36d05db9
No related branches found
No related tags found
No related merge requests found
...@@ -20,12 +20,15 @@ mesh_resolution = 51 ...@@ -20,12 +20,15 @@ mesh_resolution = 51
# ------------------- TIME ---------------------------------------------------# # ------------------- TIME ---------------------------------------------------#
# ----------------------------------------------------------------------------# # ----------------------------------------------------------------------------#
timestep_size = 0.01 timestep_size = 0.01
number_of_timesteps = 160 number_of_timesteps = 150
# decide how many timesteps you want analysed. Analysed means, that we write # decide how many timesteps you want analysed. Analysed means, that we write
# out subsequent errors of the L-iteration within the timestep. # out subsequent errors of the L-iteration within the timestep.
number_of_timesteps_to_analyse = 11 number_of_timesteps_to_analyse = 10
starttime = 0 starttime = 0
# Lw = 0.25/timestep_size
Lw = 0.25
# ----------------------------------------------------------------------------# # ----------------------------------------------------------------------------#
# ------------------- Domain and Interface -----------------------------------# # ------------------- Domain and Interface -----------------------------------#
...@@ -192,13 +195,14 @@ porosity = { ...@@ -192,13 +195,14 @@ porosity = {
5: 1 5: 1
} }
# subdom_num : subdomain L for L-scheme # subdom_num : subdomain L for L-scheme
L = { L = {
1: {'wetting': 0.25}, 1: {'wetting': Lw},
2: {'wetting': 0.25}, 2: {'wetting': Lw},
3: {'wetting': 0.25}, 3: {'wetting': Lw},
4: {'wetting': 0.25}, 4: {'wetting': Lw},
5: {'wetting': 0.25} 5: {'wetting': Lw}
} }
lamdal_w = 32 lamdal_w = 32
...@@ -374,14 +378,88 @@ sat_pressure_relationship = { ...@@ -374,14 +378,88 @@ sat_pressure_relationship = {
x, y = sym.symbols('x[0], x[1]') # needed by UFL x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True) t = sym.symbols('t', positive=True)
# cutoff function
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)
epsilon = 0.5
# *(1-(1-y)/epsilon)**2*(1-(1-x)/epsilon)**2*(1-(-1 - y)/epsilon)**2*(1-(-1 - x)/epsilon)**2
# p_e_sym = {
# 1: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
# 2: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
# 3: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
# 4: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2},
# 5: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*((-1+sym.cos(1-y))*(-1+sym.cos(1+y))*(-1+sym.cos(1-x))*(-1+sym.cos(1+x)))**2}
# }
p_e_sym = { p_e_sym = {
1: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)}, 1: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*cutoff},
2: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x)}, 2: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*cutoff},
3: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x)}, 3: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*cutoff},
4: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x)}, 4: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x))*cutoff},
5: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)} 5: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y))*cutoff}
} }
pc_e_sym = { pc_e_sym = {
1: -1*p_e_sym[1]['wetting'], 1: -1*p_e_sym[1]['wetting'],
2: -1*p_e_sym[2]['wetting'], 2: -1*p_e_sym[2]['wetting'],
...@@ -570,7 +648,7 @@ write_to_file = { ...@@ -570,7 +648,7 @@ write_to_file = {
# initialise LDD simulation class # initialise LDD simulation class
simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-7) simulation = ldd.LDDsimulation(tol=1E-14, debug=False, LDDsolver_tol=1E-7)
simulation.set_parameters(output_dir="./output/", simulation.set_parameters(output_dir="./output/with_cutoff_function",
subdomain_def_points=subdomain_def_points, subdomain_def_points=subdomain_def_points,
isRichards=isRichards, isRichards=isRichards,
interface_def_points=interface_def_points, interface_def_points=interface_def_points,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment