From 61e9cd605580aa3ad34b9c5579b378d786ab5b50 Mon Sep 17 00:00:00 2001
From: David <forenkram@gmx.de>
Date: Tue, 30 Jun 2020 17:52:09 +0200
Subject: [PATCH] make it possible to give initial iteration a timestep number
---
.../layered_soil/TP-R-layered_soil.py | 2 +-
...ed_soil_with_inner_patch-all-params-one.py | 8 +-
...layered_soil_with_inner_patch-realistic.py | 8 +-
...h_inner_patch-all-params-one-mesh-study.py | 8 +-
...oil_with_inner_patch-pure-dd-mesh-study.py | 8 +-
...layered_soil_with_inner_patch-realistic.py | 81 +-
.../TP-TP-layered_soil_with_inner_patch.py | 803 +++---------------
...ayered_soil_with_inner_patch_mesh_study.py | 802 +++--------------
8 files changed, 298 insertions(+), 1422 deletions(-)
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 880eb3a..4c4e798 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
@@ -57,7 +57,7 @@ resolutions = {
# 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]
+starttimes = {0: 0.0}
timestep_size = 0.001
number_of_timesteps = 5
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py
index 55bc777..8998a94 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-all-params-one.py
@@ -55,7 +55,7 @@ resolutions = {
# 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]
+starttimes = {0: 0.0}
timestep_size = 0.01
number_of_timesteps = 5
@@ -407,7 +407,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -450,7 +449,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
index aeb8c73..e4aebcd 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -55,7 +55,7 @@ resolutions = {
# 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]
+starttimes = {0: 0.0}
timestep_size = 0.001
number_of_timesteps = 10
@@ -405,7 +405,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -448,7 +447,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py
index 7fc2778..f5c0cfc 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-all-params-one-mesh-study.py
@@ -55,7 +55,7 @@ resolutions = {
# 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]
+starttimes = {0: 0.0}
timestep_size = 0.01
number_of_timesteps = 5
@@ -407,7 +407,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -450,7 +449,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py
index 3286af7..890fa5c 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-pure-dd-mesh-study.py
@@ -54,7 +54,7 @@ resolutions = {
# 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]
+starttimes = {0: 0.0}
timestep_size = 0.01
number_of_timesteps = 5
@@ -404,7 +404,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -447,7 +446,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-realistic.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-realistic.py
index adc2ca6..9af0be5 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/mesh_study/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -35,85 +35,85 @@ thisfile = "TP-R-layered_soil_with_inner_patch-realistic.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
-max_iter_num = 700
+max_iter_num = 1000
FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
-mesh_study = False
+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: 2e-6,
- # 64: 2e-6,
- # 128: 2e-6
+ 1: 5e-5, # h=2
+ 2: 5e-5, # h=1.1180
+ 4: 5e-5, # h=0.5590
+ 8: 5e-5, # h=0.2814
+ 16: 3e-5, # h=0.1412
+ 32: 5e-6,
+ 64: 3e-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.01
-number_of_timesteps = 5
+starttimes = {0: 0.0}
+timestep_size = 0.0025
+number_of_timesteps = 400
# LDD scheme parameters ######################################################
-Lw1 = 0.025 # /timestep_size
-Lnw1 = Lw1
+Lw1 = 0.5 # /timestep_size
+Lnw1 = 0.5
-Lw2 = 0.025 # /timestep_size
-Lnw2 = Lw2
+Lw2 = 0.5 # /timestep_size
+Lnw2 = 0.5
-Lw3 = 0.025 # /timestep_size
-Lnw3 = Lw3
+Lw3 = 0.5 # /timestep_size
+Lnw3 = 0.5
-Lw4 = 0.025 # /timestep_size
-Lnw4 = Lw4
+Lw4 = 0.5 # /timestep_size
+Lnw4 = 0.5
-Lw5 = 0.025 # /timestep_size
-Lnw5 = Lw5
+Lw5 = 0.5 # /timestep_size
+Lnw5 = 0.5
-Lw6 = 0.025 # /timestep_size
-Lnw6 = Lw6
+Lw6 = 0.5 # /timestep_size
+Lnw6 = 0.5
-lambda12_w = 4
+lambda12_w = 40
lambda12_nw = 4
-lambda23_w = 4
+lambda23_w = 40
lambda23_nw = 4
-lambda24_w = 4
+lambda24_w = 40
lambda24_nw= 4
-lambda25_w= 4
+lambda25_w= 40
lambda25_nw= 4
-lambda34_w = 4
+lambda34_w = 40
lambda34_nw = 4
-lambda36_w = 4
+lambda36_w = 40
lambda36_nw = 4
-lambda45_w = 4
+lambda45_w = 40
lambda45_nw = 4
-lambda46_w = 4
+lambda46_w = 40
lambda46_nw = 4
-lambda56_w = 4
+lambda56_w = 40
lambda56_nw = 4
include_gravity = True
-debugflag = True
+debugflag = False
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 = 4
+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 = 5
@@ -159,7 +159,6 @@ output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
-
# DOMAIN AND INTERFACES #######################################################
substructuring = dss.layeredSoilInnerPatch()
interface_def_points = substructuring.interface_def_points
@@ -405,7 +404,6 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
# MAIN ########################################################################
if __name__ == '__main__':
# dictionary of simualation parameters to pass to the run function.
@@ -448,7 +446,10 @@ if __name__ == '__main__':
"write_to_file": write_to_file,
"analyse_condition": analyse_condition
}
- for starttime in starttimes:
+ for number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
@@ -461,7 +462,7 @@ if __name__ == '__main__':
)
)
LDDsim.start()
- LDDsim.join()
+ # LDDsim.join()
# hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
diff --git a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
index 7be6b7a..5891b65 100755
--- a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
+++ b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/TP-TP-layered_soil_with_inner_patch.py
@@ -3,16 +3,15 @@
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()
@@ -56,7 +55,7 @@ resolutions = {
# 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.3]
+starttimes = {0: 0.3}
timestep_size = 0.001
number_of_timesteps = 5
@@ -160,216 +159,12 @@ 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
-}
+# DOMAIN AND INTERFACES ######################################################
+substructuring = dss.layeredSoilInnerPatch()
+interface_def_points = substructuring.interface_def_points
+adjacent_subdomains = substructuring.adjacent_subdomains
+subdomain_def_points = substructuring.subdomain_def_points
+outer_boundary_def_points = substructuring.outer_boundary_def_points
# MODEL CONFIGURATION #########################################################
@@ -494,359 +289,40 @@ intrinsic_permeability = {
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
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 5: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 6: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-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)
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
+
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 1}},
+ 3: {"testSpc": {"index": 2}},
+ 4: {"testSpc": {"index": 2}},
+ 5: {"testSpc": {"index": 2}},
+ 6: {"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 #
@@ -870,15 +346,10 @@ p_e_sym = {
'nonwetting': (-1 -t*(1.0 + x**2) - sym.sin(2+t**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()})
-
+pc_e_sym = hlp.generate_exact_symbolic_pc(
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym
+ )
symbols = {"x": x,
"y": y,
@@ -894,6 +365,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,
@@ -905,34 +377,17 @@ 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]}
- )
-
+# Dictionary of dirichlet boundary conditions. If an exact solution case is
+# used, use the hlp.generate_exact_DirichletBC() method to generate the
+# Dirichlet Boundary conditions from the exact solution.
+dirichletBC = hlp.generate_exact_DirichletBC(
+ isRichards=isRichards,
+ outer_boundary_def_points=outer_boundary_def_points,
+ exact_solution=exact_solution
+ )
+# If no exact solution is provided you need to provide a dictionary of boundary
+# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see
+# the structure.
# LOG FILE OUTPUT #############################################################
# read this file and print it to std out. This way the simulation can produce a
@@ -941,89 +396,67 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
-# RUN #########################################################################
-for starttime in starttimes:
- for mesh_resolution, solver_tol in resolutions.items():
- # initialise LDD simulation class
- simulation = ldd.LDDsimulation(
- tol=1E-14,
- LDDsolver_tol=solver_tol,
- debug=debugflag,
- max_iter_num=max_iter_num,
- FEM_Lagrange_degree=FEM_Lagrange_degree,
- mesh_study=mesh_study
- )
-
- simulation.set_parameters(
- use_case=use_case,
- output_dir=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,
- saturation=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,
- sources=source_expression,
- initial_conditions=initial_condition,
- dirichletBC_expression_strings=dirichletBC,
- exact_solution=exact_solution,
- densities=densities,
- include_gravity=include_gravity,
- gravity_acceleration=gravity_acceleration,
- write2file=write_to_file,
- )
-
- simulation.initialise()
- output_dir = simulation.output_dir
- # simulation.write_exact_solution_to_xdmf()
- output = simulation.run(analyse_condition=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.csv".format(phase)
- # 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
+# MAIN ########################################################################
+if __name__ == '__main__':
+ # 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,
+ "intrinsic_permeability": intrinsic_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 number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
+ for mesh_resolution, solver_tol in resolutions.items():
+ simulation_parameter.update({"solver_tol": solver_tol})
+ hlp.info(simulation_parameter["use_case"])
+ LDDsim = mp.Process(
+ target=hlp.run_simulation,
+ args=(
+ simulation_parameter,
+ starttime,
+ mesh_resolution
)
- else:
- errors.to_csv(
- filename,
- sep='\t',
- encoding='utf-8',
- index=False
)
+ LDDsim.start()
+ # LDDsim.join()
+ # hlp.run_simulation(
+ # mesh_resolution=mesh_resolution,
+ # starttime=starttime,
+ # parameter=simulation_parameter
+ # )
diff --git a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/mesh_study/TP-TP-layered_soil_with_inner_patch_mesh_study.py b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/mesh_study/TP-TP-layered_soil_with_inner_patch_mesh_study.py
index 212c5f4..6426c3f 100755
--- a/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/mesh_study/TP-TP-layered_soil_with_inner_patch_mesh_study.py
+++ b/Two-phase-Two-phase/multi-patch/TP-TP-layered-soil-case-with-inner-patch/mesh_study/TP-TP-layered_soil_with_inner_patch_mesh_study.py
@@ -3,15 +3,15 @@
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()
@@ -56,7 +56,7 @@ resolutions = {
# 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.3]
+starttimes = {0: 0.3}
timestep_size = 0.001
number_of_timesteps = 5
@@ -160,216 +160,12 @@ 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
-}
+# DOMAIN AND INTERFACES ######################################################
+substructuring = dss.layeredSoilInnerPatch()
+interface_def_points = substructuring.interface_def_points
+adjacent_subdomains = substructuring.adjacent_subdomains
+subdomain_def_points = substructuring.subdomain_def_points
+outer_boundary_def_points = substructuring.outer_boundary_def_points
# MODEL CONFIGURATION #########################################################
@@ -494,359 +290,40 @@ intrinsic_permeability = {
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
+# relative permeabilties
+rel_perm_definition = {
+ 1: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 2: {"wetting": "Spow2",
+ "nonwetting": "oneMinusSpow2"},
+ 3: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 4: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 5: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
+ 6: {"wetting": "Spow3",
+ "nonwetting": "oneMinusSpow3"},
}
-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)
+rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
+relative_permeability = rel_perm_dict["ka"]
+ka_prime = rel_perm_dict["ka_prime"]
+
+# S-pc relation
+Spc_on_subdomains = {
+ 1: {"testSpc": {"index": 1}},
+ 2: {"testSpc": {"index": 1}},
+ 3: {"testSpc": {"index": 2}},
+ 4: {"testSpc": {"index": 2}},
+ 5: {"testSpc": {"index": 2}},
+ 6: {"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 #
@@ -870,15 +347,10 @@ p_e_sym = {
'nonwetting': (-1 -t*(1.0 + x**2) - sym.sin(2+t**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()})
-
+pc_e_sym = hlp.generate_exact_symbolic_pc(
+ isRichards=isRichards,
+ symbolic_pressure=p_e_sym
+ )
symbols = {"x": x,
"y": y,
@@ -894,6 +366,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,
@@ -905,34 +378,17 @@ 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]}
- )
-
+# Dictionary of dirichlet boundary conditions. If an exact solution case is
+# used, use the hlp.generate_exact_DirichletBC() method to generate the
+# Dirichlet Boundary conditions from the exact solution.
+dirichletBC = hlp.generate_exact_DirichletBC(
+ isRichards=isRichards,
+ outer_boundary_def_points=outer_boundary_def_points,
+ exact_solution=exact_solution
+ )
+# If no exact solution is provided you need to provide a dictionary of boundary
+# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see
+# the structure.
# LOG FILE OUTPUT #############################################################
# read this file and print it to std out. This way the simulation can produce a
@@ -941,89 +397,67 @@ f = open(thisfile, 'r')
print(f.read())
f.close()
-
-# RUN #########################################################################
-for starttime in starttimes:
- for mesh_resolution, solver_tol in resolutions.items():
- # initialise LDD simulation class
- simulation = ldd.LDDsimulation(
- tol=1E-14,
- LDDsolver_tol=solver_tol,
- debug=debugflag,
- max_iter_num=max_iter_num,
- FEM_Lagrange_degree=FEM_Lagrange_degree,
- mesh_study=mesh_study
- )
-
- simulation.set_parameters(
- use_case=use_case,
- output_dir=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,
- saturation=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,
- sources=source_expression,
- initial_conditions=initial_condition,
- dirichletBC_expression_strings=dirichletBC,
- exact_solution=exact_solution,
- densities=densities,
- include_gravity=include_gravity,
- gravity_acceleration=gravity_acceleration,
- write2file=write_to_file,
- )
-
- simulation.initialise()
- output_dir = simulation.output_dir
- # simulation.write_exact_solution_to_xdmf()
- output = simulation.run(analyse_condition=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.csv".format(phase)
- # 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
+# MAIN ########################################################################
+if __name__ == '__main__':
+ # 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,
+ "intrinsic_permeability": intrinsic_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 number_shift, starttime in starttimes.items():
+ simulation_parameter.update(
+ {"starttime_timestep_number_shift": number_shift}
+ )
+ for mesh_resolution, solver_tol in resolutions.items():
+ simulation_parameter.update({"solver_tol": solver_tol})
+ hlp.info(simulation_parameter["use_case"])
+ LDDsim = mp.Process(
+ target=hlp.run_simulation,
+ args=(
+ simulation_parameter,
+ starttime,
+ mesh_resolution
)
- else:
- errors.to_csv(
- filename,
- sep='\t',
- encoding='utf-8',
- index=False
)
+ LDDsim.start()
+ # LDDsim.join()
+ # hlp.run_simulation(
+ # mesh_resolution=mesh_resolution,
+ # starttime=starttime,
+ # parameter=simulation_parameter
+ # )
--
GitLab