diff --git a/LDDsimulation/functions.py b/LDDsimulation/functions.py index 430dcb06146b2d3f551076788167a1b4da04c8f2..ab69bd7e5a3632ea1b58202bab71405417ceb379 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())