add vanGenuchten S-pc relationship. carefull, needs testing

LDDsimulation/functions.py

@@ -117,6 +117,26 @@ def test_S_prime_sym(pc, index):
     # inverse capillary pressure-saturation-relationship
     return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
 
+# Van Genchuten S-pc relationship
+# We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
+# 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))
+    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))
+
+# derivative of S-pc relationship with respect to pc. This is needed for the
+# construction of a analytic solution.
+def vanGenuchten_S_prime(pc, n_index, alpha):
+    # inverse capillary pressure-saturation-relationship
+    enumerator = -1.0*(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))
+    denominator = ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
+    return enumerator/denominator
+
 
 def generate_Spc_dicts(
         Spc_on_subdomains: tp.Dict[int, tp.Dict[str, tp.Dict[str, float]]]
@@ -161,7 +181,27 @@ def generate_Spc_dicts(
                         )},
                 )
             elif Spc_type == "vanGenuchten":
-                raise(NotImplementedError())
+                output["symbolic"].update(
+                    {subdomain: ft.partial(
+                        vanGenuchten_S_sym,
+                        n_index=parameters["n"],
+                        alpha=parameters["alpha"]
+                        )},
+                )
+                output["prime_symbolic"].update(
+                    {subdomain: ft.partial(
+                        vanGenuchten_S_prime_sym,
+                        n_index=parameters["n"],
+                        alpha=parameters["alpha"]
+                        )},
+                )
+                output["dolfin"].update(
+                    {subdomain: ft.partial(
+                        vanGenuchten_S,
+                        n_index=parameters["n"],
+                        alpha=parameters["alpha"]
+                        )},
+                )
             else:
                 raise(NotImplementedError())