Skip to content
Snippets Groups Projects
Commit 14f37316 authored by David Seus's avatar David Seus
Browse files

setup another TPTP 2 patch mesh study

parent c47694cc
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-different-intrinsic-perm.py
+ 7
− 7
View file @ 14f37316
......@@ -30,7 +30,7 @@ date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
# Name of the usecase that will be printed during simulation.
use_case = "TP-TP-2P-realistic-different-intrinsic-perm-smaller-L"
use_case = "TP-TP-2P-realistic-different-intrinsic-perm-params1a-higher-tol"
# The name of this very file. Needed for creating log output.
thisfile = "TP-TP-2-patch-different-intrinsic-perm.py"
......@@ -61,14 +61,14 @@ timestep_size = 0.001
number_of_timesteps = 1000
# LDD scheme parameters ######################################################
Lw1 = 0.01 #/timestep_size
Lnw1= 0.01
Lw1 = 0.1 #/timestep_size
Lnw1= 0.1
Lw2 = 0.001 #/timestep_size
Lnw2= 0.001
Lw2 = 0.1 #/timestep_size
Lnw2= 0.1
lambda_w = 2
lambda_nw = 2
lambda_w = 40
lambda_nw = 40
include_gravity = False
debugflag = False
analyse_condition = False
......
Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-same-intrinsic-perm.py
+ 14
− 13
View file @ 14f37316
......@@ -29,7 +29,7 @@ date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")
# Name of the usecase that will be printed during simulation.
use_case = "TP-TP-2P-realistic-same-intrinsic0.1-control-for-case160620"
use_case = "TP-TP-2P-realistic-same-intrinsic"
# The name of this very file. Needed for creating log output.
thisfile = "TP-TP-2-patch-same-intrinsic-perm.py"
......@@ -46,7 +46,7 @@ resolutions = {
# 4: 1e-6,
# 8: 1e-6,
# 16: 5e-6,
32: 3e-6,
32: 1e-5,
# 64: 2e-6,
# 128: 1e-6,
# 256: 1e-6,
......@@ -56,20 +56,21 @@ resolutions = {
# 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.0}
# starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9}
timestep_size = 0.001
number_of_timesteps = 1000
# LDD scheme parameters ######################################################
Lw1 = 0.25 #/timestep_size
Lnw1= 0.25
Lw1 = 0.5 #/timestep_size
Lnw1= 0.2
Lw2 = 0.25 #/timestep_size
Lnw2= 0.25
Lw2 = 0.5 #/timestep_size
Lnw2= 0.2
lambda_w = 4
lambda_nw = 4
lambda_w = 2
lambda_nw = 1
include_gravity = False
include_gravity = True
debugflag = False
analyse_condition = False
......@@ -80,7 +81,7 @@ analyse_condition = False
plot_timestep_every = 4
# Decide how many timesteps you want analysed. Analysed means, that
# subsequent errors of the L-iteration within the timestep are written out.
number_of_timesteps_to_analyse = 5
number_of_timesteps_to_analyse = 10
# 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.
......@@ -176,8 +177,8 @@ lambda_param = {#
}
intrinsic_permeability = {
1: 0.1,
2: 0.1,
1: 0.01,
2: 0.01,
}
# RELATIVE PEMRMEABILITIES
......@@ -339,7 +340,7 @@ if __name__ == '__main__':
# parameter=simulation_parameter
# )
LDDsim.join()
# LDDsim.join()
if mesh_study:
simulation_output_dir = processQueue.get()
hlp.merge_spacetime_errornorms(isRichards=isRichards,
......
Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py
+ 215
− 407
View file @ 14f37316
#!/usr/bin/python3
"""TP-TP 2 patch soil simulation.
This program sets up an LDD simulation
"""
import dolfin as df
import mshr
import numpy as np
import sympy as sym
import typing as tp
import domainPatch as dp
import functions as fts
import LDDsimulation as ldd
import functools as ft
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()
# 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")
#import ufl as ufl
# init sympy session
sym.init_printing()
# Name of the usecase that will be printed during simulation.
use_case = "TP-TP-2-patch-realistic"
# solver_tol = 6E-7
max_iter_num = 300
# The name of this very file. Needed for creating log output.
thisfile = "TP-TP-2-patch-mesh-study.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
max_iter_num = 400
FEM_Lagrange_degree = 1
mesh_study = False
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = True
resolutions = {
# 1: 1e-6,
# 2: 1e-6,
# 4: 1e-6,
# 8: 1e-6,
# 16: 1e-6,
32: 1e-6,
# 64: 1e-6,
# 128: 1e-6,
1: 1e-5,
2: 1e-5,
4: 1e-5,
8: 1e-5,
16: 5e-6,
32: 3e-6,
64: 1e-6,
128: 1e-6,
# 256: 1e-6,
}
############ GRID #######################
# mesh_resolution = 20
# 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.0}
timestep_size = 0.001
number_of_timesteps = 20
plot_timestep_every = 2
# decide how many timesteps you want analysed. Analysed means, that we write out
# subsequent errors of the L-iteration within the timestep.
number_of_timesteps_to_analyse = 6
starttimes = [0.0]
number_of_timesteps = 1000
Lw = 0.025 #/timestep_size
Lnw=Lw
# LDD scheme parameters ######################################################
Lw1 = 0.01 #/timestep_size
Lnw1= 0.01
lambda_w = 4
lambda_nw = 4
Lw2 = 0.001 #/timestep_size
Lnw2= 0.001
include_gravity = True
debugflag = True
analyse_condition = True
lambda_w = 2
lambda_nw = 2
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)
include_gravity = False
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 = 2
# 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 = 8
# toggle what should be written to files
# fine grained control over data to be written to disk in the mesh study case
# as well as for a regular simuation for a fixed grid.
if mesh_study:
write_to_file = {
# output the relative errornorm (integration in space) w.r.t. an exact
# solution for each timestep into a csv file.
'space_errornorms': True,
# save the mesh and marker functions to disk
'meshes_and_markers': True,
'L_iterations_per_timestep': 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:
......@@ -86,70 +117,19 @@ else:
'subsequent_errors': True
}
# OUTPUT FILE STRING #########################################################
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
# DOMAIN AND INTERFACE #######################################################
substructuring = dss.twoSoilLayers()
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
##### Domain and Interface ####
# global simulation domain domain
sub_domain0_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)]
# interface between subdomain1 and subdomain2
interface12_vertices = [df.Point(-1.0, 0.0),
df.Point(1.0, 0.0) ]
# subdomain1.
sub_domain1_vertices = [interface12_vertices[0],
interface12_vertices[1],
sub_domain0_vertices[2],
sub_domain0_vertices[3]]
# 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: [interface12_vertices[1], #
sub_domain0_vertices[2],
sub_domain0_vertices[3], #
interface12_vertices[0]]
}
# subdomain2
sub_domain2_vertices = [sub_domain0_vertices[0],
sub_domain0_vertices[1],
interface12_vertices[1],
interface12_vertices[0] ]
subdomain2_outer_boundary_verts = {
0: [interface12_vertices[0], #
sub_domain0_vertices[0],
sub_domain0_vertices[1],
interface12_vertices[1]]
}
# list of subdomains given by the boundary polygon vertices.
# Subdomains are given as a list of dolfin points forming
# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
# to create the subdomain. subdomain_def_points[0] contains the
# vertices of the global simulation domain and subdomain_def_points[i] contains the
# vertices of the subdomain i.
subdomain_def_points = [sub_domain0_vertices,#
sub_domain1_vertices,#
sub_domain2_vertices]
# in the below list, index 0 corresponds to the 12 interface which has index 1
interface_def_points = [interface12_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
}
# adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]]
# MODEL CONFIGURATION #########################################################
isRichards = {
1: False, #
2: False
......@@ -167,7 +147,7 @@ viscosity = {#
porosity = {#
# subdom_num : porosity
1: 0.22,#
2 : 0.0022
2: 0.022
}
# Dict of the form: { subdom_num : density }
......@@ -175,206 +155,53 @@ densities = {
1: {'wetting': 997,
'nonwetting': 1.225},
2: {'wetting': 997,
'nonwetting': 1.225},
}
intrinsic_permeability = {
1: {"wetting": 1,
"nonwetting": 1},
2: {"wetting": 1,
"nonwetting": 1},
'nonwetting': 1.225}
}
gravity_acceleration = 9.81
L = {#
# subdom_num : subdomain L for L-scheme
1 : {'wetting' :Lw,
'nonwetting': Lnw},#
2 : {'wetting' :Lw,
'nonwetting': Lnw}
1 : {'wetting' :Lw1,
'nonwetting': Lnw1},#
2 : {'wetting' :Lw2,
'nonwetting': Lnw2}
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :lambda_w,
0 : {'wetting' :lambda_w,
'nonwetting': lambda_nw},#
2 : {'wetting' :lambda_w,
'nonwetting': lambda_nw}
}
## relative permeabilty functions on subdomain 1
def rel_perm1w(s):
# relative permeabilty wetting on subdomain1
return intrinsic_permeability[1]["wetting"]*s**2
def rel_perm1nw(s):
# relative permeabilty nonwetting on subdomain1
return intrinsic_permeability[1]["nonwetting"]*(1-s)**2
_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
subdomain1_rel_perm = {
'wetting': _rel_perm1w,#
'nonwetting': _rel_perm1nw
}
## relative permeabilty functions on subdomain 2
def rel_perm2w(s):
# relative permeabilty wetting on subdomain2
return intrinsic_permeability[2]["wetting"]*s**3
def rel_perm2nw(s):
# relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
return intrinsic_permeability[2]["nonwetting"]*(1-s)**3
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)
subdomain2_rel_perm = {
'wetting': _rel_perm2w,#
'nonwetting': _rel_perm2nw
}
## dictionary of relative permeabilties on all domains.
relative_permeability = {#
1: subdomain1_rel_perm,
2: subdomain2_rel_perm
}
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
def rel_perm1w_prime(s):
# relative permeabilty on subdomain1
return intrinsic_permeability[1]["wetting"]*2*s
def rel_perm1nw_prime(s):
# relative permeabilty on subdomain1
return -intrinsic_permeability[1]["nonwetting"]*2*(1-s)
# definition of the derivatives of the relative permeabilities
# relative permeabilty functions on subdomain 1
def rel_perm2w_prime(s):
# relative permeabilty on subdomain1
return intrinsic_permeability[2]["wetting"]*3*s**2
def rel_perm2nw_prime(s):
# relative permeabilty on subdomain1
return -intrinsic_permeability[2]["nonwetting"]*3*(1-s)**2
_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
subdomain1_rel_perm_prime = {
'wetting': _rel_perm1w_prime,
'nonwetting': _rel_perm1nw_prime
}
subdomain2_rel_perm_prime = {
'wetting': _rel_perm2w_prime,
'nonwetting': _rel_perm2nw_prime
intrinsic_permeability = {
1: 0.01,
2: 0.001,
}
# dictionary of relative permeabilties on all domains.
ka_prime = {
1: subdomain1_rel_perm_prime,
2: subdomain2_rel_perm_prime,
# RELATIVE PEMRMEABILITIES
rel_perm_definition = {
1: {"wetting": "Spow2",
"nonwetting": "oneMinusSpow2"},
2: {"wetting": "Spow3",
"nonwetting": "oneMinusSpow3"},
}
rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition)
relative_permeability = rel_perm_dict["ka"]
ka_prime = rel_perm_dict["ka_prime"]
# def saturation1(pc, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1)
#
# def saturation2(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
#
# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
# def saturation1_sym(pc, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return 1/((1 + pc)**(1/(subdomain_index + 1)))
#
#
# def saturation2_sym(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# #df.conditional(pc > 0,
# return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
#
#
# # derivative of S-pc relationship with respect to pc. This is needed for the
# # construction of a analytic solution.
# def saturation1_sym_prime(pc, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1))
#
#
# def saturation2_sym_prime(pc, n_index, alpha):
# # inverse capillary pressure-saturation-relationship
# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
#
# # note that the conditional definition of S-pc in the nonsymbolic part will be
# # incorporated in the construction of the exact solution below.
# S_pc_sym = {
# 1: ft.partial(saturation1_sym, subdomain_index = 1),
# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001),
# }
#
# S_pc_sym_prime = {
# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1),
# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001),
# }
#
# sat_pressure_relationship = {
# 1: ft.partial(saturation1, subdomain_index = 1),#,
# 2: ft.partial(saturation2, n_index=3, alpha=0.001),
# }
def saturation(pc, index):
# inverse capillary pressure-saturation-relationship
return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1)
def saturation_sym(pc, index):
# inverse capillary pressure-saturation-relationship
return 1/((1 + pc)**(1/(index + 1)))
# derivative of S-pc relationship with respect to pc. This is needed for the
# construction of a analytic solution.
def saturation_sym_prime(pc, index):
# inverse capillary pressure-saturation-relationship
return -1/((index+1)*(1 + pc)**((index+2)/(index+1)))
# note that the conditional definition of S-pc in the nonsymbolic part will be
# incorporated in the construction of the exact solution below.
S_pc_sym = {
1: ft.partial(saturation_sym, index=1),
2: ft.partial(saturation_sym, index=2),
# 3: ft.partial(saturation_sym, index=2),
# 4: ft.partial(saturation_sym, index=1)
}
S_pc_sym_prime = {
1: ft.partial(saturation_sym_prime, index=1),
2: ft.partial(saturation_sym_prime, index=2),
# 3: ft.partial(saturation_sym_prime, index=2),
# 4: ft.partial(saturation_sym_prime, index=1)
}
sat_pressure_relationship = {
1: ft.partial(saturation, index=1),
2: ft.partial(saturation, index=2),
# 3: ft.partial(saturation, index=2),
# 4: ft.partial(saturation, index=1)
# S-pc relation
Spc_on_subdomains = {
1: {"testSpc": {"index": 1}},
2: {"testSpc": {"index": 2}},
}
Spc = fts.generate_Spc_dicts(Spc_on_subdomains)
S_pc_sym = Spc["symbolic"]
S_pc_sym_prime = Spc["prime_symbolic"]
sat_pressure_relationship = Spc["dolfin"]
#############################################
# Manufacture source expressions with sympy #
......@@ -384,20 +211,15 @@ t = sym.symbols('t', positive=True)
p_e_sym = {
1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)),
'nonwetting': -1 -t*(1.1 + y*y)}, #*(1-x)**2*(1+x)**2*(1-y)**2},
2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2,
'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
pc_e_sym = 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()})
'nonwetting': (-1 -t*(1.1+ y*y))},
2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)),
'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)},
}
pc_e_sym = hlp.generate_exact_symbolic_pc(
isRichards=isRichards,
symbolic_pressure=p_e_sym
)
symbols = {"x": x,
"y": y,
......@@ -413,6 +235,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,
......@@ -423,111 +246,96 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']
# 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():
# if subdomain has no outer boundary, outer_boundary_def_points[subdomain] is None
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]}
# BOUNDARY CONDITIONS #########################################################
# 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
)
# def saturation(pressure, subdomain_index):
# # inverse capillary pressure-saturation-relationship
# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
#
# sa
f = open('TP-TP-2-patch-mesh-study.py', 'r')
# 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
# log file with ./TP-R-layered_soil.py | tee simulation.log
f = open(thisfile, 'r')
print(f.read())
f.close()
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
# 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}
)
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,
write2file=write_to_file,
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
hlp.info(simulation_parameter["use_case"])
processQueue = mp.Queue()
LDDsim = mp.Process(
target=hlp.run_simulation,
args=(
simulation_parameter,
processQueue,
starttime,
mesh_resolution
)
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, different_errornorms in subdomain_output['errornorm'].items():
filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
# for errortype, errornorm in different_errornorms.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 error_type, errornorms in different_errornorms.items():
data_dict.update(
{error_type: errornorms}
)
errors = pd.DataFrame(data_dict, index=[mesh_resolution])
# check if file exists
if os.path.isfile(filename) == 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)
LDDsim.start()
# LDDsim.join()
# hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
# )
LDDsim.join()
if mesh_study:
simulation_output_dir = processQueue.get()
hlp.merge_spacetime_errornorms(isRichards=isRichards,
resolutions=resolutions,
output_dir=simulation_output_dir)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment