Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
L
LDD-for-two-phase-flow-systems
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Container registry
Model registry
Operate
Environments
Monitor
Incidents
Service Desk
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
David Seus
LDD-for-two-phase-flow-systems
Commits
5e274370
Commit
5e274370
authored
5 years ago
by
David Seus
Browse files
Options
Downloads
Patches
Plain Diff
add parallelisationtest
parent
845567da
No related branches found
No related tags found
No related merge requests found
Changes
1
Hide whitespace changes
Inline
Side-by-side
Showing
1 changed file
Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic_parallel_test.py
+1134
-0
1134 additions, 0 deletions
...-layered_soil_with_inner_patch-realistic_parallel_test.py
with
1134 additions
and
0 deletions
Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic_parallel_test.py
0 → 100755
+
1134
−
0
View file @
5e274370
#!/usr/bin/python3
"""
Layered soil simulation with inner patch.
This program sets up an LDD simulation
"""
import
dolfin
as
df
import
sympy
as
sym
import
functools
as
ft
import
LDDsimulation
as
ldd
import
helpers
as
hlp
import
datetime
import
os
import
pandas
as
pd
import
multiprocessing
as
mp
# 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
"
)
else
:
print
(
"
Directory
"
,
'
./output
'
,
"
already exists. Will use as output
\
directory
"
)
date
=
datetime
.
datetime
.
now
()
datestr
=
date
.
strftime
(
"
%Y-%m-%d
"
)
# Name of the usecase that will be printed during simulation.
use_case
=
"
TP-R-layered-soil-with-inner-patch-realistic-split-interface45
"
# The name of this very file. Needed for creating log output.
thisfile
=
"
TP-R-layered_soil_with_inner_patch-realistic.py
"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
max_iter_num
=
1
FEM_Lagrange_degree
=
1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study
=
True
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: 8e-6, # h=0.1412
# 32: 5e-6,
# 64: 2e-6,
# 128: 2e-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
=
1
# LDD scheme parameters ######################################################
Lw1
=
0.5
# /timestep_size
Lnw1
=
Lw1
Lw2
=
0.5
# /timestep_size
Lnw2
=
Lw2
Lw3
=
0.5
# /timestep_size
Lnw3
=
Lw3
Lw4
=
0.5
# /timestep_size
Lnw4
=
Lw4
Lw5
=
0.5
# /timestep_size
Lnw5
=
Lw5
Lw6
=
0.5
# /timestep_size
Lnw6
=
Lw6
lambda12_w
=
4
lambda12_nw
=
4
lambda23_w
=
4
lambda23_nw
=
4
lambda24_w
=
4
lambda24_nw
=
4
lambda25_w
=
4
lambda25_nw
=
4
lambda34_w
=
4
lambda34_nw
=
4
lambda36_w
=
4
lambda36_nw
=
4
lambda45_w
=
4
lambda45_nw
=
4
lambda46_w
=
4
lambda46_nw
=
4
lambda56_w
=
4
lambda56_nw
=
4
include_gravity
=
True
debugflag
=
True
analyse_condition
=
False
# 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
=
1
# 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
=
1
# 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
:
write_to_file
=
{
'
space_errornorms
'
:
True
,
'
meshes_and_markers
'
:
True
,
'
L_iterations_per_timestep
'
:
False
,
'
solutions
'
:
True
,
'
absolute_differences
'
:
True
,
'
condition_numbers
'
:
analyse_condition
,
'
subsequent_errors
'
:
True
}
# OUTPUT FILE STRING #########################################################
output_string
=
"
./output/{}-{}_timesteps{}_P{}
"
.
format
(
datestr
,
use_case
,
number_of_timesteps
,
FEM_Lagrange_degree
)
# DOMAIN AND INTERFACES #######################################################
# global domain
subdomain0_vertices
=
[
df
.
Point
(
-
1.0
,
-
1.0
),
#
df
.
Point
(
1.0
,
-
1.0
),
#
df
.
Point
(
1.0
,
1.0
),
#
df
.
Point
(
-
1.0
,
1.0
)]
interface12_vertices
=
[
df
.
Point
(
-
1.0
,
0.8
),
df
.
Point
(
0.3
,
0.8
),
df
.
Point
(
0.5
,
0.9
),
df
.
Point
(
0.8
,
0.7
),
df
.
Point
(
1.0
,
0.65
)]
# interface23
interface23_vertices
=
[
df
.
Point
(
-
1.0
,
0.0
),
df
.
Point
(
-
0.35
,
0.0
),
# df.Point(6.5, 4.5),
df
.
Point
(
0.0
,
0.0
)]
interface24_vertices
=
[
interface23_vertices
[
2
],
df
.
Point
(
0.6
,
0.0
),
]
interface25_vertices
=
[
interface24_vertices
[
1
],
df
.
Point
(
1.0
,
0.0
)
]
interface32_vertices
=
[
interface23_vertices
[
2
],
interface23_vertices
[
1
],
interface23_vertices
[
0
]]
interface36_vertices
=
[
df
.
Point
(
-
1.0
,
-
0.6
),
df
.
Point
(
-
0.6
,
-
0.45
)]
interface46_vertices
=
[
interface36_vertices
[
1
],
df
.
Point
(
0.3
,
-
0.25
)]
interface56_vertices
=
[
interface46_vertices
[
1
],
df
.
Point
(
0.65
,
-
0.6
),
df
.
Point
(
1.0
,
-
0.7
)]
interface34_vertices
=
[
interface36_vertices
[
1
],
interface23_vertices
[
2
]]
# Interface 45 needs to be split, because of the shape. There can be triangles
# with two facets on the interface and this creates a rogue dof type error when
# integrating over that particular interface. Accordingly, the lambda_param
# dictionary has two entries for that interface.
interface45_vertices_a
=
[
interface56_vertices
[
0
],
df
.
Point
(
0.7
,
-
0.2
),
#df.Point(0.7, -0.2),
]
interface45_vertices_b
=
[
df
.
Point
(
0.7
,
-
0.2
),
#df.Point(0.7, -0.2),
interface25_vertices
[
0
]
]
# interface_vertices introduces a global numbering of interfaces.
interface_def_points
=
[
interface12_vertices
,
interface23_vertices
,
interface24_vertices
,
interface25_vertices
,
interface34_vertices
,
interface36_vertices
,
interface45_vertices_a
,
interface45_vertices_b
,
interface46_vertices
,
interface56_vertices
,
]
adjacent_subdomains
=
[[
1
,
2
],
[
2
,
3
],
[
2
,
4
],
[
2
,
5
],
[
3
,
4
],
[
3
,
6
],
[
4
,
5
],
[
4
,
5
],
[
4
,
6
],
[
5
,
6
]
]
# subdomain1.
subdomain1_vertices
=
[
interface12_vertices
[
0
],
interface12_vertices
[
1
],
interface12_vertices
[
2
],
interface12_vertices
[
3
],
interface12_vertices
[
4
],
# southern boundary, 12 interface
subdomain0_vertices
[
2
],
# eastern boundary, outer boundary
subdomain0_vertices
[
3
]]
# northern boundary, outer on_boundary
# vertex coordinates of the outer boundaries. If it can not be specified as a
# polygon, use an entry per boundary polygon. This information is used for defining
# the Dirichlet boundary conditions. If a domain is completely internal, the
# dictionary entry should be 0: None
subdomain1_outer_boundary_verts
=
{
0
:
[
subdomain1_vertices
[
4
],
#
subdomain1_vertices
[
5
],
# eastern boundary, outer boundary
subdomain1_vertices
[
6
],
subdomain1_vertices
[
0
]]
}
#subdomain1
subdomain2_vertices
=
[
interface23_vertices
[
0
],
interface23_vertices
[
1
],
interface23_vertices
[
2
],
interface24_vertices
[
1
],
interface25_vertices
[
1
],
# southern boundary, 23 interface
subdomain1_vertices
[
4
],
# eastern boundary, outer boundary
subdomain1_vertices
[
3
],
subdomain1_vertices
[
2
],
subdomain1_vertices
[
1
],
subdomain1_vertices
[
0
]
]
# northern boundary, 12 interface
subdomain2_outer_boundary_verts
=
{
0
:
[
subdomain2_vertices
[
9
],
subdomain2_vertices
[
0
]],
1
:
[
subdomain2_vertices
[
4
],
subdomain2_vertices
[
5
]]
}
subdomain3_vertices
=
[
interface36_vertices
[
0
],
interface36_vertices
[
1
],
# interface34_vertices[0],
interface34_vertices
[
1
],
# interface32_vertices[0],
interface32_vertices
[
1
],
interface32_vertices
[
2
]
]
subdomain3_outer_boundary_verts
=
{
0
:
[
subdomain3_vertices
[
4
],
subdomain3_vertices
[
0
]]
}
# subdomain3
subdomain4_vertices
=
[
interface46_vertices
[
0
],
interface46_vertices
[
1
],
# interface45_vertices[1],
interface45_vertices_a
[
1
],
interface24_vertices
[
1
],
interface24_vertices
[
0
],
interface34_vertices
[
1
]
]
subdomain4_outer_boundary_verts
=
None
subdomain5_vertices
=
[
interface56_vertices
[
0
],
interface56_vertices
[
1
],
interface56_vertices
[
2
],
interface25_vertices
[
1
],
interface25_vertices
[
0
],
interface45_vertices_b
[
1
],
interface45_vertices_b
[
0
]
]
subdomain5_outer_boundary_verts
=
{
0
:
[
subdomain5_vertices
[
2
],
subdomain5_vertices
[
3
]]
}
subdomain6_vertices
=
[
subdomain0_vertices
[
0
],
subdomain0_vertices
[
1
],
# southern boundary, outer boundary
interface56_vertices
[
2
],
interface56_vertices
[
1
],
interface56_vertices
[
0
],
interface36_vertices
[
1
],
interface36_vertices
[
0
]
]
subdomain6_outer_boundary_verts
=
{
0
:
[
subdomain6_vertices
[
6
],
subdomain6_vertices
[
0
],
subdomain6_vertices
[
1
],
subdomain6_vertices
[
2
]]
}
subdomain_def_points
=
[
subdomain0_vertices
,
#
subdomain1_vertices
,
#
subdomain2_vertices
,
#
subdomain3_vertices
,
#
subdomain4_vertices
,
subdomain5_vertices
,
subdomain6_vertices
]
# if a subdomain has no outer boundary write None instead, i.e.
# i: None
# if i is the index of the inner subdomain.
outer_boundary_def_points
=
{
# subdomain number
1
:
subdomain1_outer_boundary_verts
,
2
:
subdomain2_outer_boundary_verts
,
3
:
subdomain3_outer_boundary_verts
,
4
:
subdomain4_outer_boundary_verts
,
5
:
subdomain5_outer_boundary_verts
,
6
:
subdomain6_outer_boundary_verts
}
# MODEL CONFIGURATION #########################################################
isRichards
=
{
1
:
True
,
2
:
True
,
3
:
False
,
4
:
False
,
5
:
False
,
6
:
False
}
# isRichards = {
# 1: True,
# 2: True,
# 3: True,
# 4: True,
# 5: True,
# 6: True
# }
# Dict of the form: { subdom_num : viscosity }
viscosity
=
{
1
:
{
'
wetting
'
:
1
,
'
nonwetting
'
:
1
/
50
},
2
:
{
'
wetting
'
:
1
,
'
nonwetting
'
:
1
/
50
},
3
:
{
'
wetting
'
:
1
,
'
nonwetting
'
:
1
/
50
},
4
:
{
'
wetting
'
:
1
,
'
nonwetting
'
:
1
/
50
},
5
:
{
'
wetting
'
:
1
,
'
nonwetting
'
:
1
/
50
},
6
:
{
'
wetting
'
:
1
,
'
nonwetting
'
:
1
/
50
},
}
# Dict of the form: { subdom_num : density }
densities
=
{
1
:
{
'
wetting
'
:
997.0
,
#997
'
nonwetting
'
:
1.225
},
#1}, #1.225},
2
:
{
'
wetting
'
:
997.0
,
#997
'
nonwetting
'
:
1.225
},
#1.225},
3
:
{
'
wetting
'
:
997.0
,
#997
'
nonwetting
'
:
1.225
},
#1.225},
4
:
{
'
wetting
'
:
997.0
,
#997
'
nonwetting
'
:
1.225
},
#1.225}
5
:
{
'
wetting
'
:
997.0
,
#997
'
nonwetting
'
:
1.225
},
#1.225},
6
:
{
'
wetting
'
:
997.0
,
#997
'
nonwetting
'
:
1.225
}
#1.225}
}
gravity_acceleration
=
9.81
# porosities taken from
# https://www.geotechdata.info/parameter/soil-porosity.html
# Dict of the form: { subdom_num : porosity }
porosity
=
{
1
:
0.2
,
#0.2, # Clayey gravels, clayey sandy gravels
2
:
0.2
,
#0.22, # Silty gravels, silty sandy gravels
3
:
0.2
,
#0.37, # Clayey sands
4
:
0.2
,
#0.2 # Silty or sandy clay
5
:
0.2
,
#
6
:
0.2
,
#
}
# subdom_num : subdomain L for L-scheme
L
=
{
1
:
{
'
wetting
'
:
Lw1
,
'
nonwetting
'
:
Lnw1
},
2
:
{
'
wetting
'
:
Lw2
,
'
nonwetting
'
:
Lnw2
},
3
:
{
'
wetting
'
:
Lw3
,
'
nonwetting
'
:
Lnw3
},
4
:
{
'
wetting
'
:
Lw4
,
'
nonwetting
'
:
Lnw4
},
5
:
{
'
wetting
'
:
Lw5
,
'
nonwetting
'
:
Lnw5
},
6
:
{
'
wetting
'
:
Lw6
,
'
nonwetting
'
:
Lnw6
}
}
# 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
# We have defined above as interfaces
# # interface_vertices introduces a global numbering of interfaces.
# interface_def_points = [interface12_vertices,
# interface23_vertices,
# interface24_vertices,
# interface25_vertices,
# interface34_vertices,
# interface36_vertices,
# interface45_vertices_a,
# interface45_vertices_b,
# interface46_vertices,
# interface56_vertices,
# ]
lambda_param
=
{
0
:
{
'
wetting
'
:
lambda12_w
,
'
nonwetting
'
:
lambda12_nw
},
1
:
{
'
wetting
'
:
lambda23_w
,
'
nonwetting
'
:
lambda23_nw
},
2
:
{
'
wetting
'
:
lambda24_w
,
'
nonwetting
'
:
lambda24_nw
},
3
:
{
'
wetting
'
:
lambda25_w
,
'
nonwetting
'
:
lambda25_nw
},
4
:
{
'
wetting
'
:
lambda34_w
,
'
nonwetting
'
:
lambda34_nw
},
5
:
{
'
wetting
'
:
lambda36_w
,
'
nonwetting
'
:
lambda36_nw
},
6
:
{
'
wetting
'
:
lambda45_w
,
'
nonwetting
'
:
lambda45_nw
},
7
:
{
'
wetting
'
:
lambda45_w
,
'
nonwetting
'
:
lambda45_nw
},
8
:
{
'
wetting
'
:
lambda46_w
,
'
nonwetting
'
:
lambda46_nw
},
9
:
{
'
wetting
'
:
lambda56_w
,
'
nonwetting
'
:
lambda56_nw
},
}
# after Lewis, see pdf file
intrinsic_permeability
=
{
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
5
:
0.01
,
#10e-2, # clay has a range
6
:
0.01
,
#10e-3
}
# relative permeabilty functions on subdomain 1
def
rel_perm1w
(
s
):
# relative permeabilty wetting on subdomain1
return
intrinsic_permeability
[
1
]
*
s
**
2
def
rel_perm1nw
(
s
):
# relative permeabilty nonwetting on subdomain1
return
intrinsic_permeability
[
1
]
*
(
1
-
s
)
**
2
# relative permeabilty functions on subdomain 2
def
rel_perm2w
(
s
):
# relative permeabilty wetting on subdomain2
return
intrinsic_permeability
[
2
]
*
s
**
2
def
rel_perm2nw
(
s
):
# relative permeabilty nonwetting on subdomain2
return
intrinsic_permeability
[
2
]
*
(
1
-
s
)
**
2
# relative permeabilty functions on subdomain 3
def
rel_perm3w
(
s
):
# relative permeabilty wetting on subdomain3
return
intrinsic_permeability
[
3
]
*
s
**
3
def
rel_perm3nw
(
s
):
# relative permeabilty nonwetting on subdomain3
return
intrinsic_permeability
[
3
]
*
(
1
-
s
)
**
3
# relative permeabilty functions on subdomain 4
def
rel_perm4w
(
s
):
# relative permeabilty wetting on subdomain4
return
intrinsic_permeability
[
4
]
*
s
**
3
def
rel_perm4nw
(
s
):
# relative permeabilty nonwetting on subdomain4
return
intrinsic_permeability
[
4
]
*
(
1
-
s
)
**
3
# relative permeabilty functions on subdomain 5
def
rel_perm5w
(
s
):
# relative permeabilty wetting on subdomain5
return
intrinsic_permeability
[
5
]
*
s
**
3
def
rel_perm5nw
(
s
):
# relative permeabilty nonwetting on subdomain5
return
intrinsic_permeability
[
5
]
*
(
1
-
s
)
**
3
# relative permeabilty functions on subdomain 6
def
rel_perm6w
(
s
):
# relative permeabilty wetting on subdomain6
return
intrinsic_permeability
[
6
]
*
s
**
3
def
rel_perm6nw
(
s
):
# relative permeabilty nonwetting on subdomain6
return
intrinsic_permeability
[
6
]
*
(
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
)
_rel_perm5w
=
ft
.
partial
(
rel_perm5w
)
_rel_perm5nw
=
ft
.
partial
(
rel_perm5nw
)
_rel_perm6w
=
ft
.
partial
(
rel_perm6w
)
_rel_perm6nw
=
ft
.
partial
(
rel_perm6nw
)
subdomain1_rel_perm
=
{
'
wetting
'
:
_rel_perm1w
,
'
nonwetting
'
:
_rel_perm1nw
}
subdomain2_rel_perm
=
{
'
wetting
'
:
_rel_perm2w
,
'
nonwetting
'
:
_rel_perm2nw
}
subdomain3_rel_perm
=
{
'
wetting
'
:
_rel_perm3w
,
'
nonwetting
'
:
_rel_perm3nw
}
subdomain4_rel_perm
=
{
'
wetting
'
:
_rel_perm4w
,
'
nonwetting
'
:
_rel_perm4nw
}
subdomain5_rel_perm
=
{
'
wetting
'
:
_rel_perm5w
,
'
nonwetting
'
:
_rel_perm5nw
}
subdomain6_rel_perm
=
{
'
wetting
'
:
_rel_perm6w
,
'
nonwetting
'
:
_rel_perm6nw
}
# dictionary of relative permeabilties on all domains.
relative_permeability
=
{
1
:
subdomain1_rel_perm
,
2
:
subdomain2_rel_perm
,
3
:
subdomain3_rel_perm
,
4
:
subdomain4_rel_perm
,
5
:
subdomain5_rel_perm
,
6
:
subdomain6_rel_perm
}
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
def
rel_perm1w_prime
(
s
):
# relative permeabilty on subdomain1
return
intrinsic_permeability
[
1
]
*
2
*
s
def
rel_perm1nw_prime
(
s
):
# relative permeabilty on subdomain1
return
-
1
*
intrinsic_permeability
[
1
]
*
2
*
(
1
-
s
)
def
rel_perm2w_prime
(
s
):
# relative permeabilty on subdomain2
return
intrinsic_permeability
[
2
]
*
2
*
s
def
rel_perm2nw_prime
(
s
):
# relative permeabilty on subdomain2
return
-
1
*
intrinsic_permeability
[
2
]
*
2
*
(
1
-
s
)
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 3
def
rel_perm3w_prime
(
s
):
# relative permeabilty on subdomain3
return
intrinsic_permeability
[
3
]
*
3
*
s
**
2
def
rel_perm3nw_prime
(
s
):
# relative permeabilty on subdomain3
return
-
1
*
intrinsic_permeability
[
3
]
*
3
*
(
1
-
s
)
**
2
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 4
def
rel_perm4w_prime
(
s
):
# relative permeabilty on subdomain4
return
intrinsic_permeability
[
4
]
*
3
*
s
**
2
def
rel_perm4nw_prime
(
s
):
# relative permeabilty on subdomain4
return
-
1
*
intrinsic_permeability
[
4
]
*
3
*
(
1
-
s
)
**
2
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 5
def
rel_perm5w_prime
(
s
):
# relative permeabilty on subdomain5
return
intrinsic_permeability
[
5
]
*
3
*
s
**
2
def
rel_perm5nw_prime
(
s
):
# relative permeabilty on subdomain5
return
-
1
*
intrinsic_permeability
[
5
]
*
3
*
(
1
-
s
)
**
2
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 6
def
rel_perm6w_prime
(
s
):
# relative permeabilty on subdomain6
return
intrinsic_permeability
[
6
]
*
3
*
s
**
2
def
rel_perm6nw_prime
(
s
):
# relative permeabilty on subdomain6
return
-
1
*
intrinsic_permeability
[
6
]
*
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
)
_rel_perm5w_prime
=
ft
.
partial
(
rel_perm5w_prime
)
_rel_perm5nw_prime
=
ft
.
partial
(
rel_perm5nw_prime
)
_rel_perm6w_prime
=
ft
.
partial
(
rel_perm6w_prime
)
_rel_perm6nw_prime
=
ft
.
partial
(
rel_perm6nw_prime
)
subdomain1_rel_perm_prime
=
{
'
wetting
'
:
_rel_perm1w_prime
,
'
nonwetting
'
:
_rel_perm1nw_prime
}
subdomain2_rel_perm_prime
=
{
'
wetting
'
:
_rel_perm2w_prime
,
'
nonwetting
'
:
_rel_perm2nw_prime
}
subdomain3_rel_perm_prime
=
{
'
wetting
'
:
_rel_perm3w_prime
,
'
nonwetting
'
:
_rel_perm3nw_prime
}
subdomain4_rel_perm_prime
=
{
'
wetting
'
:
_rel_perm4w_prime
,
'
nonwetting
'
:
_rel_perm4nw_prime
}
subdomain5_rel_perm_prime
=
{
'
wetting
'
:
_rel_perm5w_prime
,
'
nonwetting
'
:
_rel_perm5nw_prime
}
subdomain6_rel_perm_prime
=
{
'
wetting
'
:
_rel_perm6w_prime
,
'
nonwetting
'
:
_rel_perm6nw_prime
}
# dictionary of relative permeabilties on all domains.
ka_prime
=
{
1
:
subdomain1_rel_perm_prime
,
2
:
subdomain2_rel_perm_prime
,
3
:
subdomain3_rel_perm_prime
,
4
:
subdomain4_rel_perm_prime
,
5
:
subdomain5_rel_perm_prime
,
6
:
subdomain6_rel_perm_prime
,
}
# S-pc-relation ship. We use the van Genuchten approach, i.e.
# pc = 1/alpha*(S^{-1/m} -1)^1/n, where we set alpha = 0, assume
# m = 1-1/n (see Helmig) and assume that residual saturation is Sw
# this function needs to be monotonically decreasing in the capillary pressure
# pc.
# Since in the richards case pc=-pw, this becomes as a function of pw a mono
# tonically INCREASING function like in our Richards-Richards paper. However
# since we unify the treatment in the code for Richards and two-phase, we need
# the same requierment
# for both cases, two-phase and Richards.
# def saturation(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# 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 saturation_sym(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# #df.conditional(pc > 0,
# return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
#
#
# # derivative of S-pc relationship with respect to pc. This is needed for the
# # construction of a analytic solution.
# def saturation_sym_prime(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# 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(saturation_sym, n_index=3, alpha=0.001),
# 2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
# 3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
# 4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
# 5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
# 6: ft.partial(saturation_sym, n_index=3, alpha=0.001)
# }
#
# S_pc_sym_prime = {
# 1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# 2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# 3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# 4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# 5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
# 6: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
# }
#
# sat_pressure_relationship = {
# 1: ft.partial(saturation, n_index=3, alpha=0.001),
# 2: ft.partial(saturation, n_index=3, alpha=0.001),
# 3: ft.partial(saturation, n_index=3, alpha=0.001),
# 4: ft.partial(saturation, n_index=3, alpha=0.001),
# 5: ft.partial(saturation, n_index=3, alpha=0.001),
# 6: ft.partial(saturation, n_index=3, alpha=0.001)
# }
def
saturation
(
pc
,
n_index
):
# inverse capillary pressure-saturation-relationship
return
df
.
conditional
(
pc
>
0
,
1
/
((
1
+
pc
)
**
(
1
/
(
n_index
+
1
))),
1
)
def
saturation_sym
(
pc
,
n_index
):
# inverse capillary pressure-saturation-relationship
return
1
/
((
1
+
pc
)
**
(
1
/
(
n_index
+
1
)))
def
saturation_sym_prime
(
pc
,
n_index
):
# inverse capillary pressure-saturation-relationship
return
-
1
/
((
n_index
+
1
)
*
(
1
+
pc
)
**
((
n_index
+
2
)
/
(
n_index
+
1
)))
S_pc_sym
=
{
1
:
ft
.
partial
(
saturation_sym
,
n_index
=
1
),
2
:
ft
.
partial
(
saturation_sym
,
n_index
=
1
),
3
:
ft
.
partial
(
saturation_sym
,
n_index
=
2
),
4
:
ft
.
partial
(
saturation_sym
,
n_index
=
2
),
5
:
ft
.
partial
(
saturation_sym
,
n_index
=
2
),
6
:
ft
.
partial
(
saturation_sym
,
n_index
=
2
)
}
S_pc_sym_prime
=
{
1
:
ft
.
partial
(
saturation_sym_prime
,
n_index
=
1
),
2
:
ft
.
partial
(
saturation_sym_prime
,
n_index
=
1
),
3
:
ft
.
partial
(
saturation_sym_prime
,
n_index
=
2
),
4
:
ft
.
partial
(
saturation_sym_prime
,
n_index
=
2
),
5
:
ft
.
partial
(
saturation_sym_prime
,
n_index
=
2
),
6
:
ft
.
partial
(
saturation_sym_prime
,
n_index
=
2
)
}
sat_pressure_relationship
=
{
1
:
ft
.
partial
(
saturation
,
n_index
=
1
),
2
:
ft
.
partial
(
saturation
,
n_index
=
1
),
3
:
ft
.
partial
(
saturation
,
n_index
=
2
),
4
:
ft
.
partial
(
saturation
,
n_index
=
2
),
5
:
ft
.
partial
(
saturation
,
n_index
=
2
),
6
:
ft
.
partial
(
saturation
,
n_index
=
2
)
}
###############################################################################
# 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
=
{
1
:
{
'
wetting
'
:
-
7.0
-
(
1.0
+
t
*
t
)
*
(
1.0
+
x
*
x
+
y
*
y
),
'
nonwetting
'
:
0
*
t
},
2
:
{
'
wetting
'
:
-
7.0
-
(
1.0
+
t
*
t
)
*
(
1.0
+
x
*
x
+
y
*
y
),
'
nonwetting
'
:
0
*
t
},
3
:
{
'
wetting
'
:
(
-
7.0
-
(
1.0
+
t
*
t
)
*
(
1.0
+
x
*
x
)),
'
nonwetting
'
:
(
-
1.0
-
t
*
(
1.0
+
x
**
2
)
-
sym
.
sqrt
(
2
+
t
**
2
)
**
2
)
*
y
**
2
},
4
:
{
'
wetting
'
:
(
-
7.0
-
(
1.0
+
t
*
t
)
*
(
1.0
+
x
*
x
)),
'
nonwetting
'
:
(
-
1.0
-
t
*
(
1.0
+
x
**
2
)
-
sym
.
sqrt
(
2
+
t
**
2
)
**
2
)
*
y
**
2
},
5
:
{
'
wetting
'
:
(
-
7.0
-
(
1.0
+
t
*
t
)
*
(
1.0
+
x
*
x
)),
'
nonwetting
'
:
(
-
1.0
-
t
*
(
1.0
+
x
**
2
)
-
sym
.
sqrt
(
2
+
t
**
2
)
**
2
)
*
y
**
2
},
6
:
{
'
wetting
'
:
(
-
7.0
-
(
1.0
+
t
*
t
)
*
(
1.0
+
x
*
x
)),
'
nonwetting
'
:
(
-
1.0
-
t
*
(
1.0
+
x
**
2
)
-
sym
.
sqrt
(
2
+
t
**
2
)
**
2
)
*
y
**
2
},
}
pc_e_sym
=
dict
()
for
subdomain
,
isR
in
isRichards
.
items
():
if
isR
:
pc_e_sym
.
update
({
subdomain
:
-
p_e_sym
[
subdomain
][
'
wetting
'
].
copy
()})
else
:
pc_e_sym
.
update
({
subdomain
:
p_e_sym
[
subdomain
][
'
nonwetting
'
].
copy
()
-
p_e_sym
[
subdomain
][
'
wetting
'
].
copy
()})
symbols
=
{
"
x
"
:
x
,
"
y
"
:
y
,
"
t
"
:
t
}
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
exact_solution_example
=
hlp
.
generate_exact_solution_expressions
(
symbols
=
symbols
,
isRichards
=
isRichards
,
symbolic_pressure
=
p_e_sym
,
symbolic_capillary_pressure
=
pc_e_sym
,
saturation_pressure_relationship
=
S_pc_sym
,
saturation_pressure_relationship_prime
=
S_pc_sym_prime
,
viscosity
=
viscosity
,
porosity
=
porosity
,
relative_permeability
=
relative_permeability
,
relative_permeability_prime
=
ka_prime
,
densities
=
densities
,
gravity_acceleration
=
gravity_acceleration
,
include_gravity
=
include_gravity
,
)
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
# None should be written instead of an expression.
# This is a bit of a brainfuck:
# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
# Since a domain patch can have several disjoint outer boundary parts, the
# expressions need to get an enumaration index which starts at 0.
# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of
# subdomain ind and boundary part j.
# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting']
# return the actual expression needed for the dirichlet condition for both
# phases if present.
# subdomain index: {outer boudary part index: {phase: expression}}
for
subdomain
in
isRichards
.
keys
():
# subdomain can have no outer boundary
if
outer_boundary_def_points
[
subdomain
]
is
None
:
dirichletBC
.
update
({
subdomain
:
None
})
else
:
dirichletBC
.
update
({
subdomain
:
dict
()})
# set the dirichlet conditions to be the same code as exact solution on
# the subdomain.
for
outer_boundary_ind
in
outer_boundary_def_points
[
subdomain
].
keys
():
dirichletBC
[
subdomain
].
update
(
{
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
# f = open(thisfile, 'r')
# print(f.read())
# f.close()
# RUN #########################################################################
def
run_simulation
(
starttime
:
float
=
0.0
,
mesh_resolution
:
float
=
32
,
**
kwargs
)
->
None
:
"""
Setup and run LDD simulation
Setup and run LDD simulation for starttime starttime and mesh resolution
mesh_resolution.
PARAMETERS
starttime #type float number >= 0 specifying the starting point
of the simulation.
mesh_resolution #type float mesh resolution parameter for Fenics.
"""
for
key
,
value
in
kwargs
.
items
():
print
(
"
%s == %s
"
%
(
key
,
value
))
parameter
=
kwargs
# initialise LDD simulation class
simulation
=
ldd
.
LDDsimulation
(
# tol=1E-14, #parameter["tol"],
LDDsolver_tol
=
parameter
[
'
solver_tol
'
],
debug
=
parameter
[
"
debugflag
"
],
max_iter_num
=
parameter
[
"
max_iter_num
"
],
FEM_Lagrange_degree
=
parameter
[
"
FEM_Lagrange_degree
"
],
mesh_study
=
parameter
[
"
mesh_study
"
]
)
simulation
.
set_parameters
(
use_case
=
parameter
[
"
use_case
"
],
output_dir
=
parameter
[
"
output_string
"
],
subdomain_def_points
=
parameter
[
"
subdomain_def_points
"
],
isRichards
=
parameter
[
"
isRichards
"
],
interface_def_points
=
parameter
[
"
interface_def_points
"
],
outer_boundary_def_points
=
parameter
[
"
outer_boundary_def_points
"
],
adjacent_subdomains
=
parameter
[
"
adjacent_subdomains
"
],
mesh_resolution
=
parameter
[
"
mesh_resolution
"
],
viscosity
=
parameter
[
"
viscosity
"
],
porosity
=
parameter
[
"
porosity
"
],
L
=
parameter
[
"
L
"
],
lambda_param
=
parameter
[
"
lambda_param
"
],
relative_permeability
=
parameter
[
"
relative_permeability
"
],
saturation
=
parameter
[
"
sat_pressure_relationship
"
],
starttime
=
parameter
[
"
starttime
"
],
number_of_timesteps
=
parameter
[
"
number_of_timesteps
"
],
number_of_timesteps_to_analyse
=
parameter
[
"
number_of_timesteps_to_analyse
"
],
plot_timestep_every
=
parameter
[
"
plot_timestep_every
"
],
timestep_size
=
parameter
[
"
timestep_size
"
],
sources
=
parameter
[
"
source_expression
"
],
initial_conditions
=
parameter
[
"
initial_condition
"
],
dirichletBC_expression_strings
=
parameter
[
"
dirichletBC
"
],
exact_solution
=
parameter
[
"
exact_solution
"
],
densities
=
parameter
[
"
densities
"
],
include_gravity
=
parameter
[
"
include_gravity
"
],
gravity_acceleration
=
parameter
[
"
gravity_acceleration
"
],
write2file
=
parameter
[
"
write_to_file
"
],
)
simulation
.
initialise
()
output_dir
=
simulation
.
output_dir
# simulation.write_exact_solution_to_xdmf()
output
=
simulation
.
run
(
analyse_condition
=
parameter
[
"
analyse_condition
"
])
for
subdomain_index
,
subdomain_output
in
output
.
items
():
mesh_h
=
subdomain_output
[
'
mesh_size
'
]
for
phase
,
error_dict
in
subdomain_output
[
'
errornorm
'
].
items
():
filename
=
output_dir
\
+
"
subdomain{}
"
.
format
(
subdomain_index
)
\
+
"
-space-time-errornorm-{}-phase_meshrez{}.csv
"
.
format
(
phase
,
mesh_resolution
)
# for errortype, errornorm in error_dict.items():
# eocfile = open("eoc_filename", "a")
# eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
# eocfile.close()
# if subdomain.isRichards:mesh_h
data_dict
=
{
'
mesh_parameter
'
:
mesh_resolution
,
'
mesh_h
'
:
mesh_h
,
}
for
norm_type
,
errornorm
in
error_dict
.
items
():
data_dict
.
update
(
{
norm_type
:
errornorm
}
)
errors
=
pd
.
DataFrame
(
data_dict
,
index
=
[
mesh_resolution
])
# check if file exists
if
os
.
path
.
isfile
(
filename
)
is
True
:
with
open
(
filename
,
'
a
'
)
as
f
:
errors
.
to_csv
(
f
,
header
=
False
,
sep
=
'
\t
'
,
encoding
=
'
utf-8
'
,
index
=
False
)
else
:
errors
.
to_csv
(
filename
,
sep
=
'
\t
'
,
encoding
=
'
utf-8
'
,
index
=
False
)
def
info
(
title
):
"""
show process ids of multiprocessing simulation
"""
print
(
title
)
print
(
'
module name:
'
,
__name__
)
print
(
'
parent process:
'
,
os
.
getppid
())
print
(
'
process id:
'
,
os
.
getpid
())
# MAIN ########################################################################
if
__name__
==
'
__main__
'
:
# mp.set_start_method('spawn')
# dictionary of simualation parameters to pass to the run function.
# mesh_resolution and starttime are excluded, as they get passed explicitly
# to achieve parallelisation in these parameters in these parameters for
# mesh studies etc.
simulation_parameter
=
{
"
tol
"
:
1E-14
,
"
debugflag
"
:
debugflag
,
"
max_iter_num
"
:
max_iter_num
,
"
FEM_Lagrange_degree
"
:
FEM_Lagrange_degree
,
"
mesh_study
"
:
mesh_study
,
"
use_case
"
:
use_case
,
"
output_string
"
:
output_string
,
"
subdomain_def_points
"
:
subdomain_def_points
,
"
isRichards
"
:
isRichards
,
"
interface_def_points
"
:
interface_def_points
,
"
outer_boundary_def_points
"
:
outer_boundary_def_points
,
"
adjacent_subdomains
"
:
adjacent_subdomains
,
# "mesh_resolution": mesh_resolution,
"
viscosity
"
:
viscosity
,
"
porosity
"
:
porosity
,
"
L
"
:
L
,
"
lambda_param
"
:
lambda_param
,
"
relative_permeability
"
:
relative_permeability
,
"
sat_pressure_relationship
"
:
sat_pressure_relationship
,
# "starttime": starttime,
"
number_of_timesteps
"
:
number_of_timesteps
,
"
number_of_timesteps_to_analyse
"
:
number_of_timesteps_to_analyse
,
"
plot_timestep_every
"
:
plot_timestep_every
,
"
timestep_size
"
:
timestep_size
,
"
source_expression
"
:
source_expression
,
"
initial_condition
"
:
initial_condition
,
"
dirichletBC
"
:
dirichletBC
,
"
exact_solution
"
:
exact_solution
,
"
densities
"
:
densities
,
"
include_gravity
"
:
include_gravity
,
"
gravity_acceleration
"
:
gravity_acceleration
,
"
write_to_file
"
:
write_to_file
,
"
analyse_condition
"
:
analyse_condition
}
for
starttime
in
starttimes
:
for
mesh_resolution
,
solver_tol
in
resolutions
.
items
():
simulation_parameter
.
update
(
{
"
solver_tol
"
:
solver_tol
,
"
mesh_resolution
"
:
mesh_resolution
,
"
starttime
"
:
starttime
}
)
# info(simulation_parameter["use_case"])
# LDDsim = mp.Process(
# target=run_simulation,
# args=(
# starttime,
# mesh_resolution,
# simulation_parameter
# )
# )
# LDDsim.start()
# LDDsim.join()
run_simulation
(
starttime
=
starttime
,
mesh_resolution
=
mesh_resolution
,
kwargs
=
simulation_parameter
)
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment