cleanup monomial rel perms and implement vanG-Mualem

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):
+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):
+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":
+            # 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(SpowN, N=2)}
+                        {phase: ft.partial(Monomial, phase="wetting", N=2)}
                         )
                     output["ka_prime"][subdomain].update(
-                    {phase: ft.partial(SpowN_prime, N=2)}
+                        {phase: ft.partial(Monomial_prime, phase="wetting", N=2)}
                         )
-            elif function_type == "oneMinusSpow2":
+                elif function_type is "oneMinusSpow2":
                     output["ka"][subdomain].update(
-                    {phase: ft.partial(OneMinusSpowN, N=2)}
+                        {phase: ft.partial(Monomial, phase="nonwetting", N=2)}
                         )
                     output["ka_prime"][subdomain].update(
-                    {phase: ft.partial(OneMinusSpowN_prime, N=2)}
+                        {phase: ft.partial(Monomial_prime, phase="nonwetting", N=2)}
                         )
-            elif function_type == "Spow3":
+                elif function_type is "Spow3":
                     output["ka"][subdomain].update(
-                    {phase: ft.partial(SpowN, N=3)}
+                        {phase: ft.partial(Monomial, phase="wetting", N=3)}
                         )
                     output["ka_prime"][subdomain].update(
-                    {phase: ft.partial(SpowN_prime, N=3)}
+                        {phase: ft.partial(Monomial_prime, phase="wetting", N=3)}
                         )
-            elif function_type == "oneMinusSpow3":
+                elif function_type is "oneMinusSpow3":
                     output["ka"][subdomain].update(
-                    {phase: ft.partial(OneMinusSpowN, N=3)}
+                        {phase: ft.partial(Monomial, phase="nonwetting", N=3)}
                         )
                     output["ka_prime"][subdomain].update(
-                    {phase: ft.partial(OneMinusSpowN_prime, N=3)}
+                        {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.