Skip to content
Snippets Groups Projects
Commit 1f45e155 authored by David's avatar David
Browse files

clean up TPR two-patch mesh studies at the aid of helpers and substructuring

parent aab21c2e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 495 additions and 1157 deletions
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one-but-g.py
+ 83
− 194
View file @ 1f45e155
......@@ -3,7 +3,6 @@
This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
import functools as ft
......@@ -12,6 +11,8 @@ 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()
......@@ -120,80 +121,16 @@ else:
}
# OUTPUT FILE STRING #########################################################
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
)
# DOMAIN AND INTERFACE #######################################################
# global 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]]
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
# MODEL CONFIGURATION #########################################################
isRichards = {
......@@ -430,15 +367,10 @@ p_e_sym = {
'nonwetting': (-1.0-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
pc_e_sym = 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,
......@@ -464,36 +396,18 @@ 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.
# BOUNDARY CONDITIONS #########################################################
# 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
......@@ -503,88 +417,63 @@ print(f.read())
f.close()
# RUN #########################################################################
# 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,
"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():
# 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_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
)
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
)
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
# )
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study-all-params-one.py
+ 83
− 195
View file @ 1f45e155
......@@ -3,7 +3,6 @@
This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
import functools as ft
......@@ -12,6 +11,8 @@ 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()
......@@ -120,80 +121,16 @@ else:
}
# OUTPUT FILE STRING #########################################################
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
)
# DOMAIN AND INTERFACE #######################################################
# global 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]]
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
# MODEL CONFIGURATION #########################################################
isRichards = {
......@@ -430,15 +367,10 @@ p_e_sym = {
'nonwetting': (-1.0-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
pc_e_sym = 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,
......@@ -464,36 +396,18 @@ 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.
# BOUNDARY CONDITIONS #########################################################
# 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
......@@ -503,89 +417,63 @@ print(f.read())
f.close()
# RUN #########################################################################
# 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,
"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():
# 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_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
)
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
)
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
# )
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-mesh-study.py
+ 83
− 194
View file @ 1f45e155
......@@ -3,7 +3,6 @@
This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
import functools as ft
......@@ -12,6 +11,8 @@ 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()
......@@ -120,80 +121,16 @@ else:
}
# OUTPUT FILE STRING #########################################################
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
)
# DOMAIN AND INTERFACE #######################################################
# global 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]]
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
# MODEL CONFIGURATION #########################################################
isRichards = {
......@@ -430,15 +367,10 @@ p_e_sym = {
'nonwetting': (-1-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
pc_e_sym = 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,
......@@ -464,36 +396,18 @@ 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.
# BOUNDARY CONDITIONS #########################################################
# 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
......@@ -503,88 +417,63 @@ print(f.read())
f.close()
# RUN #########################################################################
# 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,
"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():
# 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_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
)
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
)
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
# )
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm-mesh-study.py
+ 83
− 195
View file @ 1f45e155
......@@ -3,7 +3,6 @@
This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
import functools as ft
......@@ -12,6 +11,8 @@ 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()
......@@ -120,82 +121,16 @@ else:
}
# OUTPUT FILE STRING #########################################################
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
)
# DOMAIN AND INTERFACE #######################################################
# global 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]]
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
# MODEL CONFIGURATION #########################################################
isRichards = {
......@@ -444,15 +379,10 @@ p_e_sym = {
'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
pc_e_sym = 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,
......@@ -478,35 +408,18 @@ 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.
# BOUNDARY CONDITIONS #########################################################
# 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
......@@ -516,88 +429,63 @@ print(f.read())
f.close()
# RUN #########################################################################
# 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,
"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():
# 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_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
)
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
)
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
# )
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm-mesh-study.py
+ 83
− 195
View file @ 1f45e155
......@@ -3,7 +3,6 @@
This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
import functools as ft
......@@ -12,6 +11,8 @@ 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()
......@@ -120,82 +121,16 @@ else:
}
# OUTPUT FILE STRING #########################################################
if mesh_study:
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
else:
for tol in resolutions.values():
solver_tol = tol
output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol
)
# DOMAIN AND INTERFACE #######################################################
# global 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]]
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
# MODEL CONFIGURATION #########################################################
isRichards = {
......@@ -444,15 +379,10 @@ p_e_sym = {
'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
pc_e_sym = 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,
......@@ -478,35 +408,18 @@ 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.
# BOUNDARY CONDITIONS #########################################################
# 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
......@@ -516,88 +429,63 @@ print(f.read())
f.close()
# RUN #########################################################################
# 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,
"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():
# 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_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
)
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
)
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
# )
Two-phase-Richards/two-patch/TP-R-two-patch-test-case/mesh_studies/TP-R-2-patch-realistic-same-intrinsic-perm.py
+ 80
− 184
View file @ 1f45e155
......@@ -3,7 +3,6 @@
This program sets up an LDD simulation
"""
import dolfin as df
import sympy as sym
import functools as ft
......@@ -12,6 +11,8 @@ 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()
......@@ -125,69 +126,11 @@ output_string = "./output/{}-{}_timesteps{}_P{}".format(
)
# 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]]
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
# MODEL CONFIGURATION #########################################################
isRichards = {
......@@ -436,15 +379,10 @@ p_e_sym = {
'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2},
} #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
pc_e_sym = 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,
......@@ -470,35 +408,18 @@ 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.
# BOUNDARY CONDITIONS #########################################################
# 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
......@@ -508,88 +429,63 @@ print(f.read())
f.close()
# RUN #########################################################################
# 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,
"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():
# 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}
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
)
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
)
LDDsim.start()
LDDsim.join()
# hlp.run_simulation(
# mesh_resolution=mesh_resolution,
# starttime=starttime,
# parameter=simulation_parameter
# )
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment