diff --git a/LDDsimulation/functions.py b/LDDsimulation/functions.py
index 5b60b2e45db446529b1136ac2b6f0ec470bfc24c..2dc8f54ca4b53be2ce4e855ffcc8a00e78716e69 100644
--- a/LDDsimulation/functions.py
+++ b/LDDsimulation/functions.py
@@ -33,18 +33,50 @@ import functools as ft
# RELATIVE PERMEABILITIES #####################################################
# Functions used as relative permeabilty functions for wetting phases
-def SpowN(S,N):
- return S**N
+def Monomial(S, phase, N):
+ if phase is "wetting":
+ return S**N
+ elif phase is "nonwetting":
+ return (1-S)**N
+ else:
+ raise(NotImplementedError())
-def SpowN_prime(S,N):
- return N*S**(N-1)
+def Monomial_prime(S, phase, N):
+ if phase is "wetting":
+ return N*S**(N-1)
+ elif phase is "nonwetting":
+ return -N*(1-S)**(N-1)
+ else:
+ raise(NotImplementedError())
-# 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 vanGenuchtenMualem(S, phase, n_index):
+ m = 1-1/n_index
+ if phase is "wetting":
+ return S**(1/2)*(1-(1-S**(1/m))**m)**2
+ elif phase is "nonwetting":
+ return (1-S)**(1/3)*(1-S**(1/m))**(2*m)
+ else:
+ raise(NotImplementedError())
+
+
+def vanGenuchtenMualem_prime(S, phase, n_index):
+ m = 1-1/n_index
+ if phase is "wetting":
+ part1 = S**(1/2)
+ part1_prime = 1/2*S**(-1/2)
+ part2 = (1-(1-S**(1/m))**m)**2
+ part2_prime = 2*(1-(1-S**(1/m))**m)*(-m*(1-S**(1/m))**(m-1))*(-1/m*S**(1/m-1))
+ return part1_prime*part2 + part1*part2_prime
+ elif phase is "nonwetting":
+ part1 = (1-S)**(1/3)
+ part1_prime = -1/3*(1-S)**(-2/3)
+ part2 = (1-S**(1/m))**(2*m)
+ part2_prime = 2*m*(1-S**(1/m))**(2*m-1)*(-1/m*S**(1/m-1))
+ return part1_prime*part2 + part1*part2_prime
+ else:
+ raise(NotImplementedError())
+
def generate_relative_permeability_dicts(
rel_perm_definition: tp.Dict[int, tp.Dict[str,str]]
@@ -68,37 +100,63 @@ def generate_relative_permeability_dicts(
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)}
- )
+ # the following isinstance handling is implemented to keep older
+ # uscase scripts running withould having to reimplement the definiton
+ # of the relative permeabilties. The more flexible way should be the
+ # one implemented in the isinstance(function_type, dict) way.
+ if isinstance(function_type, str):
+ if function_type is "Spow2":
+ output["ka"][subdomain].update(
+ {phase: ft.partial(Monomial, phase="wetting", N=2)}
+ )
+ output["ka_prime"][subdomain].update(
+ {phase: ft.partial(Monomial_prime, phase="wetting", N=2)}
+ )
+ elif function_type is "oneMinusSpow2":
+ output["ka"][subdomain].update(
+ {phase: ft.partial(Monomial, phase="nonwetting", N=2)}
+ )
+ output["ka_prime"][subdomain].update(
+ {phase: ft.partial(Monomial_prime, phase="nonwetting", N=2)}
+ )
+ elif function_type is "Spow3":
+ output["ka"][subdomain].update(
+ {phase: ft.partial(Monomial, phase="wetting", N=3)}
+ )
+ output["ka_prime"][subdomain].update(
+ {phase: ft.partial(Monomial_prime, phase="wetting", N=3)}
+ )
+ elif function_type is "oneMinusSpow3":
+ output["ka"][subdomain].update(
+ {phase: ft.partial(Monomial, phase="nonwetting", N=3)}
+ )
+ output["ka_prime"][subdomain].update(
+ {phase: ft.partial(Monomial_prime, phase="nonwetting", N=3)}
+ )
+ else:
+ raise(NotImplementedError())
+ elif isinstance(function_type, dict):
+ # this is the way things should be implemented.
+ for type, parameter in function_type.items():
+ if type is "monomial":
+ output["ka"][subdomain].update(
+ {phase: ft.partial(Monomial, phase=phase, N=parameter["power"])}
+ )
+ output["ka_prime"][subdomain].update(
+ {phase: ft.partial(Monomial_prime, phase=phase, N=parameter["power"])}
+ )
+ elif type is "vanGenuchtenMualem":
+ output["ka"][subdomain].update(
+ {phase: ft.partial(vanGenuchtenMualem, phase=phase, n_index=parameter["n"])}
+ )
+ output["ka_prime"][subdomain].update(
+ {phase: ft.partial(vanGenuchtenMualem_prime, phase=phase, n_index=parameter["n"])}
+ )
+ else:
+ raise(NotImplementedError())
+
else:
raise(NotImplementedError())
-
return output
# S-Pc RELATIONSHIPS ##########################################################
@@ -145,12 +203,12 @@ def test_S_prime_sym(pc, index):
# assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw=0
def vanGenuchten_S(pc, n_index, alpha):
# inverse capillary pressure-saturation-relationship for the simulation
- vanG = 1.0/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
+ vanG = (1.0/((1 + (alpha*pc)**n_index))**((n_index - 1)/n_index))
return df.conditional(pc > 0, vanG, 1)
def vanGenuchten_S_sym(pc, n_index, alpha):
# inverse capillary pressure-saturation-relationship
- return 1.0/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
+ return (1.0/((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.