diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
index 468ebb8bc0691caacc5c7261a5e19767607d8df5..3b054bfd98e27d85545e780d47ab1d8b48f7a16b 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py
@@ -1,23 +1,24 @@
#!/usr/bin/python3
-"""Layered soil simulation.
+""" TP-R Layered soil simulation.
This program sets up an LDD simulation
"""
import dolfin as df
-# import mshr
-# import numpy as np
import sympy as sym
-# import typing as tp
import functools as ft
-# import domainPatch as dp
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
import pandas as pd
-# check if output directory exists
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
@@ -25,42 +26,51 @@ else:
print("Directory ", './output', " already exists. Will use as output \
directory")
-
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
-# init sympy session
-sym.init_printing()
-# solver_tol = 6E-7
-use_case = "TP-R-layered-soil-all-params-set-one-new-lambda"
+
+# Name of the usecase that will be printed during simulation.
+use_case = "TP-R-layered-soil-all-params-one"
+# The name of this very file. Needed for creating log output.
thisfile = "TP-R-layered_soil-all-params-one.py"
-max_iter_num = 500
+# GENERAL SOLVER CONFIG ######################################################
+# maximal iteration per timestep
+max_iter_num = 300
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False
resolutions = {
- # 1: 1e-7, # h=2
- # 2: 2e-5, # h=1.1180
- # 4: 1e-6, # h=0.5590
- # 8: 1e-6, # h=0.2814
- # 16: 1e-6, # h=0.1412
- 32: 1e-6,
- # 64: 5e-7,
- # 128: 5e-7
+ # 1: 1e-6,
+ # 2: 1e-6,
+ # 4: 1e-6,
+ # 8: 1e-6,
+ 16: 5e-6,
+ # 32: 5e-6,
+ # 64: 2e-6,
+ # 128: 1e-6,
+ # 256: 1e-6,
}
-# GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.001
-number_of_timesteps = 5
-plot_timestep_every = 4
-# decide how many timesteps you want analysed. Analysed means, that we write
-# out subsequent errors of the L-iteration within the timestep.
-number_of_timesteps_to_analyse = 6
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
starttimes = [0.0]
+timestep_size = 0.001
+number_of_timesteps = 20
+
-Lw = 0.025 # /timestep_size
-Lnw = Lw
+# LDD scheme parameters ######################################################
+Lw1 = 0.025 # /timestep_size
+Lnw1 = Lw1
+Lw2 = 0.025 # /timestep_size
+Lnw2 = Lw2
+Lw3 = 0.025 # /timestep_size
+Lnw3 = Lw3
+Lw4 = 0.025 # /timestep_size
+Lnw4 = Lw4
lambda12_w = 40
lambda12_nw = 40
@@ -69,31 +79,42 @@ lambda23_nw = 40
lambda34_w = 40
lambda34_nw = 40
-include_gravity = True
+include_gravity = False
debugflag = False
analyse_condition = True
-if mesh_study:
- output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-else:
- for tol in resolutions.values():
- solver_tol = tol
- output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
- )
-
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
'space_errornorms': True,
+ # save the mesh and marker functions to disk
'meshes_and_markers': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True
}
else:
@@ -107,7 +128,20 @@ else:
'subsequent_errors': True
}
+# OUTPUT FILE STRING #########################################################
+if mesh_study:
+ output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+else:
+ for tol in resolutions.values():
+ solver_tol = tol
+ output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
+ )
+
+# DOMAIN AND INTERFACE #######################################################
# global domain
subdomain0_vertices = [
df.Point(-1.0, -1.0),
@@ -257,6 +291,8 @@ outer_boundary_def_points = {
4: subdomain4_outer_boundary_verts
}
+
+# MODEL CONFIGURATION #########################################################
isRichards = {
1: True,
2: True,
@@ -273,53 +309,51 @@ isRichards = {
# Dict of the form: { subdom_num : viscosity }
viscosity = {
- 1: {'wetting': 1,
- 'nonwetting': 1},
- 2: {'wetting': 1,
- 'nonwetting': 1},
- 3: {'wetting': 1,
- 'nonwetting': 1},
- 4: {'wetting': 1,
- 'nonwetting': 1},
+ 1: {'wetting': 1.0,
+ 'nonwetting': 1.0},
+ 2: {'wetting': 1.0,
+ 'nonwetting': 1.0},
+ 3: {'wetting': 1.0,
+ 'nonwetting': 1.0},
+ 4: {'wetting': 1.0,
+ 'nonwetting': 1.0},
}
-# Dict of the form: { subdom_num : density }
+
densities = {
- 1: {'wetting': 1, # 997
- 'nonwetting': 1}, # 1.225}},
- 2: {'wetting': 1, # 997
- 'nonwetting': 1}, # 1.225}},
- 3: {'wetting': 1, # 997
- 'nonwetting': 1}, # 1.225}},
- 4: {'wetting': 1, # 997
- 'nonwetting': 1}, # 1.225}}
+ 1: {'wetting': 1.0, # 997
+ 'nonwetting': 1.0}, # 1.225}},
+ 2: {'wetting': 1.0, # 997
+ 'nonwetting': 1.0}, # 1.225}},
+ 3: {'wetting': 1.0, # 997
+ 'nonwetting': 1.0}, # 1.225}},
+ 4: {'wetting': 1.0, # 997
+ 'nonwetting': 1.0}, # 1.225}}
}
-
-gravity_acceleration = 1
+gravity_acceleration = 1.0
# porosities taken from
# https://www.geotechdata.info/parameter/soil-porosity.html
# Dict of the form: { subdom_num : porosity }
porosity = {
- 1: 1, # 0.2, # Clayey gravels, clayey sandy gravels
- 2: 1, # 0.22, # Silty gravels, silty sandy gravels
- 3: 1, # 0.37, # Clayey sands
- 4: 1, # 0.2 # Silty or sandy clay
+ 1: 1.0, #0.37, # 0.2, # Clayey gravels, clayey sandy gravels
+ 2: 1.0, #0.22, # 0.22, # Silty gravels, silty sandy gravels
+ 3: 1.0, #0.2, # 0.37, # Clayey sands
+ 4: 1.0, #0.22, # 0.2 # Silty or sandy clay
}
# subdom_num : subdomain L for L-scheme
L = {
- 1: {'wetting': Lw,
- 'nonwetting': Lnw},
- 2: {'wetting': Lw,
- 'nonwetting': Lnw},
- 3: {'wetting': Lw,
- 'nonwetting': Lnw},
- 4: {'wetting': Lw,
- 'nonwetting': Lnw}
+ 1: {'wetting': Lw1,
+ 'nonwetting': Lnw1},
+ 2: {'wetting': Lw2,
+ 'nonwetting': Lnw2},
+ 3: {'wetting': Lw3,
+ 'nonwetting': Lnw3},
+ 4: {'wetting': Lw4,
+ 'nonwetting': Lnw4}
}
-# subdom_num : lambda parameter for the L-scheme
# interface_num : lambda parameter for the L-scheme on that interface.
# Note that interfaces are numbered starting from 0, because
# adjacent_subdomains is a list and not a dict. Historic fuckup, I know
@@ -332,11 +366,13 @@ lambda_param = {
'nonwetting': lambda34_nw},
}
+
+# after Lewis, see pdf file
intrinsic_permeability = {
- 1: 1,
- 2: 1,
- 3: 1,
- 4: 1
+ 1: 1.0, #0.01, # sand
+ 2: 1.0, #0.01, # sand, there is a range
+ 3: 1.0, #0.01, #10e-2, # clay has a range
+ 4: 1.0, #0.01, #10e-3
}
@@ -354,19 +390,46 @@ def rel_perm1nw(s):
# relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
- return intrinsic_permeability[2]*s**3
+ return intrinsic_permeability[2]*s**2
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdomain2
- return intrinsic_permeability[2]*(1-s)**3
+ return intrinsic_permeability[2]*(1-s)**2
+
+
+# relative permeabilty functions on subdomain 3
+def rel_perm3w(s):
+ # relative permeabilty wetting on subdomain3
+ return intrinsic_permeability[3]*s**3
+
+
+def rel_perm3nw(s):
+ # relative permeabilty nonwetting on subdomain3
+ return intrinsic_permeability[3]*(1-s)**3
+# relative permeabilty functions on subdomain 4
+def rel_perm4w(s):
+ # relative permeabilty wetting on subdomain4
+ return intrinsic_permeability[4]*s**3
+
+
+def rel_perm4nw(s):
+ # relative permeabilty nonwetting on subdomain4
+ return intrinsic_permeability[4]*(1-s)**3
+
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)
+_rel_perm3w = ft.partial(rel_perm3w)
+_rel_perm3nw = ft.partial(rel_perm3nw)
+_rel_perm4w = ft.partial(rel_perm4w)
+_rel_perm4nw = ft.partial(rel_perm4nw)
+
+
subdomain1_rel_perm = {
'wetting': _rel_perm1w,
'nonwetting': _rel_perm1nw
@@ -377,18 +440,28 @@ subdomain2_rel_perm = {
'nonwetting': _rel_perm2nw
}
-# _rel_perm3 = ft.partial(rel_perm2)
-# subdomain3_rel_perm = subdomain2_rel_perm.copy()
-#
-# _rel_perm4 = ft.partial(rel_perm1)
-# subdomain4_rel_perm = subdomain1_rel_perm.copy()
+subdomain3_rel_perm = {
+ 'wetting': _rel_perm3w,
+ 'nonwetting': _rel_perm3nw
+}
+
+subdomain4_rel_perm = {
+ 'wetting': _rel_perm4w,
+ 'nonwetting': _rel_perm4nw
+}
# dictionary of relative permeabilties on all domains.
+# relative_permeability = {
+# 1: subdomain1_rel_perm,
+# 2: subdomain1_rel_perm,
+# 3: subdomain2_rel_perm,
+# 4: subdomain2_rel_perm
+# }
relative_permeability = {
1: subdomain1_rel_perm,
- 2: subdomain1_rel_perm,
- 3: subdomain2_rel_perm,
- 4: subdomain2_rel_perm
+ 2: subdomain2_rel_perm,
+ 3: subdomain3_rel_perm,
+ 4: subdomain4_rel_perm
}
@@ -404,22 +477,48 @@ def rel_perm1nw_prime(s):
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 subdomain1
- return intrinsic_permeability[2]*3*s**2
+ # relative permeabilty on subdomain2
+ return intrinsic_permeability[2]*2*s
def rel_perm2nw_prime(s):
- # relative permeabilty on subdomain1
- return -1*intrinsic_permeability[2]*3*(1-s)**2
+ # relative permeabilty on subdomain2
+ return -1*intrinsic_permeability[2]*2*(1-s)
+
+
+# definition of the derivatives of the relative permeabilities
+# relative permeabilty functions on subdomain 3
+def rel_perm3w_prime(s):
+ # relative permeabilty on subdomain3
+ return intrinsic_permeability[3]*3*s**2
+
+
+def rel_perm3nw_prime(s):
+ # relative permeabilty on subdomain3
+ return -1*intrinsic_permeability[3]*3*(1-s)**2
+
+
+# definition of the derivatives of the relative permeabilities
+# relative permeabilty functions on subdomain 4
+def rel_perm4w_prime(s):
+ # relative permeabilty on subdomain4
+ return intrinsic_permeability[4]*3*s**2
+
+
+def rel_perm4nw_prime(s):
+ # relative permeabilty on subdomain4
+ return -1*intrinsic_permeability[4]*3*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
+_rel_perm3w_prime = ft.partial(rel_perm3w_prime)
+_rel_perm3nw_prime = ft.partial(rel_perm3nw_prime)
+_rel_perm4w_prime = ft.partial(rel_perm4w_prime)
+_rel_perm4nw_prime = ft.partial(rel_perm4nw_prime)
subdomain1_rel_perm_prime = {
'wetting': _rel_perm1w_prime,
@@ -432,15 +531,32 @@ subdomain2_rel_perm_prime = {
'nonwetting': _rel_perm2nw_prime
}
+subdomain3_rel_perm_prime = {
+ 'wetting': _rel_perm3w_prime,
+ 'nonwetting': _rel_perm3nw_prime
+}
+
+
+subdomain4_rel_perm_prime = {
+ 'wetting': _rel_perm4w_prime,
+ 'nonwetting': _rel_perm4nw_prime
+}
+
+
# dictionary of relative permeabilties on all domains.
+# ka_prime = {
+# 1: subdomain1_rel_perm_prime,
+# 2: subdomain1_rel_perm_prime,
+# 3: subdomain2_rel_perm_prime,
+# 4: subdomain2_rel_perm_prime
+# }
ka_prime = {
1: subdomain1_rel_perm_prime,
- 2: subdomain1_rel_perm_prime,
- 3: subdomain2_rel_perm_prime,
- 4: subdomain2_rel_perm_prime
+ 2: subdomain2_rel_perm_prime,
+ 3: subdomain3_rel_perm_prime,
+ 4: subdomain4_rel_perm_prime
}
-
# S-pc-relation ship. We use the van Genuchten approach, i.e.
# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume m = 1-1/n
# (see Helmig) and assume that residual saturation is Sw
@@ -500,17 +616,18 @@ sat_pressure_relationship = {
}
-#############################################
+###############################################################################
# Manufacture source expressions with sympy #
-#############################################
+###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym_2patch = {
- 1: {'wetting': -7 - (1+t*t)*(1 + x*x + y*y),
+ 1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
'nonwetting': 0.0*t}, # -1-t*(1.1 + y + x**2)**2},
2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x),
- 'nonwetting': (-1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2))*y**2},
+ 'nonwetting': (-2-t*(1.1+y + x**2))*y**2},
+ # 'nonwetting': (-1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2))*y**2},
}
p_e_sym = {
@@ -583,6 +700,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
@@ -612,6 +730,7 @@ for subdomain in isRichards.keys():
{outer_boundary_ind: exact_solution[subdomain]}
)
+
# LOG FILE OUTPUT #############################################################
# read this file and print it to std out. This way the simulation can produce a
# log file with ./TP-R-layered_soil.py | tee simulation.log
diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
index 4392f91e5ab4cc23c377e40a0e28a9f0b5b9b066..49141a97ba720dc323f76c9f28e4a35df115f8c2 100755
--- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
+++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py
@@ -1,23 +1,24 @@
#!/usr/bin/python3
-"""Layered soil simulation.
+""" TP-R Layered soil simulation.
This program sets up an LDD simulation
"""
import dolfin as df
-# import mshr
-# import numpy as np
import sympy as sym
-# import typing as tp
import functools as ft
-# import domainPatch as dp
import LDDsimulation as ldd
import helpers as hlp
import datetime
import os
import pandas as pd
-# check if output directory exists
+# init sympy session
+sym.init_printing()
+
+# PREREQUISITS ###############################################################
+# check if output directory "./output" exists. This will be used in
+# the generation of the output string.
if not os.path.exists('./output'):
os.mkdir('./output')
print("Directory ", './output', " created ")
@@ -25,40 +26,43 @@ else:
print("Directory ", './output', " already exists. Will use as output \
directory")
-
date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
-# init sympy session
-sym.init_printing()
-# solver_tol = 6E-7
+
+# Name of the usecase that will be printed during simulation.
use_case = "TP-R-layered-soil-realistic"
-# name of this very file. Needed for log output.
+# The name of this very file. Needed for creating log output.
thisfile = "TP-R-layered_soil.py"
+
+# GENERAL SOLVER CONFIG ######################################################
+# maximal iteration per timestep
max_iter_num = 300
FEM_Lagrange_degree = 1
+
+# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False
resolutions = {
- # 1: 2e-6, # h=2
- # 2: 2e-6, # h=1.1180
- # 4: 2e-6, # h=0.5590
- # 8: 2e-6, # h=0.2814
- # 16: 2e-6, # h=0.1412
- 32: 2e-6,
- # 64: 2e-6,
- # 128: 2e-6
+ # 1: 1e-6,
+ # 2: 1e-6,
+ # 4: 1e-6,
+ # 8: 1e-6,
+ # 16: 5e-6,
+ # 32: 5e-6,
+ 64: 2e-6,
+ # 128: 1e-6,
+ # 256: 1e-6,
}
-# GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.001
-number_of_timesteps = 1000
-plot_timestep_every = 4
-# decide how many timesteps you want analysed. Analysed means, that we write
-# out subsequent errors of the L-iteration within the timestep.
-number_of_timesteps_to_analyse = 4
+# starttimes gives a list of starttimes to run the simulation from.
+# The list is looped over and a simulation is run with t_0 as initial time
+# for each element t_0 in starttimes.
starttimes = [0.0]
+timestep_size = 0.001
+number_of_timesteps = 20
+
+# LDD scheme parameters ######################################################
Lw1 = 0.025 # /timestep_size
Lnw1 = Lw1
Lw2 = 0.025 # /timestep_size
@@ -79,27 +83,38 @@ include_gravity = False
debugflag = False
analyse_condition = True
-if mesh_study:
- output_string = "./output/{}-{}_timesteps{}_P{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
- )
-else:
- for tol in resolutions.values():
- solver_tol = tol
- output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
- datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
- )
-
+# I/O CONFIG #################################################################
+# when number_of_timesteps is high, it might take a long time to write all
+# timesteps to disk. Therefore, you can choose to only write data of every
+# plot_timestep_every timestep to disk.
+plot_timestep_every = 4
+# Decide how many timesteps you want analysed. Analysed means, that
+# subsequent errors of the L-iteration within the timestep are written out.
+number_of_timesteps_to_analyse = 5
-# toggle what should be written to files
+# fine grained control over data to be written to disk in the mesh study case
+# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
+ # output the relative errornorm (integration in space) w.r.t. an exact
+ # solution for each timestep into a csv file.
'space_errornorms': True,
+ # save the mesh and marker functions to disk
'meshes_and_markers': True,
+ # save xdmf/h5 data for each LDD iteration for timesteps determined by
+ # number_of_timesteps_to_analyse. I/O intensive!
'L_iterations_per_timestep': False,
+ # save solution to xdmf/h5.
'solutions': True,
+ # save absolute differences w.r.t an exact solution to xdmf/h5 file
+ # to monitor where on the domains errors happen
'absolute_differences': True,
+ # analyise condition numbers for timesteps determined by
+ # number_of_timesteps_to_analyse and save them over time to csv.
'condition_numbers': analyse_condition,
+ # output subsequent iteration errors measured in L^2 to csv for
+ # timesteps determined by number_of_timesteps_to_analyse.
+ # Usefull to monitor convergence of the acutal LDD solver.
'subsequent_errors': True
}
else:
@@ -113,7 +128,20 @@ else:
'subsequent_errors': True
}
+# OUTPUT FILE STRING #########################################################
+if mesh_study:
+ output_string = "./output/{}-{}_timesteps{}_P{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
+ )
+else:
+ for tol in resolutions.values():
+ solver_tol = tol
+ output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
+ datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
+ )
+
+# DOMAIN AND INTERFACE #######################################################
# global domain
subdomain0_vertices = [
df.Point(-1.0, -1.0),
@@ -263,6 +291,8 @@ outer_boundary_def_points = {
4: subdomain4_outer_boundary_verts
}
+
+# MODEL CONFIGURATION #########################################################
isRichards = {
1: True,
2: True,
@@ -346,23 +376,14 @@ lambda_param = {
2: {'wetting': lambda34_w,
'nonwetting': lambda34_nw},
}
-# lambda_param = {
-# 1: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 2: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 3: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# 4: {'wetting': lambda_w,
-# 'nonwetting': lambda_nw},
-# }
+
# after Lewis, see pdf file
intrinsic_permeability = {
- 1: 0.1, # sand
- 2: 0.1, # sand, there is a range
- 3: 0.001, #10e-2, # clay has a range
- 4: 0.001, #10e-3
+ 1: 0.01, # sand
+ 2: 0.01, # sand, there is a range
+ 3: 0.01, #10e-2, # clay has a range
+ 4: 0.01, #10e-3
}
@@ -606,17 +627,17 @@ sat_pressure_relationship = {
}
-#############################################
+###############################################################################
# Manufacture source expressions with sympy #
-#############################################
+###############################################################################
x, y = sym.symbols('x[0], x[1]') # needed by UFL
t = sym.symbols('t', positive=True)
p_e_sym_2patch = {
- 1: {'wetting': -6 - (1+t*t)*(1 + x*x + y*y),
+ 1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),
'nonwetting': 0.0*t}, # -1-t*(1.1 + y + x**2)**2},
- 2: {'wetting': -6.0 - (1.0 + t*t)*(1.0 + x*x),
- 'nonwetting': (-1-t*(1.1+y + x**2))*y**2},
+ 2: {'wetting': -7.0 - (1.0 + t*t)*(1.0 + x*x),
+ 'nonwetting': (-2-t*(1.1+y + x**2))*y**2},
# 'nonwetting': (-1-t*(1.1 + x**2)**2 - sym.sqrt(5+t**2))*y**2},
}
@@ -690,6 +711,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
index 62a0b3c8479b4f4dd87e0074eedc8395e78d7f7a..5c439052ea5a8235fd6ed909ae5f74426c57c088 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py
@@ -464,6 +464,7 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
+# BOUNDARY CONDITIONS #########################################################
# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
@@ -478,7 +479,6 @@ dirichletBC = dict()
# return the actual expression needed for the dirichlet condition for both
# phases if present.
-# BOUNDARY CONDITIONS #########################################################
# subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys():
# subdomain can have no outer boundary