diff --git a/LDDsimulation/domainSubstructuring.py b/LDDsimulation/domainSubstructuring.py
index 9a43ef4b166540a77782659d46b181dc0df786a6..6b0578dd335918876009580d546240d98ff4dec9 100644
--- a/LDDsimulation/domainSubstructuring.py
+++ b/LDDsimulation/domainSubstructuring.py
@@ -35,6 +35,67 @@ class domainSubstructuring(object):
raise(NotImplementedError())
+class globalDomain(domainSubstructuring):
+ """layered soil substructuring with inner patch."""
+
+ def __init__(self):
+ """Layered soil case with inner patch."""
+ super().__init__()
+ hlp.print_once("\n Layered Soil with inner Patch:\n")
+ # global domain
+ self.__subdomain0_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)
+ ]
+
+ self.__interface_def_points()
+ self.__adjacent_subdomains()
+ self.__subdomain_def_points()
+ self.__outer_boundary_def_points()
+
+ def __interface_def_points(self):
+ """Set self._interface_def_points."""
+
+ # interface_vertices introduces a global numbering of interfaces.
+ self.interface_def_points = None
+
+ def __adjacent_subdomains(self):
+ """Set self._adjacent_subdomains."""
+ self.adjacent_subdomains = None
+
+ def __subdomain_def_points(self):
+ """Set self._subdomain_def_points."""
+
+ self.subdomain_def_points = [
+ self.__subdomain0_vertices
+ ]
+
+ def __outer_boundary_def_points(self):
+ """Set self._outer_boundary_def_points."""
+ # 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
+ self.__subdomain0_outer_boundary_verts = {
+ 0: [self.__subdomain0_vertices[0],
+ self.__subdomain0_vertices[1],
+ self.__subdomain0_vertices[2],
+ self.__subdomain0_vertices[3],
+ self.__subdomain0_vertices[0]]
+ }
+
+
+ # if a subdomain has no outer boundary write None instead, i.e.
+ # i: None
+ # if i is the index of the inner subdomain.
+ self.outer_boundary_def_points = {
+ # subdomain number
+ 0: self.__subdomain0_outer_boundary_verts
+ }
+
class twoSoilLayers(domainSubstructuring):
"""layered soil substructuring with inner patch."""
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-alterantive.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-alterantive.py
similarity index 100%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-alterantive.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch-alterantive.py
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-linear-koefficients.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-linear-koefficients.py
similarity index 100%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-linear-koefficients.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch-linear-koefficients.py
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/R-one-patch-mesh-study-fixed-timestep.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-new-gravity-test-realistic.py
old mode 100755
new mode 100644
similarity index 79%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/R-one-patch-mesh-study-fixed-timestep.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch-new-gravity-test-realistic.py
index 14677c93cdbe98edb0217fbb0021084a48e7e232..a778552e24323c8364dcc12ba41e0c34d5367140
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/R-one-patch-mesh-study-fixed-timestep.py
+++ b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-new-gravity-test-realistic.py
@@ -19,56 +19,42 @@ datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
-use_case = "R-one-patch-mesh-study-fixed-timestep-new-errornorm"
-# solver_tol = 5E-9
-max_iter_num = 1000
+use_case = "TP-one-patch-new-gravity-test"
+# solver_tol = 6E-7
+max_iter_num = 300
FEM_Lagrange_degree = 1
-mesh_study = True
-# resolutions = {128: 1e-8} #[1,2,3,4,5,10,20,40,75,100]
+mesh_study = False
resolutions = {
- 1: 1e-8,
- 2: 1e-8,
- 4: 1e-8,
- 8: 1e-8,
- 16: 1e-8,
- 32: 1e-8,
- 64: 1e-8,
- # 128: 1e-8,
- # 256: 1e-8,
- # 512: 1e-8,
+ # 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.012
-number_of_timesteps = 1
+timestep_size = 0.001
+number_of_timesteps = 100
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 = 1
+number_of_timesteps_to_analyse = 5
starttimes = [0.0]
-# 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
+Lw = 40 #/timestep_size
Lnw=Lw
lambda_w = 0
lambda_nw = 0
-include_gravity = False
+include_gravity = True
debugflag = True
-analyse_condition = False
+analyse_condition = True
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
@@ -77,6 +63,7 @@ else:
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 = {
@@ -99,6 +86,7 @@ else:
'subsequent_errors': True
}
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0, -1.0), #
@@ -136,24 +124,24 @@ outer_boundary_def_points = {
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = None
isRichards = {
- 0: True, #
+ 0: False, #
}
viscosity = {#
# subdom_num : viscosity
0 : {'wetting' :1,
- 'nonwetting': 1}, #
+ 'nonwetting': 1/50}, #
}
porosity = {#
# subdom_num : porosity
- 0: 1,#
+ 0: 0.22,#
}
# Dict of the form: { subdom_num : density }
densities = {
- 0: {'wetting': 1, #997,
- 'nonwetting': 1}, #1225}
+ 0: {'wetting': 997,
+ 'nonwetting': 1.225}
}
gravity_acceleration = 9.81
@@ -271,64 +259,64 @@ sat_pressure_relationship = {
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)
+# 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):
@@ -349,10 +337,29 @@ cutoff = gaussian/(gaussian + zero_on_shrinking)
# sym.Piecewise((0, is_inside))
p_e_sym = {
- 0: {'wetting': (-7 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
+ 0: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff},
}
+# p_e_sym = {
+# 0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2,
+# 'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
+# }
+
+
+print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
+# # pw_sym_x*pw_sym_y
+# p_e_sym = {
+# 0: {'wetting': -3*pw_sym2d_x + 0*t,
+# 'nonwetting': -1*pw_sym_x*pw_sym_y+ 0*t},
+# }
+
+# p_e_sym = {
+# 0: {'wetting': -3*cutoff + 0*t,
+# 'nonwetting': -1*zero_on_shrinking+ 0*t},
+# }
+
+
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
@@ -361,7 +368,6 @@ for subdomain, isR in isRichards.items():
pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
- p_e_sym[subdomain]['wetting']})
-
symbols = {"x": x,
"y": y,
"t": t}
@@ -415,11 +421,11 @@ for subdomain in isRichards.keys():
)
-# 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-one-patch-new-gravity-test.py', 'r')
+print(f.read())
+f.close()
+
for starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items():
# initialise LDD simulation class
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/R-one-patch-mesh-study.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-new-gravity-test.py
similarity index 80%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/R-one-patch-mesh-study.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch-new-gravity-test.py
index ff81ca563e67101e3b0f2b6804c3e2717eaf2fda..f3408bdad19c07b4592d7fc99c7b3f8a82410537 100755
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/R-one-patch-mesh-study.py
+++ b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-new-gravity-test.py
@@ -19,56 +19,42 @@ datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
-use_case = "R-one-patch-mesh-study"
-# solver_tol = 5E-9
-max_iter_num = 1000
+use_case = "TP-one-patch-new-gravity-test"
+# solver_tol = 6E-7
+max_iter_num = 300
FEM_Lagrange_degree = 1
-mesh_study = True
-# resolutions = {128: 1e-7} #[1,2,3,4,5,10,20,40,75,100]
+mesh_study = False
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,
+ # 1: 1e-6,
+ # 2: 1e-6,
+ # 4: 1e-6,
+ # 8: 1e-6,
+ 16: 2e-6,
+ # 32: 1e-6,
+ # 64: 1e-6,
+ # 128: 1e-6,
+ # 256: 1e-6,
}
############ GRID #######################
# mesh_resolution = 20
-timestep_size = 0.01
-number_of_timesteps = 70
+timestep_size = 0.001
+number_of_timesteps = 100
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
-# }
+starttimes = [0.0]
-Lw = 0.025 #/timestep_size
+Lw = 0.04 #/timestep_size
Lnw=Lw
lambda_w = 0
lambda_nw = 0
include_gravity = True
-debugflag = False
-analyse_condition = False
+debugflag = True
+analyse_condition = True
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
@@ -77,6 +63,7 @@ else:
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 = {
@@ -99,6 +86,7 @@ else:
'subsequent_errors': True
}
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0, -1.0), #
@@ -136,7 +124,7 @@ outer_boundary_def_points = {
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = None
isRichards = {
- 0: True, #
+ 0: False, #
}
viscosity = {#
@@ -152,11 +140,11 @@ porosity = {#
# Dict of the form: { subdom_num : density }
densities = {
- 0: {'wetting': 1, #997,
- 'nonwetting': 1}, #1225}
+ 0: {'wetting': 1, #997,
+ 'nonwetting': 1} #1.225}
}
-gravity_acceleration = 9.81
+gravity_acceleration = 1
L = {#
# subdom_num : subdomain L for L-scheme
@@ -271,64 +259,64 @@ sat_pressure_relationship = {
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)
+# 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):
@@ -349,10 +337,29 @@ cutoff = gaussian/(gaussian + zero_on_shrinking)
# sym.Piecewise((0, is_inside))
p_e_sym = {
- 0: {'wetting': (-7 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
+ 0: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff},
}
+# p_e_sym = {
+# 0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2,
+# 'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
+# }
+
+
+print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
+# # pw_sym_x*pw_sym_y
+# p_e_sym = {
+# 0: {'wetting': -3*pw_sym2d_x + 0*t,
+# 'nonwetting': -1*pw_sym_x*pw_sym_y+ 0*t},
+# }
+
+# p_e_sym = {
+# 0: {'wetting': -3*cutoff + 0*t,
+# 'nonwetting': -1*zero_on_shrinking+ 0*t},
+# }
+
+
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
@@ -361,7 +368,6 @@ for subdomain, isR in isRichards.items():
pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
- p_e_sym[subdomain]['wetting']})
-
symbols = {"x": x,
"y": y,
"t": t}
@@ -415,11 +421,11 @@ for subdomain in isRichards.keys():
)
-# 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-one-patch-new-gravity-test.py', 'r')
+print(f.read())
+f.close()
+
for starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items():
# initialise LDD simulation class
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/R-one-patch-mesh-study-alternative.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-no-exact-injection.py
similarity index 51%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/R-one-patch-mesh-study-alternative.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch-no-exact-injection.py
index 825595390f4b32d71239cbf6439c8f01ea3f35ea..9f405c160fa7904d7789480e9b1b237344204e16 100755
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/R-one-patch-mesh-study-alternative.py
+++ b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-no-exact-injection.py
@@ -19,55 +19,41 @@ datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
-use_case = "R-one-patch-mesh-study"
-# solver_tol = 5E-9
-max_iter_num = 1000
+use_case = "TP-one-patch-no-exact-injection"
+# solver_tol = 6E-7
+max_iter_num = 500
FEM_Lagrange_degree = 1
-mesh_study = True
-# resolutions = {128: 1e-7} #[1,2,3,4,5,10,20,40,75,100]
+mesh_study = False
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,
+ # 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 = 10
+number_of_timesteps = 20
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.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
-# }
+starttimes = [0.0]
-Lw = 0.5 #/timestep_size
+Lw = 1 #/timestep_size
Lnw=Lw
lambda_w = 0
lambda_nw = 0
include_gravity = False
-debugflag = False
+debugflag = True
analyse_condition = True
if mesh_study:
@@ -82,11 +68,11 @@ if mesh_study:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
- 'L_iterations_per_timestep': True,
- 'solutions': True,
- 'absolute_differences': True,
+ 'L_iterations_per_timestep': False,
+ 'solutions': False,
+ 'absolute_differences': False,
'condition_numbers': analyse_condition,
- 'subsequent_errors': True
+ 'subsequent_errors': False
}
else:
write_to_file = {
@@ -99,6 +85,7 @@ else:
'subsequent_errors': True
}
+
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_vertices = [df.Point(-1.0, -1.0), #
@@ -136,7 +123,7 @@ outer_boundary_def_points = {
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = None
isRichards = {
- 0: True, #
+ 0: False, #
}
viscosity = {#
@@ -156,7 +143,7 @@ densities = {
'nonwetting': 1}, #1225}
}
-gravity_acceleration = 9.81
+gravity_acceleration = 1
L = {#
# subdom_num : subdomain L for L-scheme
@@ -169,15 +156,19 @@ lambda_param = {#
0: {'wetting' :lambda_w,
'nonwetting': lambda_nw},#
}
+intrinsic_permeability = {
+ 0: {"wetting": 1,
+ "nonwetting": 1},
+}
## relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
- return s**2
+ return intrinsic_permeability[0]["wetting"]*s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
- return (1-s)**2
+ return intrinsic_permeability[0]["nonwetting"]*(1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
@@ -186,30 +177,74 @@ 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**2
+# 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)**2
+#
+# _rel_perm2w = ft.partial(rel_perm2w)
+# _rel_perm2nw = ft.partial(rel_perm2nw)
+
+# subdomain2_rel_perm = {
+# 'wetting': _rel_perm2w,#
+# 'nonwetting': _rel_perm2nw
+# }
+#
+# subdomain2_rel_perm = {
+# 'wetting': _rel_perm1w,#
+# 'nonwetting': _rel_perm1nw
+# }
## dictionary of relative permeabilties on all domains.
relative_permeability = {#
- 0: subdomain1_rel_perm,
+ 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
+ return intrinsic_permeability[0]["wetting"]*2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return -2*(1-s)
+ return -1*intrinsic_permeability[0]["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"]*2*s
+#
+# def rel_perm2nw_prime(s):
+# # relative permeabilty on subdomain1
+# return -1*intrinsic_permeability[2]["nonwetting"]*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
+# }
+#
+# subdomain2_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,
@@ -221,6 +256,7 @@ 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)))
@@ -233,37 +269,78 @@ def saturation_sym_prime(pc, index):
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),
+ 0: 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 = {
- 0: ft.partial(saturation_sym_prime, index=1),
+ 0: 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 = {
- 0: ft.partial(saturation, index=1),
+ 0: ft.partial(saturation, index=1)
+ # 2: ft.partial(saturation, index=2),
+ # 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),
+# # 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 #
@@ -271,120 +348,132 @@ sat_pressure_relationship = {
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
+initial_condition = {
+ 0: {'wetting': sym.printing.ccode(-6*(1-x*x)*(1-y*y)), #*cutoff,
+ 'nonwetting': sym.printing.ccode(-(1-x*x)*(1-y*y))} #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+ # 2: {'wetting': sym.printing.ccode(-6*(1-x*x)*(1-y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
+ # 'nonwetting': sym.printing.ccode(-(1-x*x)*(1-y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+}
-def mollifier(x, epsilon):
+### constructing source experessions.
+injection_coord = [-0.65, -0.6]
+extraction_coord = [0.75, 0.75]
+injection_radius = 0.1
+extraction_radius = 0.1
+injection_rate = 0.05
+extraction_rate = 0.05
+# 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, height):
""" one d mollifier """
- out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1)
+ out_expr = height*sym.exp(-1/(1-(x**2 + y**2)/epsilon**2))
return out_expr
-mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner)
+mollifier2d_handle_i = ft.partial(mollifier2d, epsilon=injection_radius, height=injection_rate)
-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)
+source_in = sym.Piecewise(
+ ((0.01/(1 + 4*t**2/timestep_size))*mollifier2d_handle_i(x, y), (x-injection_coord[0])**2 + (y-injection_coord[1])**2 < injection_radius**2),
+ (0*t, 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_e = ft.partial(mollifier2d, epsilon=extraction_radius, height=extraction_rate)
-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)
+source_ext = sym.Piecewise(
+ (-0.01*(1/(1 + 4*t**2/timestep_size))*mollifier2d_handle_e(x, y), (x-extraction_coord[0])**2 + (y-extraction_coord[1])**2 < extraction_radius**2),
+ (0*t, 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),
-)
+extraction_water_ratio = 0.7
+injection_water_ratio = 0.7
-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),
-)
+# "wetting": sym.printing.ccode(extraction_water_ratio*source_ext),
+# "nonwetting": sym.printing.ccode((1-extraction_water_ratio)*source_ext)
-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),
-)
+source_expression = {
+ # 0: {"wetting": sym.printing.ccode(0*t),
+ # "nonwetting": sym.printing.ccode(0*t)},
+ 0: {"wetting": sym.printing.ccode(injection_water_ratio*(source_in + source_ext)),
+ "nonwetting": sym.printing.ccode((1-injection_water_ratio)*(source_in + source_ext))}
+}
-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)
+exact_solution = None
#
+# 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)
-# 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*(1 + x + y)), #*cutoff,
- 'nonwetting': (-1 -1*t*(1.1+y + x))}, #*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']})
+# 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']
+# # 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()
@@ -411,15 +500,20 @@ for subdomain in isRichards.keys():
# the subdomain.
for outer_boundary_ind in outer_boundary_def_points[subdomain].keys():
dirichletBC[subdomain].update(
- {outer_boundary_ind: exact_solution[subdomain]}
+ # {outer_boundary_ind: exact_solution[subdomain]}
+ {
+ outer_boundary_ind: {
+ "wetting": sym.printing.ccode(0*t),
+ "nonwetting": sym.printing.ccode(0*t)
+ }
+ }
)
-# 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-one-patch-no-exact-injection.py', 'r')
+print(f.read())
+f.close()
+
for starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items():
# initialise LDD simulation class
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/TP-one-patch-mesh-study.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-no-exact-solution.py
similarity index 73%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/TP-one-patch-mesh-study.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch-no-exact-solution.py
index bed62b609b08a817ee764588d230195c31e6a9d2..a7bb078143e3a6603348fa83a7c1f8283ff2ef05 100755
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/TP-one-patch-mesh-study.py
+++ b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-no-exact-solution.py
@@ -12,6 +12,11 @@ import datetime
import os
import pandas as pd
+#import ufl as ufl
+
+# init sympy session
+sym.init_printing()
+
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
#import ufl as ufl
@@ -19,43 +24,50 @@ datestr = date.strftime("%Y-%m-%d")
# init sympy session
sym.init_printing()
-use_case = "TP-one-patch"
-# solver_tol = 5E-9
-max_iter_num = 500
+use_case = "TP-one-patch-no-exact-solution-test"
+# solver_tol = 5E-7
+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}
+mesh_study = False
+resolutions = {
+ # 1: 1e-7, # h=2
+ # 2: 2e-5, # h=1.1180
+ # 4: 1e-6, # h=0.5590
+ # 8: 1e-6, # h=0.2814
+ # 16: 5e-7, # h=0.1412
+ 32: 1e-6,
+ # 64: 5e-7,
+ # 128: 5e-7
+ }
+
############ GRID #######################
# mesh_resolution = 20
-timestep_size = 0.01
-number_of_timesteps = 80
+timestep_size = 0.005
+number_of_timesteps = 20
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 = 4
+number_of_timesteps_to_analyse = 0
starttime = 0.0
-Lw = 0.025 #/timestep_size
+Lw = 0.25 #/timestep_size
Lnw=Lw
lambda_w = 40
lambda_nw = 40
include_gravity = False
-debugflag = False
+debugflag = True
analyse_condition = False
-output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
+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:
@@ -63,22 +75,23 @@ if mesh_study:
'space_errornorms': True,
'meshes_and_markers': True,
'L_iterations_per_timestep': False,
- 'solutions': True,
+ 'solutions': False,
'absolute_differences': False,
'condition_numbers': analyse_condition,
- 'subsequent_errors': True
+ 'subsequent_errors': False
}
else:
write_to_file = {
'space_errornorms': True,
'meshes_and_markers': True,
- 'L_iterations_per_timestep': False,
+ 'L_iterations_per_timestep': True,
'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), #
@@ -329,10 +342,29 @@ cutoff = gaussian/(gaussian + zero_on_shrinking)
# sym.Piecewise((0, is_inside))
p_e_sym = {
- 0: {'wetting': (-7 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
+ 0: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), #*cutoff,
'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff},
}
+# p_e_sym = {
+# 0: {'wetting': -(sym.cos(2*t-x - 2*y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2,
+# 'nonwetting': -6*(sym.cos(t-x -y)*sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
+# }
+
+
+print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
+# # pw_sym_x*pw_sym_y
+# p_e_sym = {
+# 0: {'wetting': -3*pw_sym2d_x + 0*t,
+# 'nonwetting': -1*pw_sym_x*pw_sym_y+ 0*t},
+# }
+
+# p_e_sym = {
+# 0: {'wetting': -3*cutoff + 0*t,
+# 'nonwetting': -1*zero_on_shrinking+ 0*t},
+# }
+
+
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
@@ -341,30 +373,90 @@ for subdomain, isR in isRichards.items():
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']
+# # 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']
+### constructing source experessions.
+injection_coord = [-0.65, -0.6]
+extraction_coord = [0.75, 0.7]
+injection_radius = 0.1
+extraction_radius = 0.075
+# 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 = 0.05*sym.exp(-1/(1-(x**2 + y**2)/epsilon**2))
+ return out_expr
+
+mollifier2d_handle_i = ft.partial(mollifier2d, epsilon=injection_radius)
+
+source_in = sym.Piecewise(
+ (-(1/(1 + t**2))*mollifier2d_handle_i(x, y), (x-injection_coord[0])**2 + (y-injection_coord[1])**2 < injection_radius**2),
+ (0*t, True)
+)
+
+mollifier2d_handle_e = ft.partial(mollifier2d, epsilon=extraction_radius)
+
+source_ext = sym.Piecewise(
+ (-0.01*(1/(1 + t**2))*mollifier2d_handle_e(x, y), (x-extraction_coord[0])**2 + (y-extraction_coord[1])**2 < extraction_radius**2),
+ (0*t, True)
+)
+
+extraction_water_ratio = 0.7
+injection_water_ratio = 0.7
+
+# "wetting": sym.printing.ccode(extraction_water_ratio*source_ext),
+# "nonwetting": sym.printing.ccode((1-extraction_water_ratio)*source_ext)
+
+source_expression = {
+ 0: {"wetting": sym.printing.ccode(0*t),
+ "nonwetting": sym.printing.ccode(0*t)},
+ # 2: {"wetting": sym.printing.ccode(injection_water_ratio*source_in),
+ # "nonwetting": sym.printing.ccode((1-injection_water_ratio)*source_in)}
+}
+exact_solution = None
+
+initial_condition = {
+ 0: {'wetting': sym.printing.ccode(-6*(1-x*x)*(1-y*y)), #*cutoff,
+ 'nonwetting': sym.printing.ccode(-(1-x*x)*(1-y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+ # 2: {'wetting': sym.printing.ccode(-6*(1-x*x)*(1-y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
+ # 'nonwetting': sym.printing.ccode(-(1-x*x)*(1-y*y))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+}
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
@@ -391,10 +483,17 @@ for subdomain in isRichards.keys():
# the subdomain.
for outer_boundary_ind in outer_boundary_def_points[subdomain].keys():
dirichletBC[subdomain].update(
- {outer_boundary_ind: exact_solution[subdomain]}
+ # {outer_boundary_ind: exact_solution[subdomain]}
+ {
+ outer_boundary_ind: {
+ "wetting": sym.printing.ccode(0*t),
+ "nonwetting": sym.printing.ccode(0*t)
+ }
+ }
)
+
# def saturation(pressure, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-purely-postive-pc.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch-purely-postive-pc.py
similarity index 100%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch-purely-postive-pc.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch-purely-postive-pc.py
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch.py b/Two-phase-Two-phase/one-patch/Archive/TP-one-patch.py
similarity index 100%
rename from Two-phase-Two-phase/one-patch/TP-one-patch/TP-one-patch.py
rename to Two-phase-Two-phase/one-patch/Archive/TP-one-patch.py
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch.py b/Two-phase-Two-phase/one-patch/TP-one-patch.py
new file mode 100755
index 0000000000000000000000000000000000000000..a4eb720790f1aee3ad5d1b170277ee7bf35744f0
--- /dev/null
+++ b/Two-phase-Two-phase/one-patch/TP-one-patch.py
@@ -0,0 +1,318 @@
+#!/usr/bin/python3
+"""TP one patch soil simulation.
+
+This program sets up an L-scheme based simulation for two-phase
+"""
+import dolfin as df
+import sympy as sym
+import functions as fts
+import LDDsimulation as ldd
+import helpers as hlp
+import datetime
+import os
+import multiprocessing as mp
+import domainSubstructuring as dss
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
+if not os.path.exists('./output'):
+ os.mkdir('./output')
+ print("Directory ", './output', " created ")
+else:
+ print("Directory ", './output', " already exists. Will use as output \
+ directory")
+
+date = datetime.datetime.now()
+datestr = date.strftime("%Y-%m-%d")
+
+# Name of the usecase that will be printed during simulation.
+use_case = "TP-one-patch-realistic-same-intrinsic"
+# The name of this very file. Needed for creating log output.
+thisfile = "TP-one-patch.py"
+
+# GENERAL SOLVER CONFIG ######################################################
+# maximal iteration per timestep
+max_iter_num = 1000
+FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
+mesh_study = False
+resolutions = {
+ # 1: 1e-6,
+ # 2: 1e-6,
+ # 4: 1e-6,
+ # 8: 1e-6,
+ # 16: 5e-6,
+ 32: 1e-6,
+ # 64: 2e-6,
+ # 128: 1e-6,
+ # 256: 1e-6,
+ }
+
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
+starttimes = {0: 0.0}
+# starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
+timestep_size = 0.001
+number_of_timesteps = 2000
+
+# LDD scheme parameters ######################################################
+Lw = 0.01 #/timestep_size
+Lnw= 0.01
+
+include_gravity = False
+debugflag = False
+analyse_condition = False
+
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 10
+
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
+if mesh_study:
+ write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
+ 'space_errornorms': True,
+ # save the mesh and marker functions to disk
+ 'meshes_and_markers': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
+ 'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
+ 'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
+ 'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
+ 'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
+ 'subsequent_errors': True
+ }
+else:
+ 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
+ }
+
+# OUTPUT FILE STRING #########################################################
+output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+
+# DOMAIN AND INTERFACE #######################################################
+substructuring = dss.globalDomain()
+interface_def_points = substructuring.interface_def_points
+adjacent_subdomains = substructuring.adjacent_subdomains
+subdomain_def_points = substructuring.subdomain_def_points
+outer_boundary_def_points = substructuring.outer_boundary_def_points
+
+# MODEL CONFIGURATION #########################################################
+isRichards = {
+ 0: False #
+ }
+
+
+viscosity = {#
+# subdom_num : viscosity
+ 0: {'wetting' :1.0,
+ 'nonwetting': 1/50} #
+}
+
+porosity = {#
+# subdom_num : porosity
+ 0: 0.2
+}
+
+# Dict of the form: { subdom_num : density }
+densities = {
+ 0: {'wetting': 997.0,
+ 'nonwetting': 1.225}
+}
+
+gravity_acceleration = 9.81
+
+L = {#
+# subdom_num : subdomain L for L-scheme
+ 0 : {'wetting' :Lw,
+ 'nonwetting': Lnw}
+}
+
+lambda_param = None
+
+intrinsic_permeability = {
+ 0: 0.01
+}
+
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 0: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"}
+}
+
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
+
+# S-pc relation
+Spc_on_subdomains = {
+ 0: {"testSpc": {"index": 1}}
+}
+
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
+
+###############################################################################
+# 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 = {
+ 0: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
+ 'nonwetting': (-1 -t*(1.0 + x**2) - sym.sin(2+t**2)**2*y**2) }
+}
+
+pc_e_sym = hlp.generate_exact_symbolic_pc(
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym
+ )
+
+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,
+ intrinsic_permeability=intrinsic_permeability,
+ 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']
+
+# BOUNDARY CONDITIONS #########################################################
+# Dictionary of dirichlet boundary conditions. If an exact solution case is
+# used, use the hlp.generate_exact_DirichletBC() method to generate the
+# Dirichlet Boundary conditions from the exact solution.
+dirichletBC = hlp.generate_exact_DirichletBC(
+ isRichards=isRichards,
+ outer_boundary_def_points=outer_boundary_def_points,
+ exact_solution=exact_solution
+ )
+# If no exact solution is provided you need to provide a dictionary of boundary
+# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see
+# the structure.
+
+# LOG FILE OUTPUT #############################################################
+# read this file and print it to std out. This way the simulation can produce a
+# log file with ./TP-R-layered_soil.py | tee simulation.log
+f = open(thisfile, 'r')
+print(f.read())
+f.close()
+
+# MAIN ########################################################################
+if __name__ == '__main__':
+ # dictionary of simualation parameters to pass to the run function.
+ # mesh_resolution and starttime are excluded, as they get passed explicitly
+ # to achieve parallelisation in these parameters in these parameters for
+ # mesh studies etc.
+ simulation_parameter = {
+ "tol": 1E-14,
+ "debugflag": debugflag,
+ "max_iter_num": max_iter_num,
+ "FEM_Lagrange_degree": FEM_Lagrange_degree,
+ "mesh_study": mesh_study,
+ "use_case": use_case,
+ "output_string": 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,
+ "intrinsic_permeability": intrinsic_permeability,
+ "sat_pressure_relationship": 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,
+ "source_expression": source_expression,
+ "initial_condition": initial_condition,
+ "dirichletBC": dirichletBC,
+ "exact_solution": exact_solution,
+ "densities": densities,
+ "include_gravity": include_gravity,
+ "gravity_acceleration": gravity_acceleration,
+ "write_to_file": write_to_file,
+ "analyse_condition": analyse_condition
+ }
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
+ for mesh_resolution, solver_tol in resolutions.items():
+ simulation_parameter.update({"solver_tol": solver_tol})
+ hlp.info(simulation_parameter["use_case"])
+ processQueue = mp.Queue()
+ LDDsim = mp.Process(
+ target=hlp.run_simulation,
+ args=(
+ simulation_parameter,
+ processQueue,
+ starttime,
+ mesh_resolution
+ )
+ )
+ LDDsim.start()
+ # LDDsim.join()
+ # hlp.run_simulation(
+ # mesh_resolution=mesh_resolution,
+ # starttime=starttime,
+ # parameter=simulation_parameter
+ # )
+
+ # LDDsim.join()
+ if mesh_study:
+ simulation_output_dir = processQueue.get()
+ hlp.merge_spacetime_errornorms(isRichards=isRichards,
+ resolutions=resolutions,
+ output_dir=simulation_output_dir)
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/debug_tests/R-one-patch-const-in-time.py b/Two-phase-Two-phase/one-patch/TP-one-patch/debug_tests/R-one-patch-const-in-time.py
deleted file mode 100755
index fb619ab2b354234d0768bbfbbf9fdaefdcc68bdf..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/debug_tests/R-one-patch-const-in-time.py
+++ /dev/null
@@ -1,522 +0,0 @@
-#!/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-const-in-time"
-# 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: 5e-7,
- # 2: 5e-7,
- # 4: 5e-7,
- # 8: 5e-7,
- # 16: 5e-7,
- # 32: 5e-7,
- 64: 5e-7,
- # 128: 5e-7,
- # 256: 1e-10,
- # 512: 1e-10,
- }
-
-############ GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.005
-number_of_timesteps = 1
-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 = 1
-starttimes = [0.0, 0.5, 1]
-
-# 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 = False
-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)))
-
-
-# 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),
-}
-
-# # 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_handle = {
-# 0: ft.partial(saturation_sym, index=1),
-# }
-#
-# S_pc_sym_prime_handle = {
-# 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': (-3 - (1+t*t)*(1 + x*x + y*y))*cutoff,
-# 'nonwetting': (-1 -t*(1+y + x**2)**2)*cutoff},
-# }
-
-p_e_sym = {
- 0: {'wetting': -6 -(1 + x*x + y*y) + 0*t}
- # 'nonwetting': -1 -(sym.sin(3*(1+y)/2*sym.pi)*sym.sin(5*(1+x)/2*sym.pi))**2},
-}
-print(f"\n\n\nsymbolic type is {type(p_e_sym[0]['wetting'])}\n\n\n")
-# # pw_sym_x*pw_sym_y
-# p_e_sym = {
-# 0: {'wetting': -3*pw_sym2d_x + 0*t,
-# 'nonwetting': -1*pw_sym_x*pw_sym_y+ 0*t},
-# }
-
-# p_e_sym = {
-# 0: {'wetting': -3*cutoff + 0*t,
-# 'nonwetting': -1*zero_on_shrinking+ 0*t},
-# }
-
-
-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']})
-
-
-
-# S_pc_sym = {
-# 0: sym.Piecewise(
-# (1, pc_e_sym[0]<= 0),
-# (S_pc_sym_handle[0](pc_e_sym[0]), ((0<pc_e_sym[0])& (pc_e_sym[0] < 1))),
-# (0, True)
-# )
-# }
-#
-# S_pc_sym_prime = {
-# 0: sym.Piecewise(
-# (S_pc_sym_prime_handle[0](pc_e_sym[0]), ((pc_e_sym[0] > 0)& (pc_e_sym[0] < 1))),
-# (0, True)
-# )
-# }
-
-
-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,
- symbolic_S_pc_relationship=S_pc_sym,
- symbolic_S_pc_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]}
- )
-
-
-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)
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/run-simulation b/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/run-simulation
deleted file mode 100755
index 0eb497502a082a0fec07a5449b1fe946d59c8cc7..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study/run-simulation
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/bin/bash
-
-[ $# -eq 0 ] && { echo "Usage: $0 simulation_file [logfile_name]"; exit 1; }
-
-SIMULATION_FILE=$1
-SIMULATION=${SIMULATION_FILE%.py}
-LOGFILE_DEFAULT="$SIMULATION.log"
-
-DATE=$(date -I)
-LOGFILE=${2:-$DATE-$LOGFILE_DEFAULT}
-
-GREETING="Simulation $SIMULATION is run on $DATE by $USER"
-
-echo $GREETING
-echo "running $SIMULATION_FILE | tee $LOGFILE"
-./$SIMULATION_FILE | tee $LOGFILE
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-constant-pressures.py b/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-constant-pressures.py
deleted file mode 100755
index 3816aa6041dafdc822e600be7ba2ee2f13e2c3dc..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-constant-pressures.py
+++ /dev/null
@@ -1,490 +0,0 @@
-#!/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-one-patch-mesh-study-fixed-timestep-constant-pressures"
-# solver_tol = 5E-9
-max_iter_num = 2000
-FEM_Lagrange_degree = 1
-mesh_study = True
-# resolutions = {128: 1e-7} #[1,2,3,4,5,10,20,40,75,100]
-resolutions = {
- 1: 1e-10,
- 2: 1e-10,
- 4: 1e-10,
- 8: 1e-10,
- 16: 1e-10,
- 32: 1e-10,
- 64: 1e-10,
- 128: 1e-10,
- 256: 1e-10,
- 512: 1e-10,
- }
-
-############ GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.01
-number_of_timesteps = 1
-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 = 1
-starttimes = [0.0, 0.05, 0.1, 0.7, 1.3]
-
-# 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.05 #/timestep_size
-Lnw=Lw
-
-lambda_w = 0
-lambda_nw = 0
-
-include_gravity = False
-debugflag = True
-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: False, #
- }
-
-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': -3 +0.0*t, #*cutoff,
- 'nonwetting': -1 +0.0*t}, #*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)
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-nonwetting0.py b/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-nonwetting0.py
deleted file mode 100755
index f15efcf437c5a960dff1b9133ba6c4f36b30f844..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-nonwetting0.py
+++ /dev/null
@@ -1,491 +0,0 @@
-#!/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-one-patch-mesh-study-fixed-timestep-nonwetting0"
-# solver_tol = 5E-9
-max_iter_num = 2000
-FEM_Lagrange_degree = 1
-mesh_study = True
-# resolutions = {128: 1e-7} #[1,2,3,4,5,10,20,40,75,100]
-resolutions = {
- 1: 1e-10,
- 2: 1e-10,
- 4: 1e-10,
- 8: 1e-10,
- 16: 1e-10,
- 32: 1e-10,
- 64: 1e-10,
- 128: 1e-10,
- 256: 1e-10,
- # 512: 1e-10,
- }
-
-############ GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.01
-number_of_timesteps = 1
-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 = 1
-# starttimes = [0.0, 0.05, 0.1, 0.7, 1.3]
-starttimes = [0.7]
-
-# 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.05 #/timestep_size
-Lnw=Lw
-
-lambda_w = 0
-lambda_nw = 0
-
-include_gravity = False
-debugflag = True
-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: False, #
- }
-
-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': 0.0*t}, #*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)
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-wetting0.py b/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-wetting0.py
deleted file mode 100755
index 9821788e1557ed8c15233d9efe9b1941cb129e7b..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep-wetting0.py
+++ /dev/null
@@ -1,490 +0,0 @@
-#!/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-one-patch-mesh-study-fixed-timestep-wetting-constantexi"
-# solver_tol = 5E-9
-max_iter_num = 2000
-FEM_Lagrange_degree = 1
-mesh_study = True
-# resolutions = {128: 1e-7} #[1,2,3,4,5,10,20,40,75,100]
-resolutions = {
- 1: 1e-10,
- 2: 1e-10,
- 4: 1e-10,
- 8: 1e-10,
- 16: 1e-10,
- 32: 1e-10,
- 64: 1e-10,
- 128: 1e-10,
- 256: 1e-10,
- 512: 1e-10,
- }
-
-############ GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.01
-number_of_timesteps = 1
-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 = 1
-starttimes = [0.0, 0.05, 0.1, 0.7, 1.3]
-
-# 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.05 #/timestep_size
-Lnw=Lw
-
-lambda_w = 0
-lambda_nw = 0
-
-include_gravity = False
-debugflag = True
-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: False, #
- }
-
-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': -10+0*t, #*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)
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep.py b/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep.py
deleted file mode 100755
index e172f97fda09fce04ade6693b1ca2562fc370e44..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/TP-one-patch-mesh-study-fixed-timestep.py
+++ /dev/null
@@ -1,490 +0,0 @@
-#!/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-one-patch-mesh-study-fixed-timestep"
-# 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: 5e-7,
- 2: 5e-7,
- 4: 5e-7,
- 8: 5e-7,
- 16: 5e-7,
- 32: 5e-7,
- 64: 5e-7,
- 128: 5e-7,
- 256: 5e-7,
- # 512: 1e-10,
- }
-
-############ GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.0025
-number_of_timesteps = 1
-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 = 1
-starttimes = [0.0, 0.05, 0.1, 0.7, 1.3]
-
-# 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 = False
-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: False, #
- }
-
-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)
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/run-simulation b/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/run-simulation
deleted file mode 100755
index 0eb497502a082a0fec07a5449b1fe946d59c8cc7..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/mesh_study_for_fixed_timestep/run-simulation
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/bin/bash
-
-[ $# -eq 0 ] && { echo "Usage: $0 simulation_file [logfile_name]"; exit 1; }
-
-SIMULATION_FILE=$1
-SIMULATION=${SIMULATION_FILE%.py}
-LOGFILE_DEFAULT="$SIMULATION.log"
-
-DATE=$(date -I)
-LOGFILE=${2:-$DATE-$LOGFILE_DEFAULT}
-
-GREETING="Simulation $SIMULATION is run on $DATE by $USER"
-
-echo $GREETING
-echo "running $SIMULATION_FILE | tee $LOGFILE"
-./$SIMULATION_FILE | tee $LOGFILE
diff --git a/Two-phase-Two-phase/one-patch/TP-one-patch/run-simulation b/Two-phase-Two-phase/one-patch/TP-one-patch/run-simulation
deleted file mode 100755
index 0eb497502a082a0fec07a5449b1fe946d59c8cc7..0000000000000000000000000000000000000000
--- a/Two-phase-Two-phase/one-patch/TP-one-patch/run-simulation
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/bin/bash
-
-[ $# -eq 0 ] && { echo "Usage: $0 simulation_file [logfile_name]"; exit 1; }
-
-SIMULATION_FILE=$1
-SIMULATION=${SIMULATION_FILE%.py}
-LOGFILE_DEFAULT="$SIMULATION.log"
-
-DATE=$(date -I)
-LOGFILE=${2:-$DATE-$LOGFILE_DEFAULT}
-
-GREETING="Simulation $SIMULATION is run on $DATE by $USER"
-
-echo $GREETING
-echo "running $SIMULATION_FILE | tee $LOGFILE"
-./$SIMULATION_FILE | tee $LOGFILE