diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
index d0d78747e81577d7178e04fce5e472cf1a19b1a2..7b96ff47ebfc5f148caa738f9bd8a771df071986 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
@@ -2,15 +2,13 @@
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
@@ -187,176 +185,28 @@ intrinsic_permeability = {
"nonwetting": 1},
}
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return intrinsic_permeability[1]["wetting"]*s**2
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return intrinsic_permeability[1]["nonwetting"]*(1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-# ## relative permeabilty functions on subdomain 2
-# def rel_perm2w(s):
-# # relative permeabilty wetting on subdomain2
-# return s**3
-# def rel_perm2nw(s):
-# # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
-# return (1-s)**3
-#
-# _rel_perm2w = ft.partial(rel_perm2w)
-# _rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
- 1: subdomain1_rel_perm,
- 2: subdomain2_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]["wetting"]*2*s
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[1]["nonwetting"]*2*(1-s)
-
-# # definition of the derivatives of the relative permeabilities
-# # relative permeabilty functions on subdomain 1
-# def rel_perm2w_prime(s):
-# # relative permeabilty on subdomain1
-# return 3*s**2
-#
-# def rel_perm2nw_prime(s):
-# # relative permeabilty on subdomain1
-# return -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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=1),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=1),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=1),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 1}},
}
+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 #
@@ -397,6 +247,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,
@@ -453,6 +304,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
index 7e26fbfc5c077d05f0d79255d3291cb91764a39d..2b9d2cd8c4f2a30e39d2872a2aba5b2f7e6ccc25 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
@@ -2,15 +2,13 @@
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
@@ -166,9 +164,9 @@ date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
# Name of the usecase that will be printed during simulation.
-use_case = "TPR-2-patch-realistic"
+use_case = "TPR-2-patch-params-one-pure-dd"
# The name of this very file. Needed for creating log output.
-thisfile = "TP-R-2-patch-realistic.py"
+thisfile = "TP-R-2-patch-pure-dd.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
@@ -315,176 +313,35 @@ lambda_param = {#
'nonwetting': lambda_nw}
}
-## relative permeabilty functions on subdomain 1
-def rel_perm1w(s):
- # relative permeabilty wetting on subdomain1
- return s**2
-
-def rel_perm1nw(s):
- # relative permeabilty nonwetting on subdomain1
- return (1-s)**2
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-# ## relative permeabilty functions on subdomain 2
-# def rel_perm2w(s):
-# # relative permeabilty wetting on subdomain2
-# return s**3
-# def rel_perm2nw(s):
-# # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
-# return (1-s)**3
-#
-# _rel_perm2w = ft.partial(rel_perm2w)
-# _rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm
+intrinsic_permeability = {
+ 1: {"wetting": 1,
+ "nonwetting": 1},
+ 2: {"wetting": 1,
+ "nonwetting": 1},
}
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm1w_prime(s):
- # relative permeabilty on subdomain1
- return 2*s
-
-def rel_perm1nw_prime(s):
- # relative permeabilty on subdomain1
- return -2*(1-s)
-
-# # definition of the derivatives of the relative permeabilities
-# # relative permeabilty functions on subdomain 1
-# def rel_perm2w_prime(s):
-# # relative permeabilty on subdomain1
-# return 3*s**2
-#
-# def rel_perm2nw_prime(s):
-# # relative permeabilty on subdomain1
-# return -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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=1),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=1),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=1),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 1}},
}
+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 #
@@ -524,6 +381,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,
@@ -580,6 +438,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
index 10586745f8a82d2ad21f22b857751715f9b70352..fca9fe24f8fe9047647cd5abd6d052ef68ad0476 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py
@@ -3,19 +3,16 @@
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()
@@ -185,189 +182,28 @@ intrinsic_permeability = {
2: 0.01,
}
-
-# 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
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_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)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -381,7 +217,6 @@ p_e_sym = {
'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -401,6 +236,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,
@@ -457,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
index f2f48ffab2ff83c1f1a1cbc3c9fd1223302d289e..a2929e89bc6d38d9ea8d257936d208a3d7a804f5 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-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
@@ -183,189 +182,28 @@ intrinsic_permeability = {
2: 0.0001,
}
-
-# 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
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_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)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -398,6 +236,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,
@@ -454,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
index 92b123ffc013ce529113510ea5c37adcf068c34d..1b40691474c2305a209859a9a8288b6f227bf446 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-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
@@ -59,7 +58,7 @@ resolutions = {
# for each element t_0 in starttimes.
starttimes = [0.0]
timestep_size = 0.01
-number_of_timesteps = 800
+number_of_timesteps = 8
# LDD scheme parameters ######################################################
Lw1 = 0.025
@@ -72,7 +71,7 @@ lambda_w = 4
lambda_nw = 4
include_gravity = False
-debugflag = False
+debugflag = True
analyse_condition = False
# I/O CONFIG #################################################################
@@ -183,189 +182,28 @@ intrinsic_permeability = {
2: 0.01,
}
-
-# 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
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_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)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -398,6 +236,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,
@@ -454,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
index c95a5d15f6709b4617b40d15fd8f8e6bd2c287ef..6b16a1c022ef2f9256a99c46912a7a36f6097a4e 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py
@@ -3,15 +3,13 @@
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
@@ -184,177 +182,28 @@ intrinsic_permeability = {
2: 0.01,
}
-## 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
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,#
- 'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-# 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)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -368,7 +217,6 @@ p_e_sym = {
'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
-
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
@@ -388,6 +236,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,
@@ -444,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
index 113d0b8f29aed511d2f629a86387ea0bcbee74fe..d6bc8fe1f6d45036e486fe85f2ecba960f562bee 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.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
@@ -183,177 +182,28 @@ intrinsic_permeability = {
2: 1.0,
}
-## 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
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,#
- 'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-# 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)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -386,6 +236,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,
@@ -442,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
index cd3c85924cb8ca413ebf3374878a55c1c81b652a..c613e0bfc7add32c4b989f2f1d18b51365f88255 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.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
@@ -183,177 +182,28 @@ intrinsic_permeability = {
2: 1.0,
}
-## 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
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,#
- 'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-# 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)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -386,6 +236,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,
@@ -442,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
index b69793dce5c9b206822d34cad90d0be4d32f94a9..ea48fc1c03542a282e2059bdd7b96aafbdaca244 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-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
@@ -183,177 +182,28 @@ intrinsic_permeability = {
2: 1,
}
-## 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
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,#
- 'nonwetting': _rel_perm1nw
-}
-## relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,#
- 'nonwetting': _rel_perm2nw
-}
-
-## dictionary of relative permeabilties on all domains.
-relative_permeability = {#
- 1: subdomain1_rel_perm,
- 2: subdomain2_rel_perm
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-# 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)
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
-}
-
-
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -386,6 +236,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,
@@ -442,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
index 2862c098d3667e8547704f1ed13c005cf867ad80..877dae6b28bb540247178a59fa36adcd563fe8a7 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-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
@@ -183,189 +182,28 @@ intrinsic_permeability = {
2: 0.01,
}
-
-# 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
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_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)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -398,6 +236,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,
@@ -454,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
index 0b305e17438d6c65383736b4f03870973a53ab39..249ab3f371d28ad5ce167e4d17b9e4e478541583 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-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
@@ -183,189 +182,28 @@ intrinsic_permeability = {
2: 0.0001,
}
-
-# 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
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_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)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -398,6 +236,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,
@@ -454,6 +293,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,
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
index d054ce483da877240d7191dd9232c457d8e05dc9..54910805fab46bd492de45807caa0bb5c6c47156 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-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
@@ -183,189 +182,28 @@ intrinsic_permeability = {
2: 0.01,
}
-
-# 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
-
-
-_rel_perm1w = ft.partial(rel_perm1w)
-_rel_perm1nw = ft.partial(rel_perm1nw)
-
-subdomain1_rel_perm = {
- 'wetting': _rel_perm1w,
- 'nonwetting': _rel_perm1nw
-}
-
-# relative permeabilty functions on subdomain 2
-def rel_perm2w(s):
- # relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
-
-
-def rel_perm2nw(s):
- # relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
-
-
-_rel_perm2w = ft.partial(rel_perm2w)
-_rel_perm2nw = ft.partial(rel_perm2nw)
-
-subdomain2_rel_perm = {
- 'wetting': _rel_perm2w,
- 'nonwetting': _rel_perm2nw
-}
-
-# dictionary of relative permeabilties on all domains.
-relative_permeability = {
- 1: subdomain1_rel_perm,
- 2: subdomain2_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)
-
-
-# definition of the derivatives of the relative permeabilities
-# relative permeabilty functions on subdomain 1
-def rel_perm2w_prime(s):
- # relative permeabilty on subdomain2
- return intrinsic_permeability[2]*3*s**2
-
-
-def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain2
- return -3*intrinsic_permeability[2]*(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)
-
-subdomain1_rel_perm_prime = {
- 'wetting': _rel_perm1w_prime,
- 'nonwetting': _rel_perm1nw_prime
+# RELATIVE PEMRMEABILITIES
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
-subdomain2_rel_perm_prime = {
- 'wetting': _rel_perm2w_prime,
- 'nonwetting': _rel_perm2nw_prime
-}
-
-# dictionary of relative permeabilties on all domains.
-ka_prime = {
- 1: subdomain1_rel_perm_prime,
- 2: subdomain2_rel_perm_prime,
-}
-
-
-# def saturation1(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
-#
-# def saturation2(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
-#
-# # 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 saturation1_sym(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return 1/((1 + pc)**(1/(subdomain_index + 1)))
-#
-#
-# def saturation2_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 saturation1_sym_prime(pc, subdomain_index):
-# # inverse capillary pressure-saturation-relationship
-# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
-#
-#
-# def saturation2_sym_prime(pc, n_index, alpha):
-# # inverse capillary pressure-saturation-relationship
-# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
-#
-# # 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(saturation1_sym, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
-# }
-#
-# S_pc_sym_prime = {
-# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
-# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
-# }
-#
-# sat_pressure_relationship = {
-# 1: ft.partial(saturation1, subdomain_index = 1),#,
-# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
-# }
-
-
-def saturation(pc, index):
- # inverse capillary pressure-saturation-relationship
- return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
-
-
-def saturation_sym(pc, index):
- # inverse capillary pressure-saturation-relationship
- return 1/((1 + pc)**(1/(index + 1)))
-
-
-# derivative of S-pc relationship with respect to pc. This is needed for the
-# construction of a analytic solution.
-def saturation_sym_prime(pc, index):
- # inverse capillary pressure-saturation-relationship
- return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
-
-
-# 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, index=1),
- 2: ft.partial(saturation_sym, index=2),
- # 3: ft.partial(saturation_sym, index=2),
- # 4: ft.partial(saturation_sym, index=1)
-}
-
-S_pc_sym_prime = {
- 1: ft.partial(saturation_sym_prime, index=1),
- 2: ft.partial(saturation_sym_prime, index=2),
- # 3: ft.partial(saturation_sym_prime, index=2),
- # 4: ft.partial(saturation_sym_prime, index=1)
-}
-
-sat_pressure_relationship = {
- 1: ft.partial(saturation, index=1),
- 2: ft.partial(saturation, index=2),
- # 3: ft.partial(saturation, index=2),
- # 4: ft.partial(saturation, index=1)
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 2}},
}
+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 #
@@ -398,6 +236,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,
@@ -454,6 +293,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,