From 6200e1d85ace1520c6f380181937c0b2e47debd8 Mon Sep 17 00:00:00 2001
From: David <forenkram@gmx.de>
Date: Mon, 29 Jun 2020 15:04:01 +0200
Subject: [PATCH] add vanGenuchten S-pc relationship. carefull, needs testing
---
LDDsimulation/functions.py | 42 +++++++++++++++++++++++++++++++++++++-
1 file changed, 41 insertions(+), 1 deletion(-)
diff --git a/LDDsimulation/functions.py b/LDDsimulation/functions.py
index 430dcb0..ab69bd7 100644
--- a/LDDsimulation/functions.py
+++ b/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())
--
GitLab