From e2658800b3649adb65961f7b718f51771aad6d20 Mon Sep 17 00:00:00 2001
From: David <forenkram@gmx.de>
Date: Mon, 29 Jun 2020 12:13:42 +0200
Subject: [PATCH] reimplement intrinsic permeabilty and automate generation rel
 perm dicts with new module functions

---
 LDDsimulation/LDDsimulation.py                |  43 +--
 LDDsimulation/domainPatch.py                  |  24 +-
 LDDsimulation/functions.py                    | 168 +++++++++++
 LDDsimulation/helpers.py                      |  21 +-
 .../TP-R-multi-patch-with-inner-patch.py      | 260 +++---------------
 5 files changed, 255 insertions(+), 261 deletions(-)
 create mode 100644 LDDsimulation/functions.py

diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 6370b9a..ef8e8ef 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -220,7 +220,8 @@ class LDDsimulation(object):
         porosity: tp.Dict[int, tp.Dict[str, float]],#
         L: tp.Dict[int, tp.Dict[str, float]],#
         lambda_param: tp.Dict[int, tp.Dict[str, float]],#
-        relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
+        relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],
+        intrinsic_permeability: tp.Dict[int, float],
         saturation: tp.Dict[int, tp.Callable[...,None]],#
         starttime: float,#
         timestep_size: tp.Dict[int, float],#
@@ -255,6 +256,7 @@ class LDDsimulation(object):
         self.L = L
         self.lambda_param = lambda_param
         self.relative_permeability = relative_permeability
+        self.intrinsic_permeability = intrinsic_permeability
         self.saturation = saturation
         # simulation start time and time variable
         self.starttime = starttime
@@ -1388,25 +1390,26 @@ class LDDsimulation(object):
                 subdom_densities = None
             else:
                 subdom_densities = self.densities[subdom_num]
-            self.subdomain.update(#
-                    {subdom_num : dp.DomainPatch(#
-                        subdomain_index = subdom_num,#
-                        isRichards = isR,#
-                        mesh = self.mesh_subdomain[subdom_num],#
-                        viscosity = self.viscosity[subdom_num],#
-                        porosity = self.porosity[subdom_num],#
-                        outer_boundary_def_points = self.outer_boundary_def_points[subdom_num],#
-                        interfaces = self.interface,#
-                        has_interface = interface_list,#
-                        densities = subdom_densities,#
-                        include_gravity = self.include_gravity,#
-                        gravity_acceleration = self.gravity_acceleration,#
-                        L = self.L[subdom_num],#
-                        relative_permeability = self.relative_permeability[subdom_num],#
-                        saturation = self.saturation[subdom_num],#
-                        timestep_size = self.timestep_size,#
-                        interpolation_degree = self.interpolation_degree,#
-                        tol = self.tol,
+            self.subdomain.update(
+                    {subdom_num : dp.DomainPatch(
+                        subdomain_index=subdom_num,
+                        isRichards=isR,
+                        mesh=self.mesh_subdomain[subdom_num],
+                        viscosity=self.viscosity[subdom_num],
+                        porosity=self.porosity[subdom_num],
+                        outer_boundary_def_points=self.outer_boundary_def_points[subdom_num],
+                        interfaces=self.interface,
+                        has_interface=interface_list,
+                        densities=subdom_densities,
+                        include_gravity=self.include_gravity,
+                        gravity_acceleration=self.gravity_acceleration,
+                        L=self.L[subdom_num],
+                        relative_permeability=self.relative_permeability[subdom_num],
+                        intrinsic_permeability=self.intrinsic_permeability[subdom_num],
+                        saturation=self.saturation[subdom_num],
+                        timestep_size=self.timestep_size,
+                        interpolation_degree=self.interpolation_degree,
+                        tol=self.tol,
                         degree=self.FEM_Lagrange_degree)
                     }
                 )
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 56ff411..84dc252 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -118,7 +118,8 @@ class DomainPatch(df.SubDomain):
                  interfaces: tp.List[tp.Type[bi.Interface]],#
                  has_interface: tp.List[int],#
                  L: tp.Dict[str, float],#
-                 relative_permeability: tp.Dict[str, tp.Callable[...,None]],#
+                 relative_permeability: tp.Dict[str, tp.Callable[...,None]],
+                 intrinsic_permeability: float,
                  saturation: tp.Callable[..., None],#
                  timestep_size: float,#
                  densities: tp.Dict[str, float] = None,
@@ -153,6 +154,7 @@ class DomainPatch(df.SubDomain):
         self.gravity_acceleration = gravity_acceleration
         self.L = L
         self.relative_permeability = relative_permeability
+        self.intrinsic_permeability = intrinsic_permeability
         self.saturation = saturation
         # timestep size, tau in the paper
         self.timestep_size = timestep_size
@@ -346,6 +348,7 @@ class DomainPatch(df.SubDomain):
         pa_prev_iter = self.pressure_prev_iter[phase]
         phi_a  = self.testfunction[phase]
         ka = self.relative_permeability[phase]
+        ki = self.intrinsic_permeability
         S = self.saturation
         mu_a = self.viscosity[phase]
         source_a = self.source[phase]
@@ -389,16 +392,16 @@ class DomainPatch(df.SubDomain):
         # form2 and rhs2 are different for wetting and nonwetting
         # the forms of wetting and nonwetting phase differ by a minus sign
         if phase == "wetting":
-            form2 = dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
+            form2 = dt*ki/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
             rhs2 = -(porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
             if self.include_gravity:
-                rhs_gravity = (-dt/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a)))*dx
+                rhs_gravity = (-dt*ki/mu_a*df.dot(ka(S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a)))*dx
         else:
-            form2 = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
+            form2 = dt*ki/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(pa), df.grad(phi_a))*dx
             # the non wetting phase has a + before instead of a - before the bracket
             rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
             if self.include_gravity:
-                rhs_gravity = (-dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a)))*dx
+                rhs_gravity = (-dt*ki/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a)))*dx
 
         if self.has_interface is not None:
             form = form1 + form2 + sum(interface_forms)
@@ -431,6 +434,7 @@ class DomainPatch(df.SubDomain):
         """
         subdomain = self.subdomain_index
         dt = self.timestep_size
+        ki = self.intrinsic_permeability
 
         for interf_ind in self.has_interface:
             interface = self.interface[interf_ind]
@@ -474,22 +478,22 @@ class DomainPatch(df.SubDomain):
                     # the wetting phase is always present and we always need
                     # to calculate a gl0 term.
                     if self.include_gravity:
-                        flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev + z_a)
+                        flux = -ki/mu_a*ka(S(pc_prev))*df.grad(pa_prev + z_a)
                     else:
-                        flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
+                        flux = -ki/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
                 else:
                     # phase == 'nonwetting'
                     if self.isRichards: # and phase is "nonwetting"
                         # in this case the current interface has a TP neighbour
                         if self.include_gravity:
-                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(z_a)
+                            flux = -ki/mu_a*ka(1-S(pc_prev))*df.grad(z_a)
                         else:
                             flux = 0
                     else:
                         if self.include_gravity:
-                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
+                            flux = -ki/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
                         else:
-                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+                            flux = -ki/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
 
                 flux_function = df.project(flux, self.function_space["flux"][phase])
                 flux_x, flux_y = flux_function.split()
diff --git a/LDDsimulation/functions.py b/LDDsimulation/functions.py
new file mode 100644
index 0000000..8c0b568
--- /dev/null
+++ b/LDDsimulation/functions.py
@@ -0,0 +1,168 @@
+"""provide methods for creating coefficient function dicts for LDD.
+
+This module provides definitions of functions used by LDDsimulation. Functions
+are defined as normal functions and methods are provided to assemble
+dictionaries of callables needed by the LDDsimulation class.
+"""
+import typing as tp
+import dolfin as df
+import functools as ft
+
+# Functions used as relative permeabilty functions for wetting phases
+def SpowN(S,N):
+    return S**N
+
+def SpowN_prime(S,N):
+    return N*S**(N-1)
+
+# Functions used as relative permeabilty functions for nonwetting phases
+def OneMinusSpowN(S,N):
+    return (1-S)**N
+
+def OneMinusSpowN_prime(S,N):
+    return -N*(1-S)**(N-1)
+
+def generate_relative_permeability_dicts(
+        rel_perm_definition: tp.Dict[int, tp.Dict[str,str]]
+        ) -> tp.Dict[str, tp.Dict[int,tp.Dict[str, tp.Callable]]]:
+    """Generate permeabilty dictionary from definition dict.
+
+    Generate permeabilty dictionary from input definition dictionary
+    rel_perm_definition. This dictionary contains for each subdomain a
+    dictionary which in turn contains for each phase a descriptive string
+    describing which function should be used as relative permeabilty for that
+    particular phase. The supported cases are defined by the if statements
+    below.
+    The output is a dictionary containing both the relative permeability
+    function dictionaries (output["ka"]) as well as their corresponding
+    derivative dictionaries (output["ka_prime"]) for each subdomain.
+    """
+    output = dict()
+    output.update({"ka": dict()})
+    output.update({"ka_prime": dict()})
+    for subdomain, ka_dicts in rel_perm_definition.items():
+        output["ka"].update({subdomain: dict()})
+        output["ka_prime"].update({subdomain: dict()})
+        for phase, function_type in ka_dicts.items():
+            if function_type == "Spow2":
+                output["ka"][subdomain].update(
+                    {phase: ft.partial(SpowN, N=2)}
+                    )
+                output["ka_prime"][subdomain].update(
+                    {phase: ft.partial(SpowN_prime, N=2)}
+                    )
+            elif function_type == "oneMinusSpow2":
+                output["ka"][subdomain].update(
+                    {phase: ft.partial(OneMinusSpowN, N=2)}
+                    )
+                output["ka_prime"][subdomain].update(
+                    {phase: ft.partial(OneMinusSpowN_prime, N=2)}
+                    )
+            elif function_type == "Spow3":
+                output["ka"][subdomain].update(
+                    {phase: ft.partial(SpowN, N=3)}
+                    )
+                output["ka_prime"][subdomain].update(
+                    {phase: ft.partial(SpowN_prime, N=3)}
+                    )
+            elif function_type == "oneMinusSpow3":
+                output["ka"][subdomain].update(
+                    {phase: ft.partial(OneMinusSpowN, N=3)}
+                    )
+                output["ka_prime"][subdomain].update(
+                    {phase: ft.partial(OneMinusSpowN_prime, N=3)}
+                    )
+            else:
+                raise(NotImplementedError())
+
+    return output
+# class SpcRelation(object):
+#     """provide capillary pressure saturation relationships functions."""
+#     def __init__(self, Spc_on_subdomains):
+#         """build base fuction dictionary"""
+#         self._build_base_dict()
+#         self._Spc_on_subdomains = Spc_on_subdomains
+#         self._build_callables()
+#
+#     def _build_base_dict(self):
+#         """Build base dictionary."""
+#         self.__Spc = {
+#             "testSpc": self.testSpc(),
+#         }
+#
+#     def _build_callables(self):
+#         """Build callable dictionaries."""
+#         self.symbolic = dict()
+#         self.prime_symbolic = dict()
+#         self.dolfin = dict()
+#
+#         for subdomain, Spc_dict in self._Spc_on_subdomains.items():
+#             for Spc_type, parameters in Spc_dict.items():
+#                 if Spc_type == "testSpc":
+#                     self.symbolic.update(
+#                         {subdomain: ft.partialmethod(
+#                             self.__Spc[Spc_type].S_sym,
+#                             index=parameters["index"]
+#                             )},
+#                     )
+#                     self.prime_symbolic.update(
+#                         {subdomain: ft.partial(
+#                             self.__Spc[Spc_type].S_prime_sym,
+#                             index=parameters["index"]
+#                             )},
+#                     )
+#                     self.dolfin.update(
+#                         {subdomain: ft.partial(
+#                             self.__Spc[Spc_type].S,
+#                             index=parameters["index"]
+#                             )},
+#                     )
+#                 elif Spc_type == "vanGenuchten":
+#                     raise(NotImplementedError())
+#                 else:
+#                     raise(NotImplementedError())
+#
+#     class testSpc(object):
+#         """Test S-pc relationship used in R-R paper."""
+#
+#         def __init__(self):
+#             """Construct testSpc."""
+#             print("testSpc")
+#
+#         def S(pc, index):
+#             """Inverse capillary pressure-saturation-relationship.
+#
+#             Inverse capillary pressure-saturation-relationship that will
+#             be used by the simulation class
+#             this function needs to be monotonically decreasing in the
+#             capillary_pressure. Since in the richards case pc=-pw, this
+#             becomes as a function of pw a monotonically INCREASING function
+#             like in our Richards-Richards paper.
+#             However since we unify the treatment in the
+#             code for Richards and two-phase, we need the same requierment
+#             for both cases, two-phase and Richards.
+#             """
+#             # inverse capillary pressure-saturation-relationship
+#             return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
+#
+#         def S_sym(pc, index):
+#             """Inverse capillary pressure-saturation-relationship.
+#
+#             Inverse capillary pressure-saturation-relationship as symbolic
+#             expression, that will be used by
+#             helpers.generate_exact_solution_expressions()
+#             """
+#             # 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 S_prime_sym(pc, index):
+#             """Derivative of inverse pc-S-relationship.
+#
+#             Derivative of the inverse pc-S-relationship as symbolic
+#             expression, that will be used by
+#             helpers.generate_exact_solution_expressions()
+#             """
+#             # inverse capillary pressure-saturation-relationship
+#             return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
diff --git a/LDDsimulation/helpers.py b/LDDsimulation/helpers.py
index a2a9374..eff9e51 100644
--- a/LDDsimulation/helpers.py
+++ b/LDDsimulation/helpers.py
@@ -58,6 +58,7 @@ def run_simulation(
         L=parameter["L"],
         lambda_param=parameter["lambda_param"],
         relative_permeability=parameter["relative_permeability"],
+        intrinsic_permeability=parameter["intrinsic_permeability"],
         saturation=parameter["sat_pressure_relationship"],
         starttime=starttime,  # parameter["starttime"],
         number_of_timesteps=parameter["number_of_timesteps"],
@@ -142,8 +143,9 @@ def generate_exact_solution_expressions(
         isRichards: tp.Dict[int, bool],
         symbolic_pressure: tp.Dict[int, tp.Dict[str, str]],
         symbolic_capillary_pressure: tp.Dict[int, str],
-        viscosity: tp.Dict[int, tp.List[float]],#
-        porosity: tp.Dict[int, tp.Dict[str, float]],#
+        viscosity: tp.Dict[int, tp.List[float]],
+        porosity: tp.Dict[int, tp.Dict[str, float]],
+        intrinsic_permeability: tp.Dict[int, float],
         relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
         relative_permeability_prime: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],
         saturation_pressure_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
@@ -184,6 +186,7 @@ def generate_exact_solution_expressions(
 
         # conditional for S_pc_prime
         pc = symbolic_capillary_pressure[subdomain]
+        ki = intrinsic_permeability[subdomain]
         dtpc = sym.diff(pc, t, 1)
         dxpc = sym.diff(pc, x, 1)
         dypc = sym.diff(pc, y, 1)
@@ -226,20 +229,20 @@ def generate_exact_solution_expressions(
 
             if phase == "nonwetting":
                 # x part of div(flux) for nonwetting
-                div_flux_x = -1/mu*(-1*dka(1-S)*dS*dxpc*dxpa + dxdxpa*ka(1-S))
+                div_flux_x = -ki/mu*(-1*dka(1-S)*dS*dxpc*dxpa + dxdxpa*ka(1-S))
                 # y part of div(flux) for nonwetting
                 if include_gravity:
-                    div_flux_y = -1/mu*(-1*dka(1-S)*dS*dypc*(dypa + rho*g) + dydypa*ka(1-S))
+                    div_flux_y = -ki/mu*(-1*dka(1-S)*dS*dypc*(dypa + rho*g) + dydypa*ka(1-S))
                 else:
-                    div_flux_y = -1/mu*(-1*dka(1-S)*dS*dypc*dypa + dydypa*ka(1-S))
+                    div_flux_y = -ki/mu*(-1*dka(1-S)*dS*dypc*dypa + dydypa*ka(1-S))
             else:
                 # x part of div(flux) for wetting
-                div_flux_x = -1/mu*(dka(S)*dS*dxpc*dxpa + dxdxpa*ka(S))
+                div_flux_x = -ki/mu*(dka(S)*dS*dxpc*dxpa + dxdxpa*ka(S))
                 # y part of div(flux) for wetting
                 if include_gravity:
-                    div_flux_y = -1/mu*(dka(S)*dS*dypc*(dypa + rho*g) + dydypa*ka(S))
+                    div_flux_y = -ki/mu*(dka(S)*dS*dypc*(dypa + rho*g) + dydypa*ka(S))
                 else:
-                    div_flux_y = -1/mu*(dka(S)*dS*dypc*(dypa) + dydypa*ka(S))
+                    div_flux_y = -ki/mu*(dka(S)*dS*dypc*(dypa) + dydypa*ka(S))
 
             div_flux[subdomain].update({phase: div_flux_x + div_flux_y})
             contructed_rhs = dtS[subdomain][phase] + div_flux[subdomain][phase]
@@ -291,7 +294,7 @@ def generate_exact_symbolic_pc(
                 symbolic_pressure: tp.Dict[int, tp.Dict[str, str]]
             ) -> tp.Dict[str, tp.Dict[str, str]]:
     """generate symbolic capillary pressure from exact solution"""
-    
+
     symbolic_pc = dict()
     for subdomain, isR in isRichards.items():
         if isR:
diff --git a/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py b/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py
index 2d33afd..eaedb8e 100755
--- a/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py
+++ b/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py
@@ -6,6 +6,7 @@ This program sets up an LDD simulation
 import dolfin as df
 import sympy as sym
 import functools as ft
+import functions as fts
 import LDDsimulation as ldd
 import helpers as hlp
 import datetime
@@ -46,8 +47,8 @@ resolutions = {
                 # 1: 1e-6,
                 # 2: 1e-6,
                 # 4: 1e-6,
-                8: 1e-5,
-                # 16: 5e-6,
+                # 8: 5e-5,
+                16: 1e-5,
                 # 32: 5e-6,
                 # 64: 2e-6,
                 # 128: 1e-6,
@@ -62,19 +63,19 @@ timestep_size = 0.001
 number_of_timesteps = 5
 
 # LDD scheme parameters  ######################################################
-Lw1 = 0.5  # /timestep_size
+Lw1 = 0.025  # /timestep_size
 Lnw1 = Lw1
 
-Lw2 = 0.5  # /timestep_size
+Lw2 = 0.025  # /timestep_size
 Lnw2 = Lw2
 
-Lw3 = 0.5  # /timestep_size
+Lw3 = 0.025  # /timestep_size
 Lnw3 = Lw3
 
-Lw4 = 0.5  # /timestep_size
+Lw4 = 0.025  # /timestep_size
 Lnw4 = Lw4
 
-Lw5 = 0.5  # /timestep_size
+Lw5 = 0.025  # /timestep_size
 Lnw5 = Lw5
 
 
@@ -100,7 +101,7 @@ lambda15_w = 4
 lambda15_nw = 4
 
 
-include_gravity = True
+include_gravity = False
 debugflag = True
 analyse_condition = False
 
@@ -273,222 +274,35 @@ intrinsic_permeability = {
     6: 0.01,  #10e-3
 }
 
-
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
-    # relative permeabilty wetting on subdomain1
-    return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
-    # relative permeabilty nonwetting on subdomain1
-    return intrinsic_permeability[1]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
-    # relative permeabilty wetting on subdomain2
-    return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
-    # relative permeabilty nonwetting on subdomain2
-    return intrinsic_permeability[2]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 3
-def rel_perm3w(s):
-    # relative permeabilty wetting on subdomain3
-    return intrinsic_permeability[3]*s**3
-
-
-def rel_perm3nw(s):
-    # relative permeabilty nonwetting on subdomain3
-    return intrinsic_permeability[3]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 4
-def rel_perm4w(s):
-    # relative permeabilty wetting on subdomain4
-    return intrinsic_permeability[4]*s**3
-
-
-def rel_perm4nw(s):
-    # relative permeabilty nonwetting on subdomain4
-    return intrinsic_permeability[4]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 5
-def rel_perm5w(s):
-    # relative permeabilty wetting on subdomain5
-    return intrinsic_permeability[5]*s**2
-
-
-def rel_perm5nw(s):
-    # relative permeabilty nonwetting on subdomain5
-    return intrinsic_permeability[5]*(1-s)**2
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-_rel_perm3w = ft.partial(rel_perm3w)
-_rel_perm3nw = ft.partial(rel_perm3nw)
-
-_rel_perm4w = ft.partial(rel_perm4w)
-_rel_perm4nw = ft.partial(rel_perm4nw)
-
-_rel_perm5w = ft.partial(rel_perm5w)
-_rel_perm5nw = ft.partial(rel_perm5nw)
-
-
-subdomain1_rel_perm = {
-    'wetting': _rel_perm1w,
-    'nonwetting': _rel_perm1nw
-}
-
-subdomain2_rel_perm = {
-    'wetting': _rel_perm2w,
-    'nonwetting': _rel_perm2nw
-}
-
-subdomain3_rel_perm = {
-    'wetting': _rel_perm3w,
-    'nonwetting': _rel_perm3nw
-}
-
-subdomain4_rel_perm = {
-    'wetting': _rel_perm4w,
-    'nonwetting': _rel_perm4nw
-}
-
-subdomain5_rel_perm = {
-    'wetting': _rel_perm5w,
-    'nonwetting': _rel_perm5nw
-}
-
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
-    1: subdomain1_rel_perm,
-    2: subdomain2_rel_perm,
-    3: subdomain3_rel_perm,
-    4: subdomain4_rel_perm,
-    5: subdomain5_rel_perm
+# rel_perm_generator = func.relative_permeability()
+rel_perm_definition = {
+    1: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"},
+    2: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
+    3: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
+    4: {"wetting": "Spow3",
+        "nonwetting": "oneMinusSpow3"},
+    5: {"wetting": "Spow2",
+        "nonwetting": "oneMinusSpow2"}
 }
 
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
-    # relative permeabilty on subdomain1
-    return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
-    # relative permeabilty on subdomain1
-    return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-def rel_perm2w_prime(s):
-    # relative permeabilty on subdomain2
-    return intrinsic_permeability[2]*3.0*s**2
-
-
-def rel_perm2nw_prime(s):
-    # relative permeabilty on subdomain2
-    return -1*intrinsic_permeability[2]*3.0*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 3
-def rel_perm3w_prime(s):
-    # relative permeabilty on subdomain3
-    return intrinsic_permeability[3]*3.0*s**2
-
-
-def rel_perm3nw_prime(s):
-    # relative permeabilty on subdomain3
-    return -1*intrinsic_permeability[3]*3.0*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 4
-def rel_perm4w_prime(s):
-    # relative permeabilty on subdomain4
-    return intrinsic_permeability[4]*3.0*s**2
-
-
-def rel_perm4nw_prime(s):
-    # relative permeabilty on subdomain4
-    return -1*intrinsic_permeability[4]*3.0*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 5
-def rel_perm5w_prime(s):
-    # relative permeabilty on subdomain5
-    return intrinsic_permeability[5]*2.0*s
-
-
-def rel_perm5nw_prime(s):
-    # relative permeabilty on subdomain5
-    return -1*intrinsic_permeability[5]*2.0*(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)
-_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
-_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
-_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
-_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
-_rel_perm5w_prime = ft.partial(rel_perm5w_prime)
-_rel_perm5nw_prime = ft.partial(rel_perm5nw_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
-}
-
-subdomain3_rel_perm_prime = {
-    'wetting': _rel_perm3w_prime,
-    'nonwetting': _rel_perm3nw_prime
-}
-
-
-subdomain4_rel_perm_prime = {
-    'wetting': _rel_perm4w_prime,
-    'nonwetting': _rel_perm4nw_prime
-}
-
-subdomain5_rel_perm_prime = {
-    'wetting': _rel_perm5w_prime,
-    'nonwetting': _rel_perm5nw_prime
-}
-
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
-    1: subdomain1_rel_perm_prime,
-    2: subdomain2_rel_perm_prime,
-    3: subdomain3_rel_perm_prime,
-    4: subdomain4_rel_perm_prime,
-    5: subdomain5_rel_perm_prime
-}
-
-
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
+
+# Spc_on_subdomains = {
+#     1: {"testSpc": {"index": 1}},
+#     2: {"testSpc": {"index": 2}},
+#     3: {"testSpc": {"index": 2}},
+#     4: {"testSpc": {"index": 2}},
+#     5: {"testSpc": {"index": 1}}
+# }
+# Spc = fts.SpcRelation(Spc_on_subdomains)
+# S_pc_sym = Spc.symbolic
+# S_pc_sym_prime = Spc.prime_symbolic
+# sat_pressure_relationship = Spc.dolfin
 
 # this function needs to be monotonically decreasing in the capillary_pressure.
 # since in the richards case pc=-pw, this becomes as a function of pw a mono
@@ -578,6 +392,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
                         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,
@@ -634,6 +449,7 @@ if __name__ == '__main__':
         "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,
-- 
GitLab