diff --git a/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py b/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py
index a12fded8f694e1c46d51fd781f2a931f19a4b712..1a9a2c1ebb78ae2bd02e22d1ac3e8e7b1446aa0e 100755
--- a/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py
+++ b/TP-TP-2-patch-pure-dd/TP-TP-2-patch-pure-dd.py
@@ -7,11 +7,36 @@ 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 = 1E-6
+
+############ GRID #######################ü
+mesh_resolution = 30
+timestep_size = 0.0001
+number_of_timesteps = 11000
+# 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 = 10
+starttime = 0
+
+Lw = 10 #/timestep_size
+Lnw=Lw
+
+l_param_w = 20
+l_param_nw = 20
+
+include_gravity = True
+debugflag = False
+analyse_condition = False
+
+output_string = "./output/nondirichlet_number_of_timesteps{}_".format(number_of_timesteps)
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0,-1.0), #
@@ -88,17 +113,6 @@ isRichards = {
}
-solver_tol = 1E-8
-
-############ GRID #######################ü
-mesh_resolution = 31
-timestep_size = 0.0001
-number_of_timesteps = 1500
-# 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 = 11
-starttime = 0
-
viscosity = {#
# subdom_num : viscosity
1 : {'wetting' :1,
@@ -123,8 +137,7 @@ densities = {
gravity_acceleration = 9.81
-Lw = 10/timestep_size
-Lnw = Lw
+
L = {#
# subdom_num : subdomain L for L-scheme
1 : {'wetting' :Lw,
@@ -133,13 +146,13 @@ L = {#
'nonwetting': Lnw}
}
-l_param = 10
+
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :l_param,
- 'nonwetting': l_param},#
- 2 : {'wetting' :l_param,
- 'nonwetting': l_param}
+ 1 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw},#
+ 2 : {'wetting' :l_param_w,
+ 'nonwetting': l_param_nw}
}
## relative permeabilty functions on subdomain 1
@@ -189,7 +202,7 @@ def rel_perm1w_prime(s):
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
# # definition of the derivatives of the relative permeabilities
# # relative permeabilty functions on subdomain 1
@@ -199,7 +212,7 @@ def rel_perm2w_prime(s):
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -316,7 +329,7 @@ sat_pressure_relationship = {
#
# sat_pressure_relationship = {
# 1: ft.partial(saturation, n_index=3, alpha=0.001),
-# 2: ft.partial(saturation, n_index=6, 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),
@@ -331,106 +344,125 @@ sat_pressure_relationship = {
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
-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)
- }
+symbols = { "x": x,
+ "y": y,
+ "t": 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))
- }
+# 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},
+# }
-p1w = 1 - (1+t*t)*(1 + x*x + y*y)
-p2w = p1w
p_e_sym = {
- 1: {'wetting': p1w,
- 'nonwetting': p1w + Spc[1]},
- 2: {'wetting': 1 - (1+t*t)*(1 + x*x + y*y),
- 'nonwetting': p2w + Spc[2]},
-}
-
-pc_e_sym = {
- 1: p_e_sym[1]['nonwetting'] - p_e_sym[1]['wetting'],
- 2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting'],
+ 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 = {
-# 1: -1*p_e_sym[1]['wetting'],
-# 2: -1*p_e_sym[2]['wetting'],
-# }
-
-# turn above symbolic code into exact solution for dolphin and
-# construct the rhs that matches the above exact solution.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
+pc_e_sym = dict()
for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
if isR:
- subdomain_has_phases = ["wetting"]
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+ 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()
@@ -474,33 +506,34 @@ write_to_file = {
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol = 1E-14, LDDsolver_tol=solver_tol, debug = False)
-simulation.set_parameters(output_dir = "./output/",#
- 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=True,
- write2file = write_to_file,#
- )
+simulation = ldd.LDDsimulation(tol=1E-14, LDDsolver_tol=solver_tol, debug=debugflag)
+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()
+simulation.run(analyse_condition=analyse_condition)