diff --git a/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py b/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py
index 74af2066898977981e19f9b5e7671834901a10a3..9c526f83cc418c5cc433876d6cebeae6a4be2adb 100755
--- a/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py
+++ b/TP-multi-patch-plus-gravity-with-same-wetting-phase-as-RR/TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py
@@ -7,30 +7,34 @@ import sympy as sym
# import domainPatch as dp
import LDDsimulation as ldd
import functools as ft
+import helpers as hlp
# import ufl as ufl
# init sympy session
sym.init_printing()
-# ----------------------------------------------------------------------------#
-# ------------------- MESH ---------------------------------------------------#
-# ----------------------------------------------------------------------------#
-mesh_resolution = 51
-# ----------------------------------------:-------------------------------------#
-# ------------------- TIME ---------------------------------------------------#
-# ----------------------------------------------------------------------------#
-timestep_size = 0.005
-number_of_timesteps = 160
-# decide how many timesteps you want analysed. Analysed means, that we write
-# out subsequent errors of the L-iteration within the timestep.
-number_of_timesteps_to_analyse = 11
+solver_tol = 5E-7
+
+############ GRID #######################ΓΌ
+mesh_resolution = 30
+timestep_size = 0.001
+number_of_timesteps = 1000
+# decide how many timesteps you want analysed. Analysed means, that we write out
+# subsequent errors of the L-iteration within the timestep.
+number_of_timesteps_to_analyse = 10
starttime = 0
-Lw = 1000
-Lnw = Lw
+Lw = 1 #/timestep_size
+Lnw=Lw
+
+l_param_w = 40
+l_param_nw = 40
+
+include_gravity = True
+debugflag = False
+analyse_condition = True
-l_param_w = 80
-l_param_nw = 80
+output_string = "./output/like_RR_number_of_timesteps{}_".format(number_of_timesteps)
# ----------------------------------------------------------------------------#
# ------------------- Domain and Interface -----------------------------------#
@@ -258,7 +262,7 @@ def rel_perm1w_prime(s):
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
- return 2*(1-s)
+ return -2*(1-s)
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
@@ -268,7 +272,7 @@ def rel_perm2w_prime(s):
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
- return 3*(1-s)**2
+ return -3*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -354,13 +358,13 @@ x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym = {
- 1: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 1: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
- 2: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x),
+ 2: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
- 3: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x),
+ 3: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t},
- 4: {'wetting': 1.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
+ 4: {'wetting': (1.0 - (1.0 + t*t)*(1.0 + x*x + y*y)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
'nonwetting': 0.0*t}
}
@@ -374,80 +378,35 @@ p_e_sym = {
pc_e_sym = dict()
for subdomain, isR in isRichards.items():
if isR:
- pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()})
+ pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
- pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - p_e_sym[subdomain]['wetting'].copy()})
+ 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.
-dtS = dict()
-div_flux = dict()
-source_expression = dict()
-exact_solution = dict()
-initial_condition = dict()
-for subdomain, isR in isRichards.items():
- dtS.update({subdomain: dict()})
- div_flux.update({subdomain: dict()})
- source_expression.update({subdomain: dict()})
- exact_solution.update({subdomain: dict()})
- initial_condition.update({subdomain: dict()})
- if isR:
- subdomain_has_phases = ["wetting"]
- else:
- subdomain_has_phases = ["wetting", "nonwetting"]
-
- # conditional for S_pc_prime
- pc = pc_e_sym[subdomain]
- dtpc = sym.diff(pc, t, 1)
- dxpc = sym.diff(pc, x, 1)
- dypc = sym.diff(pc, y, 1)
- S = sym.Piecewise((S_pc_sym[subdomain](pc), pc > 0), (1, True))
- dS = sym.Piecewise((S_pc_sym_prime[subdomain](pc), pc > 0), (0, True))
- for phase in subdomain_has_phases:
- # Turn above symbolic expression for exact solution into c code
- exact_solution[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase])}
- )
- # save the c code for initial conditions
- initial_condition[subdomain].update(
- {phase: sym.printing.ccode(p_e_sym[subdomain][phase].subs(t, 0))}
- )
- if phase == "nonwetting":
- dtS[subdomain].update(
- {phase: -porosity[subdomain]*dS*dtpc}
- )
- else:
- dtS[subdomain].update(
- {phase: porosity[subdomain]*dS*dtpc}
- )
- pa = p_e_sym[subdomain][phase]
- dxpa = sym.diff(pa, x, 1)
- dxdxpa = sym.diff(pa, x, 2)
- dypa = sym.diff(pa, y, 1)
- dydypa = sym.diff(pa, y, 2)
- mu = viscosity[subdomain][phase]
- ka = relative_permeability[subdomain][phase]
- dka = ka_prime[subdomain][phase]
- rho = densities[subdomain][phase]
- g = gravity_acceleration
-
- if phase == "nonwetting":
- # x part of div(flux) for nonwetting
- dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
- # y part of div(flux) for nonwetting
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) \
- + 1/mu*dydypa*ka(1-S)
- else:
- # x part of div(flux) for wetting
- dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
- # y part of div(flux) for wetting
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
- div_flux[subdomain].update({phase: dxdxflux + dydyflux})
- contructed_rhs = dtS[subdomain][phase] - div_flux[subdomain][phase]
- source_expression[subdomain].update(
- {phase: sym.printing.ccode(contructed_rhs)}
- )
- # print(f"source_expression[{subdomain}][{phase}] =", source_expression[subdomain][phase])
+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()
@@ -485,8 +444,8 @@ write_to_file = {
}
# initialise LDD simulation class
-simulation = ldd.LDDsimulation(tol=1E-14, debug=True, LDDsolver_tol=1E-6)
-simulation.set_parameters(output_dir="./output/",
+simulation = ldd.LDDsimulation(tol=1E-14, LDDsolver_tol=solver_tol, debug=debugflag)
+simulation.set_parameters(output_dir=output_string,
subdomain_def_points=subdomain_def_points,
isRichards=isRichards,
interface_def_points=interface_def_points,
@@ -508,12 +467,12 @@ simulation.set_parameters(output_dir="./output/",
dirichletBC_expression_strings=dirichletBC,
exact_solution=exact_solution,
densities=densities,
- include_gravity=True,
+ include_gravity=include_gravity,
write2file=write_to_file,
)
simulation.initialise()
# print(simulation.__dict__)
-simulation.run()
+simulation.run(analyse_condition=analyse_condition)
# simulation.LDDsolver(time=0, debug=True, analyse_timestep=True)
# df.info(parameters, True)
diff --git a/layered-soil-case/layered_soil.py b/layered-soil-case/layered_soil.py
new file mode 100755
index 0000000000000000000000000000000000000000..952d6e66d70c1f04df2727eb3a549fbbb6431658
--- /dev/null
+++ b/layered-soil-case/layered_soil.py
@@ -0,0 +1,134 @@
+#!/usr/bin/python3
+"""This program sets up a domain together with a decomposition into subdomains
+modelling layered soil. This is used for our LDD article with tp-tp and tp-r
+coupling.
+
+Along with the subdomains and the mesh domain markers are set upself.
+The resulting mesh is saved into files for later use.
+"""
+
+#!/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
+
+# global domain
+subdomain0_vertices = [df.Point(0.0,0.0), #
+ df.Point(13.0,0.0),#
+ df.Point(13.0,8.0),#
+ df.Point(0.0,8.0)]
+
+interface12_vertices = [df.Point(0.0, 7.0),
+ df.Point(9.0, 7.0),
+ df.Point(10.5, 7.5),
+ df.Point(12.0, 7.0),
+ df.Point(13.0, 6.5)]
+# subdomain1.
+subdomain1_vertices = [interface12_vertices[0],
+ interface12_vertices[1],
+ interface12_vertices[2],
+ interface12_vertices[3],
+ interface12_vertices[4], # southern boundary, 12 interface
+ df.Point(13.0, 8.0), # eastern boundary, outer boundary
+ df.Point(0.0, 8.0) ] # northern boundary, outer on_boundary
+
+# interface23
+interface23_vertices = [df.Point(0.0, 5.0),
+ df.Point(3.0, 5.0),
+ df.Point(6.5, 4.5),
+ df.Point(9.5, 5.0),
+ df.Point(11.5, 3.5),
+ df.Point(13.0, 3)]
+
+#subdomain1
+subdomain2_vertices = [interface23_vertices[0],
+ interface23_vertices[1],
+ interface23_vertices[2],
+ interface23_vertices[3],
+ interface23_vertices[4],
+ interface23_vertices[5], # southern boundary, 23 interface
+ subdomain1_vertices[4], # eastern boundary, outer boundary
+ subdomain1_vertices[3],
+ subdomain1_vertices[2],
+ subdomain1_vertices[1],
+ subdomain1_vertices[0] ] # northern boundary, 12 interface
+
+# interface34
+interface34_vertices = [df.Point(0.0, 2.0),
+ df.Point(4.0, 2.0),
+ df.Point(9.0, 2.5),
+ df.Point(10.5, 2.0),
+ df.Point(13.0, 1.5)]
+
+# subdomain3
+subdomain3_vertices = [interface34_vertices[0],
+ interface34_vertices[1],
+ interface34_vertices[2],
+ interface34_vertices[3],
+ interface34_vertices[4], # southern boundary, 34 interface
+ subdomain2_vertices[5], # eastern boundary, outer boundary
+ subdomain2_vertices[4],
+ subdomain2_vertices[3],
+ subdomain2_vertices[2],
+ subdomain2_vertices[1],
+ subdomain2_vertices[0] ] # northern boundary, 23 interface
+
+# subdomain4
+subdomain4_vertices = [df.Point(0.0, 0.0),
+ df.Point(13.0, 0.0), # southern boundary, outer boundary
+ subdomain3_vertices[4],# eastern boundary, outer boundary
+ subdomain3_vertices[3],
+ subdomain3_vertices[2],
+ subdomain3_vertices[1],
+ subdomain3_vertices[0] ] # northern boundary, 34 interface
+
+subdomain_vertices = [subdomain0_vertices,#
+ subdomain1_vertices,#
+ subdomain2_vertices,#
+ subdomain3_vertices,#
+ subdomain4_vertices]
+# subdomain_vertices = [subdomain0_vertices,#
+# subdomain1_vertices,#
+# subdomain2_vertices]
+
+
+# interface_vertices introduces a global numbering of interfaces.
+interface_vertices = [interface12_vertices, interface23_vertices, interface34_vertices]
+adjacent_subdomains = [[1,2], [2,3], [3,4]]
+# adjacent_subdomains = [[1,2]]
+# interface_vertices = [interface12_vertices]
+# initialise LDD simulation class
+simulation = ldd.LDDsimulation()
+simulation._init_meshes_and_markers(subdomain_vertices, mesh_resolution=2)
+# subdomain marker functions
+domain_marker = simulation.domain_marker
+mesh_subdomain = simulation.mesh_subdomain
+simulation._init_interfaces(interface_vertices, adjacent_subdomains)
+
+interface = simulation.interface
+interface_marker = simulation.interface_marker
+
+
+# Save mesh to file
+df.File('./domain_layered_soil.xml.gz') << mesh_subdomain[0]
+df.File('./global_interface_marker.pvd') << interface_marker
+#df.File('./subdomain1.xml.gz') << mesh_subdomain
+df.File('./domain_markers.pvd') << domain_marker
+#df.File('./subdomain_boundary_markers.pvd') << subdomain_boundary_marker
+
+
+
+# Save sub domains to file
+#file = File("subdomains_layered_soil.xml")
+#file << subdomains
+
+#file_double = File("subdomains_double.xml")
+#file_double << subdomains_double
+
+# Save sub domains to VTK files
+#file = File("subdomains_layered_soil.pvd")
+#file << subdomains
diff --git a/old_geometry/old_geometry.py b/old_geometry/old_geometry.py
new file mode 100644
index 0000000000000000000000000000000000000000..e19248f4dc950c0cc737acef7a1e71a93325999e
--- /dev/null
+++ b/old_geometry/old_geometry.py
@@ -0,0 +1,115 @@
+# global domain
+subdomain0_vertices = [df.Point(0.0,0.0), #
+ df.Point(13.0,0.0),#
+ df.Point(13.0,8.0),#
+ df.Point(0.0,8.0)]
+
+interface12_vertices = [df.Point(0.0, 7.0),
+ df.Point(9.0, 7.0),
+ df.Point(10.5, 7.5),
+ df.Point(12.0, 7.0),
+ df.Point(13.0, 6.5)]
+# subdomain1.
+subdomain1_vertices = [interface12_vertices[0],
+ interface12_vertices[1],
+ interface12_vertices[2],
+ interface12_vertices[3],
+ interface12_vertices[4], # southern boundary, 12 interface
+ subdomain0_vertices[2], # eastern boundary, outer boundary
+ subdomain0_vertices[3]] # northern boundary, outer on_boundary
+
+# vertex coordinates of the outer boundaries. If it can not be specified as a
+# polygon, use an entry per boundary polygon. This information is used for defining
+# the Dirichlet boundary conditions. If a domain is completely internal, the
+# dictionary entry should be 0: None
+subdomain1_outer_boundary_verts = {
+ 0: [interface12_vertices[4], #
+ subdomain0_vertices[2], # eastern boundary, outer boundary
+ subdomain0_vertices[3],
+ interface12_vertices[0]]
+}
+
+
+# interface23
+interface23_vertices = [df.Point(0.0, 5.0),
+ df.Point(3.0, 5.0),
+ # df.Point(6.5, 4.5),
+ df.Point(6.5, 5.0),
+ df.Point(9.5, 5.0),
+ # df.Point(11.5, 3.5),
+ # df.Point(13.0, 3)
+ df.Point(11.5, 5.0),
+ df.Point(13.0, 5.0)
+ ]
+
+#subdomain1
+subdomain2_vertices = [interface23_vertices[0],
+ interface23_vertices[1],
+ interface23_vertices[2],
+ interface23_vertices[3],
+ interface23_vertices[4],
+ interface23_vertices[5], # southern boundary, 23 interface
+ subdomain1_vertices[4], # eastern boundary, outer boundary
+ subdomain1_vertices[3],
+ subdomain1_vertices[2],
+ subdomain1_vertices[1],
+ subdomain1_vertices[0] ] # northern boundary, 12 interface
+
+subdomain2_outer_boundary_verts = {
+ 0: [interface23_vertices[5],
+ subdomain1_vertices[4]],
+ 1: [subdomain1_vertices[0],
+ interface23_vertices[0]]
+}
+
+
+# interface34
+interface34_vertices = [df.Point(0.0, 2.0),
+ df.Point(4.0, 2.0),
+ df.Point(9.0, 2.5),
+ df.Point(10.5, 2.0),
+ df.Point(13.0, 1.5)]
+
+# subdomain3
+subdomain3_vertices = [interface34_vertices[0],
+ interface34_vertices[1],
+ interface34_vertices[2],
+ interface34_vertices[3],
+ interface34_vertices[4], # southern boundary, 34 interface
+ subdomain2_vertices[5], # eastern boundary, outer boundary
+ subdomain2_vertices[4],
+ subdomain2_vertices[3],
+ subdomain2_vertices[2],
+ subdomain2_vertices[1],
+ subdomain2_vertices[0] ] # northern boundary, 23 interface
+
+subdomain3_outer_boundary_verts = {
+ 0: [interface34_vertices[4],
+ subdomain2_vertices[5]],
+ 1: [subdomain2_vertices[0],
+ interface34_vertices[0]]
+}
+
+# subdomain4
+subdomain4_vertices = [subdomain0_vertices[0],
+ subdomain0_vertices[1], # southern boundary, outer boundary
+ subdomain3_vertices[4],# eastern boundary, outer boundary
+ subdomain3_vertices[3],
+ subdomain3_vertices[2],
+ subdomain3_vertices[1],
+ subdomain3_vertices[0] ] # northern boundary, 34 interface
+
+subdomain4_outer_boundary_verts = {
+ 0: [subdomain4_vertices[6],
+ subdomain4_vertices[0],
+ subdomain4_vertices[1],
+ subdomain4_vertices[2]]
+}
+
+
+subdomain_def_points = [subdomain0_vertices,#
+ subdomain1_vertices,#
+ subdomain2_vertices,#
+ subdomain3_vertices,#
+ subdomain4_vertices
+ ]