diff --git a/LDDsimulation.py b/LDDsimulation.py
index 4f2e9964893588e66349ca2c4bfce27791687cc5..06aeb1a76bfae38120ee6ac78c44bf096b4105a7 100644
--- a/LDDsimulation.py
+++ b/LDDsimulation.py
@@ -98,7 +98,11 @@ class LDDsimulation(object):
porosity: tp.Dict[int, float],#
L: tp.Dict[int, tp.List[float]],#
lambda_param: tp.Dict[int, tp.List[float]],#
- relative_permeability: tp.Dict[int, tp.Callable[...,None]])-> None:
+ relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
+ saturation: tp.Dict[int, tp.Callable[...,None]],#
+ timestep: tp.Dict[int, float],#
+ )-> None:
+
""" set parameters of an instance of class LDDsimulation"""
self.output_dir = output_dir
self.subdomain_def_points = subdomain_def_points
@@ -111,6 +115,8 @@ class LDDsimulation(object):
self.L = L
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
+ self.saturation = saturation
+ self.timestep = timestep
self._parameters_set = True
def initialise(self) -> None:
@@ -261,7 +267,9 @@ class LDDsimulation(object):
has_interface = interface_list,#
L = self.L[subdom_num],#
lambda_param = self.lambda_param[subdom_num],#
- relative_permeability = self.relative_permeability[subdom_num]#
+ relative_permeability = self.relative_permeability[subdom_num],#
+ saturation = self.saturation[subdom_num],#
+ timestep = self.timestep[subdom_num],#
))
diff --git a/RR-2-patch-test.py b/RR-2-patch-test.py
index 509d1a29f54e55a81a47a26d0d401875bbf49c2c..9cc714754e7500cd7ad5478dd79e58213ff0f92e 100755
--- a/RR-2-patch-test.py
+++ b/RR-2-patch-test.py
@@ -43,6 +43,17 @@ interface_def_points = [interface12_vertices]
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]]
is_Richards = [True, True]
+
+Tmax = 2
+dt1 = 0.1
+dt2 = dt1
+# the timestep is taken as a dictionary to be able to run different timesteps on
+# different subdomains.
+timestep = {
+ 1: dt1,
+ 2: dt2
+}
+
viscosity = {#
# subdom_num : viscosity
1 : [1], #
@@ -67,22 +78,44 @@ lambda_param = {#
2 : [4]
}
+## relative permeabilty functions on subdomain 1
def rel_perm1(s):
# relative permeabilty on subdomain1
return s**2
_rel_perm1 = ft.partial(rel_perm1)
+subdomain1_rel_perm = {
+ 'wetting': _rel_perm1,#
+ 'nonwetting': None
+}
+## relative permeabilty functions on subdomain 2
def rel_perm2(s):
# relative permeabilty on subdomain2
return s**3
_rel_perm2 = ft.partial(rel_perm2)
+subdomain2_rel_perm = {
+ 'wetting': _rel_perm2,#
+ 'nonwetting': None
+}
+
+## dictionary of relative permeabilties on all domains.
relative_permeability = {#
- 1: _rel_perm1,
- 2: _rel_perm2,
+ 1: subdomain1_rel_perm,
+ 2: subdomain2_rel_perm
}
+def saturation(pressure, subdomain_index):
+ # inverse capillary pressure-saturation-relationship
+ return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
+
+sat_pressure_relationship = {#
+ 1: ft.partial(saturation, subdomain_index = 1),#
+ 2: ft.partial(saturation, subdomain_index = 2)
+}
+
+
# def saturation(pressure, subdomain_index):
# # inverse capillary pressure-saturation-relationship
@@ -103,7 +136,10 @@ simulation.set_parameters(output_dir = "",#
porosity = porosity,#
L = L,#
lambda_param = lambda_param,#
- relative_permeability = relative_permeability)
+ relative_permeability = relative_permeability,#
+ saturation = sat_pressure_relationship,#
+ timestep = timestep#
+ )
simulation.initialise()
# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
diff --git a/domainPatch.py b/domainPatch.py
index fe76f08eecb0a4c5f0f700ad12420f98262034d0..a878315e64992eca9ee31ed6e1b8f19c4cd8debf 100644
--- a/domainPatch.py
+++ b/domainPatch.py
@@ -111,7 +111,9 @@ class DomainPatch(df.SubDomain):
has_interface: tp.List[int],#
L: tp.List[float],#
lambda_param: tp.List[float],#
- relative_permeability: tp.Callable[...,None]#
+ relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
+ saturation: tp.Callable[..., None],#
+ timestep: float,#
):
# because we declare our own __init__ method for the derived class BoundaryPart,
# we overwrite the __init__-method of the parent class df.SubDomain. However,
@@ -128,14 +130,36 @@ class DomainPatch(df.SubDomain):
self.L = L
self.lambda_param = lambda_param
self.relative_permeability = relative_permeability
-
+ self.saturation = saturation
+ self.timestep = timestep
+ # variable holding the previous iteration.
+ self.previous_iteration = None
+ # measures are set by self._init_measures()
+ self.dx = None
+ self.ds = None
self._init_function_space()
self._init_dof_and_vertex_maps()
self._init_boundary_markers()
+ self._init_measures()
# END constructor
#### PUBLIC METHODS
+ # def governing_form(self, phase: str = 'wetting') -> ufl.object:
+ # """ return the governing form of the subdomain"""
+ # if self.isRichards:
+ # # define measures
+ #
+ #
+ # else:
+ # if phase == "wetting":
+ # print("governing form wetting phase")
+ # elif phase == "nonwetting":
+ # print("governing form nonwetting phase")
+ # else:
+ # raise RuntimeError('missmatch for input parameter phase. \nExpected either phase == wetting or phase == nonwetting ')
+ #
+ # return form
#### PRIVATE METHODS
def _init_function_space(self) -> None:
@@ -193,6 +217,18 @@ class DomainPatch(df.SubDomain):
# each interface gets marked with the global interface index.
self.interface[glob_index].mark(self.boundary_marker, glob_index)
+ def _init_measures(self):
+ """ define measures for the governing form
+
+ This method sets the class Variables
+ self.dx
+ self.ds
+ """
+ # domain measure
+ self.dx = df.Measure('dx', domain = self.mesh)
+ # measure over the interfaces
+ self.ds = df.Measure('ds', domain = self.mesh, subdomain_data = self.boundary_marker)
+
# END is_Richards
# END OF CLASS DomainPatch
diff --git a/layered_soil.py b/layered_soil.py
index 0848f79a4db688b2caf74d8507ac5c2970fd8b69..952d6e66d70c1f04df2727eb3a549fbbb6431658 100755
--- a/layered_soil.py
+++ b/layered_soil.py
@@ -103,11 +103,11 @@ adjacent_subdomains = [[1,2], [2,3], [3,4]]
# interface_vertices = [interface12_vertices]
# initialise LDD simulation class
simulation = ldd.LDDsimulation()
-simulation.init_meshes_and_markers(subdomain_vertices, mesh_resolution=2)
+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)
+simulation._init_interfaces(interface_vertices, adjacent_subdomains)
interface = simulation.interface
interface_marker = simulation.interface_marker