diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
index 48766dab7d33ddb4d8aad5a46f5be49cb5a7998d..2db3def978f40226b2e3fa1c08f47d6681585688 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
@@ -5,14 +5,14 @@ This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
-import functools as ft
+import functions as fts
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
-import pandas as pd
import multiprocessing as mp
import domainSubstructuring as dss
+
# init sympy session
sym.init_printing()
@@ -223,246 +223,33 @@ intrinsic_permeability = {
4: 1.0, #0.01, #10e-3
}
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return intrinsic_permeability[1]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**2
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 3
-def rel_perm3w(s):
- # relative permeabilty wetting on subdomain3
- return intrinsic_permeability[3]*s**3
-
-
-def rel_perm3nw(s):
- # relative permeabilty nonwetting on subdomain3
- return intrinsic_permeability[3]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 4
-def rel_perm4w(s):
- # relative permeabilty wetting on subdomain4
- return intrinsic_permeability[4]*s**3
-
-
-def rel_perm4nw(s):
- # relative permeabilty nonwetting on subdomain4
- return intrinsic_permeability[4]*(1-s)**3
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-_rel_perm3w = ft.partial(rel_perm3w)
-_rel_perm3nw = ft.partial(rel_perm3nw)
-_rel_perm4w = ft.partial(rel_perm4w)
-_rel_perm4nw = ft.partial(rel_perm4nw)
-
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-subdomain3_rel_perm = {
- 'wetting': _rel_perm3w,
- 'nonwetting': _rel_perm3nw
-}
-
-subdomain4_rel_perm = {
- 'wetting': _rel_perm4w,
- 'nonwetting': _rel_perm4nw
-}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-# dictionary of relative permeabilties on all domains.
-# relative_permeability = {
-# 1: subdomain1_rel_perm,
-# 2: subdomain1_rel_perm,
-# 3: subdomain2_rel_perm,
-# 4: subdomain2_rel_perm
-# }
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm,
- 3: subdomain3_rel_perm,
- 4: subdomain4_rel_perm
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 2: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 3: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
+ 4: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
- # relative permeabilty on subdomain1
- return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*2*s
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -1*intrinsic_permeability[2]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 3
-def rel_perm3w_prime(s):
- # relative permeabilty on subdomain3
- return intrinsic_permeability[3]*3*s**2
-
-
-def rel_perm3nw_prime(s):
- # relative permeabilty on subdomain3
- return -1*intrinsic_permeability[3]*3*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 4
-def rel_perm4w_prime(s):
- # relative permeabilty on subdomain4
- return intrinsic_permeability[4]*3*s**2
-
-
-def rel_perm4nw_prime(s):
- # relative permeabilty on subdomain4
- return -1*intrinsic_permeability[4]*3*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
-_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
-_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
-_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-subdomain3_rel_perm_prime = {
- 'wetting': _rel_perm3w_prime,
- 'nonwetting': _rel_perm3nw_prime
-}
-
-
-subdomain4_rel_perm_prime = {
- 'wetting': _rel_perm4w_prime,
- 'nonwetting': _rel_perm4nw_prime
-}
-
-
-# dictionary of relative permeabilties on all domains.
-# ka_prime = {
-# 1: subdomain1_rel_perm_prime,
-# 2: subdomain1_rel_perm_prime,
-# 3: subdomain2_rel_perm_prime,
-# 4: subdomain2_rel_perm_prime
-# }
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
- 3: subdomain3_rel_perm_prime,
- 4: subdomain4_rel_perm_prime
-}
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-# this function needs to be monotonically decreasing in the capillary pressure
-# pc. Since in the richards case pc=-pw, this becomes as a function of pw a
-# mono tonically INCREASING function like in our Richards-Richards paper.
-# However since we unify the treatment in the code for Richards and two-phase,
-# we need the same requierment for both cases, two-phase and Richards.
-def saturation(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = df.conditional(
- pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
- return expr
-
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-def saturation_sym(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- # df.conditional(pc > 0,
- return 1/((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 saturation_sym_prime(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))\
- / ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
- return expr
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
-}
-
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
-}
-
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=3, alpha=0.001),
- 2: ft.partial(saturation, n_index=3, alpha=0.001),
- 3: ft.partial(saturation, n_index=6, alpha=0.001),
- 4: ft.partial(saturation, n_index=6, alpha=0.001)
-}
-
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
###############################################################################
# Manufacture source expressions with sympy #
@@ -472,10 +259,9 @@ t = sym.symbols('t', positive=True)
p_e_sym_2patch = {
1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
- 'nonwetting': 0.0*t}, # -1-t*(1.1 + y + x**2)**2},
+ 'nonwetting': 0.0*t},
2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x),
'nonwetting': (-2-t*(1.1+y + x**2))*y**2},
- # 'nonwetting': (-1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2))*y**2},
}
p_e_sym = {
@@ -489,30 +275,6 @@ p_e_sym = {
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
-
-# p_e_sym = {
-# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# -2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# - 2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 3: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# },
-# 4: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# }
-# }
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -532,6 +294,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
+ intrinsic_permeability=intrinsic_permeability,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
@@ -588,6 +351,7 @@ if __name__ == '__main__':
"L": L,
"lambda_param": lambda_param,
"relative_permeability": relative_permeability,
+ "intrinsic_permeability": intrinsic_permeability,
"sat_pressure_relationship": sat_pressure_relationship,
# "starttime": starttime,
"number_of_timesteps": number_of_timesteps,
@@ -618,7 +382,7 @@ if __name__ == '__main__':
)
LDDsim.start()
LDDsim.join()
- # run_simulation(
+ # hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
index fd9651e24c0a04f90e65a3f9141d284de79d83b3..269cfbcfc274b11511b3431823630b3f04893d4f 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
-import functools as ft
+import functions as fts
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
-import pandas as pd
import multiprocessing as mp
import domainSubstructuring as dss
@@ -239,246 +238,33 @@ intrinsic_permeability = {
4: 0.01, # 10e-3
}
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return intrinsic_permeability[1]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**2
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 3
-def rel_perm3w(s):
- # relative permeabilty wetting on subdomain3
- return intrinsic_permeability[3]*s**3
-
-
-def rel_perm3nw(s):
- # relative permeabilty nonwetting on subdomain3
- return intrinsic_permeability[3]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 4
-def rel_perm4w(s):
- # relative permeabilty wetting on subdomain4
- return intrinsic_permeability[4]*s**3
-
-
-def rel_perm4nw(s):
- # relative permeabilty nonwetting on subdomain4
- return intrinsic_permeability[4]*(1-s)**3
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-_rel_perm3w = ft.partial(rel_perm3w)
-_rel_perm3nw = ft.partial(rel_perm3nw)
-_rel_perm4w = ft.partial(rel_perm4w)
-_rel_perm4nw = ft.partial(rel_perm4nw)
-
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain3_rel_perm = {
- 'wetting': _rel_perm3w,
- 'nonwetting': _rel_perm3nw
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 2: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 3: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
+ 4: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
}
-
-subdomain4_rel_perm = {
- 'wetting': _rel_perm4w,
- 'nonwetting': _rel_perm4nw
-}
-
-# dictionary of relative permeabilties on all domains.
-# relative_permeability = {
-# 1: subdomain1_rel_perm,
-# 2: subdomain1_rel_perm,
-# 3: subdomain2_rel_perm,
-# 4: subdomain2_rel_perm
-# }
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm,
- 3: subdomain3_rel_perm,
- 4: subdomain4_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
- # relative permeabilty on subdomain1
- return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*2*s
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -1*intrinsic_permeability[2]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 3
-def rel_perm3w_prime(s):
- # relative permeabilty on subdomain3
- return intrinsic_permeability[3]*3*s**2
-
-
-def rel_perm3nw_prime(s):
- # relative permeabilty on subdomain3
- return -1*intrinsic_permeability[3]*3*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 4
-def rel_perm4w_prime(s):
- # relative permeabilty on subdomain4
- return intrinsic_permeability[4]*3*s**2
-
-
-def rel_perm4nw_prime(s):
- # relative permeabilty on subdomain4
- return -1*intrinsic_permeability[4]*3*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
-_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
-_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
-_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-subdomain3_rel_perm_prime = {
- 'wetting': _rel_perm3w_prime,
- 'nonwetting': _rel_perm3nw_prime
-}
-
-
-subdomain4_rel_perm_prime = {
- 'wetting': _rel_perm4w_prime,
- 'nonwetting': _rel_perm4nw_prime
-}
-
-
-# dictionary of relative permeabilties on all domains.
-# ka_prime = {
-# 1: subdomain1_rel_perm_prime,
-# 2: subdomain1_rel_perm_prime,
-# 3: subdomain2_rel_perm_prime,
-# 4: subdomain2_rel_perm_prime
-# }
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
- 3: subdomain3_rel_perm_prime,
- 4: subdomain4_rel_perm_prime
-}
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-# this function needs to be monotonically decreasing in the capillary pressure
-# pc. Since in the richards case pc=-pw, this becomes as a function of pw a
-# mono tonically INCREASING function like in our Richards-Richards paper.
-# However since we unify the treatment in the code for Richards and two-phase,
-# we need the same requierment for both cases, two-phase and Richards.
-def saturation(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = df.conditional(
- pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
- return expr
-
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-def saturation_sym(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- # df.conditional(pc > 0,
- return 1/((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 saturation_sym_prime(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))\
- / ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
- return expr
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
-}
-
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
-}
-
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=3, alpha=0.001),
- 2: ft.partial(saturation, n_index=3, alpha=0.001),
- 3: ft.partial(saturation, n_index=6, alpha=0.001),
- 4: ft.partial(saturation, n_index=6, alpha=0.001)
-}
-
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
#############################################
# Manufacture source expressions with sympy #
@@ -505,30 +291,6 @@ p_e_sym = {
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
-
-# p_e_sym = {
-# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# -2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# - 2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 3: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# },
-# 4: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# }
-# }
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -548,6 +310,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
+ intrinsic_permeability=intrinsic_permeability,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
@@ -604,6 +367,7 @@ if __name__ == '__main__':
"L": L,
"lambda_param": lambda_param,
"relative_permeability": relative_permeability,
+ "intrinsic_permeability": intrinsic_permeability,
"sat_pressure_relationship": sat_pressure_relationship,
# "starttime": starttime,
"number_of_timesteps": number_of_timesteps,
@@ -634,7 +398,7 @@ if __name__ == '__main__':
)
LDDsim.start()
LDDsim.join()
- # run_simulation(
+ # hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
index 1916d58fd6edf8a3ec64d61a6352fe41d0f76a8f..013e15beca6f60b1f7c5713e1a6e5e81d4940dee 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
-import functools as ft
+import functions as fts
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
-import pandas as pd
import multiprocessing as mp
import domainSubstructuring as dss
@@ -219,16 +218,6 @@ lambda_param = {
2: {'wetting': lambda34_w,
'nonwetting': lambda34_nw},
}
-# lambda_param = {
-# 1: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 2: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 3: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 4: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# }
# after Lewis, see pdf file
intrinsic_permeability = {
@@ -238,246 +227,33 @@ intrinsic_permeability = {
4: 0.01, # 10e-3
}
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return intrinsic_permeability[1]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**2
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 3
-def rel_perm3w(s):
- # relative permeabilty wetting on subdomain3
- return intrinsic_permeability[3]*s**3
-
-
-def rel_perm3nw(s):
- # relative permeabilty nonwetting on subdomain3
- return intrinsic_permeability[3]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 4
-def rel_perm4w(s):
- # relative permeabilty wetting on subdomain4
- return intrinsic_permeability[4]*s**3
-
-
-def rel_perm4nw(s):
- # relative permeabilty nonwetting on subdomain4
- return intrinsic_permeability[4]*(1-s)**3
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-_rel_perm3w = ft.partial(rel_perm3w)
-_rel_perm3nw = ft.partial(rel_perm3nw)
-_rel_perm4w = ft.partial(rel_perm4w)
-_rel_perm4nw = ft.partial(rel_perm4nw)
-
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-subdomain3_rel_perm = {
- 'wetting': _rel_perm3w,
- 'nonwetting': _rel_perm3nw
-}
-
-subdomain4_rel_perm = {
- 'wetting': _rel_perm4w,
- 'nonwetting': _rel_perm4nw
-}
-
-# dictionary of relative permeabilties on all domains.
-# relative_permeability = {
-# 1: subdomain1_rel_perm,
-# 2: subdomain1_rel_perm,
-# 3: subdomain2_rel_perm,
-# 4: subdomain2_rel_perm
-# }
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm,
- 3: subdomain3_rel_perm,
- 4: subdomain4_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
- # relative permeabilty on subdomain1
- return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*2*s
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -1*intrinsic_permeability[2]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 3
-def rel_perm3w_prime(s):
- # relative permeabilty on subdomain3
- return intrinsic_permeability[3]*3*s**2
-
-
-def rel_perm3nw_prime(s):
- # relative permeabilty on subdomain3
- return -1*intrinsic_permeability[3]*3*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 4
-def rel_perm4w_prime(s):
- # relative permeabilty on subdomain4
- return intrinsic_permeability[4]*3*s**2
-
-
-def rel_perm4nw_prime(s):
- # relative permeabilty on subdomain4
- return -1*intrinsic_permeability[4]*3*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
-_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
-_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
-_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-subdomain3_rel_perm_prime = {
- 'wetting': _rel_perm3w_prime,
- 'nonwetting': _rel_perm3nw_prime
-}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-
-subdomain4_rel_perm_prime = {
- 'wetting': _rel_perm4w_prime,
- 'nonwetting': _rel_perm4nw_prime
-}
-
-
-# dictionary of relative permeabilties on all domains.
-# ka_prime = {
-# 1: subdomain1_rel_perm_prime,
-# 2: subdomain1_rel_perm_prime,
-# 3: subdomain2_rel_perm_prime,
-# 4: subdomain2_rel_perm_prime
-# }
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
- 3: subdomain3_rel_perm_prime,
- 4: subdomain4_rel_perm_prime
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 2: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 3: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
+ 4: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
}
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-# this function needs to be monotonically decreasing in the capillary pressure
-# pc. Since in the richards case pc=-pw, this becomes as a function of pw a
-# mono tonically INCREASING function like in our Richards-Richards paper.
-# However since we unify the treatment in the code for Richards and two-phase,
-# we need the same requierment for both cases, two-phase and Richards.
-def saturation(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = df.conditional(
- pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
- return expr
-
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-def saturation_sym(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- # df.conditional(pc > 0,
- return 1/((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 saturation_sym_prime(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))\
- / ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
- return expr
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
-}
-
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
-}
-
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=3, alpha=0.001),
- 2: ft.partial(saturation, n_index=3, alpha=0.001),
- 3: ft.partial(saturation, n_index=6, alpha=0.001),
- 4: ft.partial(saturation, n_index=6, alpha=0.001)
-}
-
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
#############################################
# Manufacture source expressions with sympy #
@@ -504,30 +280,6 @@ p_e_sym = {
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
-
-# p_e_sym = {
-# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# -2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# - 2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 3: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# },
-# 4: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# }
-# }
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -547,6 +299,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
+ intrinsic_permeability=intrinsic_permeability,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
@@ -603,6 +356,7 @@ if __name__ == '__main__':
"L": L,
"lambda_param": lambda_param,
"relative_permeability": relative_permeability,
+ "intrinsic_permeability": intrinsic_permeability,
"sat_pressure_relationship": sat_pressure_relationship,
# "starttime": starttime,
"number_of_timesteps": number_of_timesteps,
@@ -633,7 +387,7 @@ if __name__ == '__main__':
)
LDDsim.start()
LDDsim.join()
- # run_simulation(
+ # hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
index c115fa15ed2638fa20b14bf6e466e8008c68ded7..babbf0e2f639be1a2a1131943b33dcbcb59c8681 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
-import functools as ft
+import functions as fts
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
-import pandas as pd
import multiprocessing as mp
import domainSubstructuring as dss
@@ -235,246 +234,33 @@ intrinsic_permeability = {
4: 0.01, #10e-3
}
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return intrinsic_permeability[1]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**2
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 3
-def rel_perm3w(s):
- # relative permeabilty wetting on subdomain3
- return intrinsic_permeability[3]*s**3
-
-
-def rel_perm3nw(s):
- # relative permeabilty nonwetting on subdomain3
- return intrinsic_permeability[3]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 4
-def rel_perm4w(s):
- # relative permeabilty wetting on subdomain4
- return intrinsic_permeability[4]*s**3
-
-
-def rel_perm4nw(s):
- # relative permeabilty nonwetting on subdomain4
- return intrinsic_permeability[4]*(1-s)**3
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-_rel_perm3w = ft.partial(rel_perm3w)
-_rel_perm3nw = ft.partial(rel_perm3nw)
-_rel_perm4w = ft.partial(rel_perm4w)
-_rel_perm4nw = ft.partial(rel_perm4nw)
-
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain3_rel_perm = {
- 'wetting': _rel_perm3w,
- 'nonwetting': _rel_perm3nw
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 2: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 3: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
+ 4: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
}
-
-subdomain4_rel_perm = {
- 'wetting': _rel_perm4w,
- 'nonwetting': _rel_perm4nw
-}
-
-# dictionary of relative permeabilties on all domains.
-# relative_permeability = {
-# 1: subdomain1_rel_perm,
-# 2: subdomain1_rel_perm,
-# 3: subdomain2_rel_perm,
-# 4: subdomain2_rel_perm
-# }
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm,
- 3: subdomain3_rel_perm,
- 4: subdomain4_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
- # relative permeabilty on subdomain1
- return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*2*s
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -1*intrinsic_permeability[2]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 3
-def rel_perm3w_prime(s):
- # relative permeabilty on subdomain3
- return intrinsic_permeability[3]*3*s**2
-
-
-def rel_perm3nw_prime(s):
- # relative permeabilty on subdomain3
- return -1*intrinsic_permeability[3]*3*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 4
-def rel_perm4w_prime(s):
- # relative permeabilty on subdomain4
- return intrinsic_permeability[4]*3*s**2
-
-
-def rel_perm4nw_prime(s):
- # relative permeabilty on subdomain4
- return -1*intrinsic_permeability[4]*3*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
-_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
-_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
-_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-subdomain3_rel_perm_prime = {
- 'wetting': _rel_perm3w_prime,
- 'nonwetting': _rel_perm3nw_prime
-}
-
-
-subdomain4_rel_perm_prime = {
- 'wetting': _rel_perm4w_prime,
- 'nonwetting': _rel_perm4nw_prime
-}
-
-
-# dictionary of relative permeabilties on all domains.
-# ka_prime = {
-# 1: subdomain1_rel_perm_prime,
-# 2: subdomain1_rel_perm_prime,
-# 3: subdomain2_rel_perm_prime,
-# 4: subdomain2_rel_perm_prime
-# }
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
- 3: subdomain3_rel_perm_prime,
- 4: subdomain4_rel_perm_prime
-}
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-# this function needs to be monotonically decreasing in the capillary pressure
-# pc. Since in the richards case pc=-pw, this becomes as a function of pw a
-# mono tonically INCREASING function like in our Richards-Richards paper.
-# However since we unify the treatment in the code for Richards and two-phase,
-# we need the same requierment for both cases, two-phase and Richards.
-def saturation(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = df.conditional(
- pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
- return expr
-
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-def saturation_sym(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- # df.conditional(pc > 0,
- return 1/((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 saturation_sym_prime(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))\
- / ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
- return expr
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
-}
-
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
-}
-
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=3, alpha=0.001),
- 2: ft.partial(saturation, n_index=3, alpha=0.001),
- 3: ft.partial(saturation, n_index=6, alpha=0.001),
- 4: ft.partial(saturation, n_index=6, alpha=0.001)
-}
-
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
###############################################################################
# Manufacture source expressions with sympy #
@@ -501,30 +287,6 @@ p_e_sym = {
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
-
-# p_e_sym = {
-# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# -2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# - 2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 3: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# },
-# 4: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# }
-# }
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -544,6 +306,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
+ intrinsic_permeability=intrinsic_permeability,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
@@ -600,6 +363,7 @@ if __name__ == '__main__':
"L": L,
"lambda_param": lambda_param,
"relative_permeability": relative_permeability,
+ "intrinsic_permeability": intrinsic_permeability,
"sat_pressure_relationship": sat_pressure_relationship,
# "starttime": starttime,
"number_of_timesteps": number_of_timesteps,
@@ -630,7 +394,7 @@ if __name__ == '__main__':
)
LDDsim.start()
LDDsim.join()
- # run_simulation(
+ # hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
diff --git a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.py
index 3c14df1e9159379de770ec5f8c9813516acc9aca..68119eb8a90c6fadb7553341c786f809c15b3aa3 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
-import functools as ft
+import functions as fts
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
-import pandas as pd
import multiprocessing as mp
import domainSubstructuring as dss
@@ -217,16 +216,6 @@ lambda_param = {
2: {'wetting': lambda34_w,
'nonwetting': lambda34_nw},
}
-# lambda_param = {
-# 1: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 2: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 3: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 4: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# }
# after Lewis, see pdf file
intrinsic_permeability = {
@@ -236,246 +225,33 @@ intrinsic_permeability = {
4: 0.01, # 10e-3
}
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return intrinsic_permeability[1]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**2
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 3
-def rel_perm3w(s):
- # relative permeabilty wetting on subdomain3
- return intrinsic_permeability[3]*s**3
-
-
-def rel_perm3nw(s):
- # relative permeabilty nonwetting on subdomain3
- return intrinsic_permeability[3]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 4
-def rel_perm4w(s):
- # relative permeabilty wetting on subdomain4
- return intrinsic_permeability[4]*s**3
-
-
-def rel_perm4nw(s):
- # relative permeabilty nonwetting on subdomain4
- return intrinsic_permeability[4]*(1-s)**3
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-_rel_perm3w = ft.partial(rel_perm3w)
-_rel_perm3nw = ft.partial(rel_perm3nw)
-_rel_perm4w = ft.partial(rel_perm4w)
-_rel_perm4nw = ft.partial(rel_perm4nw)
-
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-subdomain3_rel_perm = {
- 'wetting': _rel_perm3w,
- 'nonwetting': _rel_perm3nw
-}
-
-subdomain4_rel_perm = {
- 'wetting': _rel_perm4w,
- 'nonwetting': _rel_perm4nw
-}
-
-# dictionary of relative permeabilties on all domains.
-# relative_permeability = {
-# 1: subdomain1_rel_perm,
-# 2: subdomain1_rel_perm,
-# 3: subdomain2_rel_perm,
-# 4: subdomain2_rel_perm
-# }
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm,
- 3: subdomain3_rel_perm,
- 4: subdomain4_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
- # relative permeabilty on subdomain1
- return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*2*s
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -1*intrinsic_permeability[2]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 3
-def rel_perm3w_prime(s):
- # relative permeabilty on subdomain3
- return intrinsic_permeability[3]*3*s**2
-
-
-def rel_perm3nw_prime(s):
- # relative permeabilty on subdomain3
- return -1*intrinsic_permeability[3]*3*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 4
-def rel_perm4w_prime(s):
- # relative permeabilty on subdomain4
- return intrinsic_permeability[4]*3*s**2
-
-
-def rel_perm4nw_prime(s):
- # relative permeabilty on subdomain4
- return -1*intrinsic_permeability[4]*3*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
-_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
-_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
-_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-subdomain3_rel_perm_prime = {
- 'wetting': _rel_perm3w_prime,
- 'nonwetting': _rel_perm3nw_prime
-}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-
-subdomain4_rel_perm_prime = {
- 'wetting': _rel_perm4w_prime,
- 'nonwetting': _rel_perm4nw_prime
-}
-
-
-# dictionary of relative permeabilties on all domains.
-# ka_prime = {
-# 1: subdomain1_rel_perm_prime,
-# 2: subdomain1_rel_perm_prime,
-# 3: subdomain2_rel_perm_prime,
-# 4: subdomain2_rel_perm_prime
-# }
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
- 3: subdomain3_rel_perm_prime,
- 4: subdomain4_rel_perm_prime
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 2: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 3: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
+ 4: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
}
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-# this function needs to be monotonically decreasing in the capillary pressure
-# pc. Since in the richards case pc=-pw, this becomes as a function of pw a
-# mono tonically INCREASING function like in our Richards-Richards paper.
-# However since we unify the treatment in the code for Richards and two-phase,
-# we need the same requierment for both cases, two-phase and Richards.
-def saturation(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = df.conditional(
- pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
- return expr
-
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-def saturation_sym(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- # df.conditional(pc > 0,
- return 1/((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 saturation_sym_prime(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))\
- / ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
- return expr
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
-}
-
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
-}
-
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=3, alpha=0.001),
- 2: ft.partial(saturation, n_index=3, alpha=0.001),
- 3: ft.partial(saturation, n_index=6, alpha=0.001),
- 4: ft.partial(saturation, n_index=6, alpha=0.001)
-}
-
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
#############################################
# Manufacture source expressions with sympy #
@@ -485,10 +261,9 @@ t = sym.symbols('t', positive=True)
p_e_sym_2patch = {
1: {'wetting': -6 - (1+t*t)*(1 + x*x + y*y),
- 'nonwetting': 0.0*t}, # -1-t*(1.1 + y + x**2)**2},
+ 'nonwetting': 0.0*t},
2: {'wetting': -6.0 - (1.0 + t*t)*(1.0 + x*x),
'nonwetting': (-1-t*(1.1+y + x**2))*y**2},
- # 'nonwetting': (-1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2))*y**2},
}
p_e_sym = {
@@ -502,30 +277,6 @@ p_e_sym = {
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
-
-# p_e_sym = {
-# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# -2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# - 2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 3: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# },
-# 4: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# }
-# }
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -545,6 +296,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
+ intrinsic_permeability=intrinsic_permeability,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
@@ -601,6 +353,7 @@ if __name__ == '__main__':
"L": L,
"lambda_param": lambda_param,
"relative_permeability": relative_permeability,
+ "intrinsic_permeability": intrinsic_permeability,
"sat_pressure_relationship": sat_pressure_relationship,
# "starttime": starttime,
"number_of_timesteps": number_of_timesteps,
@@ -631,7 +384,7 @@ if __name__ == '__main__':
)
LDDsim.start()
LDDsim.join()
- # run_simulation(
+ # hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
diff --git a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.py
index bee65c9a8806119abdf2e3ec4c193bf51b3d52c8..a73c781d22d67a7dc238d8ec922dbb60614d8218 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.py
@@ -5,12 +5,11 @@ This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
-import functools as ft
+import functions as fts
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
-import pandas as pd
import multiprocessing as mp
import domainSubstructuring as dss
@@ -43,14 +42,14 @@ FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = True
resolutions = {
- 1: 9e-6, # h=2
- 2: 9e-6, # h=1.1180
- 4: 9e-6, # h=0.5590
- 8: 9e-6, # h=0.2814
+ # 1: 9e-6, # h=2
+ # 2: 9e-6, # h=1.1180
+ # 4: 9e-6, # h=0.5590
+ # 8: 9e-6, # h=0.2814
16: 5e-6, # h=0.1412
- 32: 4e-6,
- 64: 1e-6,
- 128: 1e-6
+ # 32: 4e-6,
+ # 64: 1e-6,
+ # 128: 1e-6
}
# starttimes gives a list of starttimes to run the simulation from.
@@ -58,7 +57,7 @@ resolutions = {
# for each element t_0 in starttimes.
starttimes = [0.0]
timestep_size = 0.001
-number_of_timesteps = 1000
+number_of_timesteps = 5
# LDD scheme parameters ######################################################
@@ -79,7 +78,7 @@ lambda34_w = 4
lambda34_nw = 4
include_gravity = False
-debugflag = False
+debugflag = True
analyse_condition = False
# I/O CONFIG #################################################################
@@ -236,246 +235,33 @@ intrinsic_permeability = {
4: 0.001, # 10e-3
}
-
-# relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return intrinsic_permeability[1]*s**2
-
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return intrinsic_permeability[1]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**2
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**2
-
-
-# relative permeabilty functions on subdomain 3
-def rel_perm3w(s):
- # relative permeabilty wetting on subdomain3
- return intrinsic_permeability[3]*s**3
-
-
-def rel_perm3nw(s):
- # relative permeabilty nonwetting on subdomain3
- return intrinsic_permeability[3]*(1-s)**3
-
-
-# relative permeabilty functions on subdomain 4
-def rel_perm4w(s):
- # relative permeabilty wetting on subdomain4
- return intrinsic_permeability[4]*s**3
-
-
-def rel_perm4nw(s):
- # relative permeabilty nonwetting on subdomain4
- return intrinsic_permeability[4]*(1-s)**3
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-_rel_perm3w = ft.partial(rel_perm3w)
-_rel_perm3nw = ft.partial(rel_perm3nw)
-_rel_perm4w = ft.partial(rel_perm4w)
-_rel_perm4nw = ft.partial(rel_perm4nw)
-
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain3_rel_perm = {
- 'wetting': _rel_perm3w,
- 'nonwetting': _rel_perm3nw
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 2: {"vanGenuchten": {"n": 3, "alpha": 0.001}},
+ 3: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
+ 4: {"vanGenuchten": {"n": 6, "alpha": 0.001}},
}
-
-subdomain4_rel_perm = {
- 'wetting': _rel_perm4w,
- 'nonwetting': _rel_perm4nw
-}
-
-# dictionary of relative permeabilties on all domains.
-# relative_permeability = {
-# 1: subdomain1_rel_perm,
-# 2: subdomain1_rel_perm,
-# 3: subdomain2_rel_perm,
-# 4: subdomain2_rel_perm
-# }
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm,
- 3: subdomain3_rel_perm,
- 4: subdomain4_rel_perm
-}
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
- # relative permeabilty on subdomain1
- return intrinsic_permeability[1]*2*s
-
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[1]*2*(1-s)
-
-
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*2*s
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -1*intrinsic_permeability[2]*2*(1-s)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 3
-def rel_perm3w_prime(s):
- # relative permeabilty on subdomain3
- return intrinsic_permeability[3]*3*s**2
-
-
-def rel_perm3nw_prime(s):
- # relative permeabilty on subdomain3
- return -1*intrinsic_permeability[3]*3*(1-s)**2
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 4
-def rel_perm4w_prime(s):
- # relative permeabilty on subdomain4
- return intrinsic_permeability[4]*3*s**2
-
-
-def rel_perm4nw_prime(s):
- # relative permeabilty on subdomain4
- return -1*intrinsic_permeability[4]*3*(1-s)**2
-
-
-_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
-_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
-_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
-_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
-_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
-_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
-_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
-_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-subdomain3_rel_perm_prime = {
- 'wetting': _rel_perm3w_prime,
- 'nonwetting': _rel_perm3nw_prime
-}
-
-
-subdomain4_rel_perm_prime = {
- 'wetting': _rel_perm4w_prime,
- 'nonwetting': _rel_perm4nw_prime
-}
-
-
-# dictionary of relative permeabilties on all domains.
-# ka_prime = {
-# 1: subdomain1_rel_perm_prime,
-# 2: subdomain1_rel_perm_prime,
-# 3: subdomain2_rel_perm_prime,
-# 4: subdomain2_rel_perm_prime
-# }
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
- 3: subdomain3_rel_perm_prime,
- 4: subdomain4_rel_perm_prime
-}
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-# this function needs to be monotonically decreasing in the capillary pressure
-# pc. Since in the richards case pc=-pw, this becomes as a function of pw a
-# mono tonically INCREASING function like in our Richards-Richards paper.
-# However since we unify the treatment in the code for Richards and two-phase,
-# we need the same requierment for both cases, two-phase and Richards.
-def saturation(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = df.conditional(
- pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
- return expr
-
-
-# S-pc-relation ship. We use the van Genuchten approach, i.e.
-# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
-# (see Helmig) and assume that residual saturation is Sw
-def saturation_sym(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- # df.conditional(pc > 0,
- return 1/((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 saturation_sym_prime(pc, n_index, alpha):
- # inverse capillary pressure-saturation-relationship
- expr = -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1))\
- / ((1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index))
- return expr
-
-
-# note that the conditional definition of S-pc in the nonsymbolic part will be
-# incorporated in the construction of the exact solution below.
-S_pc_sym = {
- 1: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym, n_index=6, alpha=0.001)
-}
-
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
- 3: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001),
- 4: ft.partial(saturation_sym_prime, n_index=6, alpha=0.001)
-}
-
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, n_index=3, alpha=0.001),
- 2: ft.partial(saturation, n_index=3, alpha=0.001),
- 3: ft.partial(saturation, n_index=6, alpha=0.001),
- 4: ft.partial(saturation, n_index=6, alpha=0.001)
-}
-
+Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
+S_pc_sym = Spc["symbolic"]
+S_pc_sym_prime = Spc["prime_symbolic"]
+sat_pressure_relationship = Spc["dolfin"]
#############################################
# Manufacture source expressions with sympy #
@@ -502,30 +288,6 @@ p_e_sym = {
'nonwetting': p_e_sym_2patch[2]['nonwetting']}
}
-
-# p_e_sym = {
-# 1: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# -2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
-# 'nonwetting':
-# - 2 - t*(1 + (y - 5.0) + x**2)**2
-# - sym.sqrt(2 + t**2)*(1 + (y - 5.0))
-# },
-# 3: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# },
-# 4: {'wetting':
-# 1.0 - (1.0 + t*t)*(10.0 + x*x + (y - 5.0)*(y - 5.0))
-# - (y - 5.0)*(y - 5.0)*3*sym.sin(-2*t + 2*x)*sym.sin(1/2*y - 1.2*t),
-# 'nonwetting': - 2 - t*(1 + x**2)**2 - sym.sqrt(2 + t**2)
-# }
-# }
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -545,6 +307,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions(
saturation_pressure_relationship_prime=S_pc_sym_prime,
viscosity=viscosity,
porosity=porosity,
+ intrinsic_permeability=intrinsic_permeability,
relative_permeability=relative_permeability,
relative_permeability_prime=ka_prime,
densities=densities,
@@ -601,6 +364,7 @@ if __name__ == '__main__':
"L": L,
"lambda_param": lambda_param,
"relative_permeability": relative_permeability,
+ "intrinsic_permeability": intrinsic_permeability,
"sat_pressure_relationship": sat_pressure_relationship,
# "starttime": starttime,
"number_of_timesteps": number_of_timesteps,
@@ -631,7 +395,7 @@ if __name__ == '__main__':
)
LDDsim.start()
LDDsim.join()
- # run_simulation(
+ # hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter