domainPatch.py

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],#
                     ))
 
 

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)

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

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