Skip to content
Snippets Groups Projects
Commit 16ea3a3e authored by David's avatar David
Browse files

clean up multi-patch with inner patch by helpers and substructuring

parent 1f45e155
No related branches found
No related tags found
No related merge requests found
...@@ -529,3 +529,166 @@ class layeredSoilInnerPatch(domainSubstructuring): ...@@ -529,3 +529,166 @@ class layeredSoilInnerPatch(domainSubstructuring):
5: self.__subdomain5_outer_boundary_verts, 5: self.__subdomain5_outer_boundary_verts,
6: self.__subdomain6_outer_boundary_verts 6: self.__subdomain6_outer_boundary_verts
} }
class chessBoardInnerPatch(domainSubstructuring):
"""layered soil substructuring with inner patch."""
def __init__(self):
"""Layered soil case with inner patch."""
super().__init__()
hlp.print_once("\n Layered Soil with inner Patch:\n")
# global domain
self.__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)
]
self.__interface_def_points()
self.__adjacent_subdomains()
self.__subdomain_def_points()
self.__outer_boundary_def_points()
def __interface_def_points(self):
"""Set self._interface_def_points."""
self.__interface23_vertices = [
df.Point(0.0, -0.6),
df.Point(0.7, 0.0)]
self.__interface12_vertices = [
self.__interface23_vertices[1],
df.Point(1.0, 0.0)]
self.__interface13_vertices = [
df.Point(0.0, 0.0),
self.__interface23_vertices[1]]
self.__interface15_vertices = [
df.Point(0.0, 0.0),
df.Point(0.0, 1.0)]
self.__interface34_vertices = [
df.Point(0.0, 0.0),
self.__interface23_vertices[0]]
self.__interface24_vertices = [
self.__interface23_vertices[0],
df.Point(0.0, -1.0)]
self.__interface45_vertices = [
df.Point(-1.0, 0.0),
df.Point(0.0, 0.0)]
# interface_vertices introduces a global numbering of interfaces.
self.interface_def_points = [
self.__interface13_vertices,
self.__interface12_vertices,
self.__interface23_vertices,
self.__interface24_vertices,
self.__interface34_vertices,
self.__interface45_vertices,
self.__interface15_vertices,
]
def __adjacent_subdomains(self):
"""Set self._adjacent_subdomains."""
self.adjacent_subdomains = [
[1, 3],
[1, 2],
[2, 3],
[2, 4],
[3, 4],
[4, 5],
[1, 5]
]
def __subdomain_def_points(self):
"""Set self._subdomain_def_points."""
# subdomain1.
self.__subdomain1_vertices = [
self.__interface23_vertices[0],
self.__interface23_vertices[1],
self.__interface12_vertices[1],
self.__subdomain0_vertices[2],
df.Point(0.0, 1.0)]
# subdomain2
self.__subdomain2_vertices = [
self.__interface24_vertices[1],
self.__subdomain0_vertices[1],
self.__interface12_vertices[1],
self.__interface23_vertices[1],
self.__interface23_vertices[0]]
self.__subdomain3_vertices = [
self.__interface23_vertices[0],
self.__interface23_vertices[1],
self.__interface13_vertices[0]]
self.__subdomain4_vertices = [
self.__subdomain0_vertices[0],
self.__interface24_vertices[1],
self.__interface34_vertices[1],
self.__interface34_vertices[0],
self.__interface45_vertices[0]]
self.__subdomain5_vertices = [
self.__interface45_vertices[0],
self.__interface15_vertices[0],
self.__interface15_vertices[1],
self.__subdomain0_vertices[3]]
self.subdomain_def_points = [
self.__subdomain0_vertices,
self.__subdomain1_vertices,
self.__subdomain2_vertices,
self.__subdomain3_vertices,
self.__subdomain4_vertices,
self.__subdomain5_vertices,
]
def __outer_boundary_def_points(self):
"""Set self._outer_boundary_def_points."""
# 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
self.__subdomain1_outer_boundary_verts = {
0: [self.__interface12_vertices[1],
self.__subdomain0_vertices[2],
df.Point(0.0, 1.0)]
}
self.__subdomain2_outer_boundary_verts = {
0: [self.__interface24_vertices[1],
self.__subdomain0_vertices[1],
self.__interface12_vertices[1]]
}
self.__subdomain3_outer_boundary_verts = None
self.__subdomain4_outer_boundary_verts = {
0: [self.__interface45_vertices[0],
self.__subdomain0_vertices[0],
self.__interface24_vertices[1]]
}
self.__subdomain5_outer_boundary_verts = {
0: [self.__interface15_vertices[1],
self.__subdomain0_vertices[3],
self.__interface45_vertices[0]]
}
# if a subdomain has no outer boundary write None instead, i.e.
# i: None
# if i is the index of the inner subdomain.
self.outer_boundary_def_points = {
1: self.__subdomain1_outer_boundary_verts,
2: self.__subdomain2_outer_boundary_verts,
3: self.__subdomain3_outer_boundary_verts,
4: self.__subdomain4_outer_boundary_verts,
5: self.__subdomain5_outer_boundary_verts,
}
...@@ -3,7 +3,6 @@ ...@@ -3,7 +3,6 @@
This program sets up an LDD simulation This program sets up an LDD simulation
""" """
import dolfin as df import dolfin as df
import sympy as sym import sympy as sym
import functools as ft import functools as ft
...@@ -12,8 +11,15 @@ import helpers as hlp ...@@ -12,8 +11,15 @@ import helpers as hlp
import datetime import datetime
import os import os
import pandas as pd import pandas as pd
import multiprocessing as mp
import domainSubstructuring as dss
# check if output directory exists # init sympy session
sym.init_printing()
# PREREQUISITS ###############################################################
# check if output directory "./output" exists. This will be used in
# the generation of the output string.
if not os.path.exists('./output'): if not os.path.exists('./output'):
os.mkdir('./output') os.mkdir('./output')
print("Directory ", './output', " created ") print("Directory ", './output', " created ")
...@@ -24,37 +30,38 @@ else: ...@@ -24,37 +30,38 @@ else:
date = datetime.datetime.now() date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d") datestr = date.strftime("%Y-%m-%d")
# Name of the usecase that will be printed during simulation.
# init sympy session
sym.init_printing()
# solver_tol = 6E-7
use_case = "TP-R-five-domain-with-inner-patch-realistic" use_case = "TP-R-five-domain-with-inner-patch-realistic"
# name of this very file. Needed for log output. # The name of this very file. Needed for creating log output.
thisfile = "TP-R-multi-patch-with-inner-patch.py" thisfile = "TP-R-multi-patch-with-inner-patch.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
max_iter_num = 700 max_iter_num = 700
FEM_Lagrange_degree = 1 FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
mesh_study = False mesh_study = False
resolutions = { resolutions = {
# 1: 2e-6, # h=2 # 1: 1e-6,
# 2: 2e-6, # h=1.1180 # 2: 1e-6,
# 4: 2e-6, # h=0.5590 # 4: 1e-6,
# 8: 2e-6, # h=0.2814 8: 1e-5,
# 16: 8e-6, # h=0.1412 # 16: 5e-6,
32: 5e-6, # 32: 5e-6,
# 64: 2e-6, # 64: 2e-6,
# 128: 2e-6 # 128: 1e-6,
# 256: 1e-6,
} }
# GRID ####################### # starttimes gives a list of starttimes to run the simulation from.
# mesh_resolution = 20 # The list is looped over and a simulation is run with t_0 as initial time
timestep_size = 0.001 # for each element t_0 in starttimes.
number_of_timesteps = 1000
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 = 8
starttimes = [0.0] starttimes = [0.0]
timestep_size = 0.001
number_of_timesteps = 5
# LDD scheme parameters ######################################################
Lw1 = 0.5 # /timestep_size Lw1 = 0.5 # /timestep_size
Lnw1 = Lw1 Lnw1 = Lw1
...@@ -94,23 +101,41 @@ lambda15_nw = 4 ...@@ -94,23 +101,41 @@ lambda15_nw = 4
include_gravity = True include_gravity = True
debugflag = False debugflag = True
analyse_condition = False analyse_condition = False
output_string = "./output/{}-{}_timesteps{}_P{}".format( # I/O CONFIG #################################################################
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree # when number_of_timesteps is high, it might take a long time to write all
) # timesteps to disk. Therefore, you can choose to only write data of every
# plot_timestep_every timestep to disk.
plot_timestep_every = 1
# toggle what should be written to files # Decide how many timesteps you want analysed. Analysed means, that
# subsequent errors of the L-iteration within the timestep are written out.
number_of_timesteps_to_analyse = 1
# fine grained control over data to be written to disk in the mesh study case
# as well as for a regular simuation for a fixed grid.
if mesh_study: if mesh_study:
write_to_file = { 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, 'space_errornorms': True,
# save the mesh and marker functions to disk
'meshes_and_markers': True, 'meshes_and_markers': True,
# save xdmf/h5 data for each LDD iteration for timesteps determined by
# number_of_timesteps_to_analyse. I/O intensive!
'L_iterations_per_timestep': False, 'L_iterations_per_timestep': False,
# save solution to xdmf/h5.
'solutions': True, '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, '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, '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 'subsequent_errors': True
} }
else: else:
...@@ -124,134 +149,115 @@ else: ...@@ -124,134 +149,115 @@ else:
'subsequent_errors': True 'subsequent_errors': True
} }
# OUTPUT FILE STRING #########################################################
output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
# ----------------------------------------------------------------------------# # DOMAIN AND INTERFACE #######################################################
# ------------------- Domain and Interface -----------------------------------# substructuring = dss.chessBoardInnerPatch()
# ----------------------------------------------------------------------------# interface_def_points = substructuring.interface_def_points
# global simulation domain domain adjacent_subdomains = substructuring.adjacent_subdomains
sub_domain0_vertices = [df.Point(-1.0, -1.0), subdomain_def_points = substructuring.subdomain_def_points
df.Point(1.0, -1.0), outer_boundary_def_points = substructuring.outer_boundary_def_points
df.Point(1.0, 1.0),
df.Point(-1.0, 1.0)] # MODEL CONFIGURATION #########################################################
# # subdomain1.
# interfaces # sub_domain1_vertices = [interface23_vertices[0],
# interface23_vertices[1],
interface23_vertices = [df.Point(0.0, -0.6), # interface12_vertices[1],
df.Point(0.7, 0.0)] # sub_domain0_vertices[2],
# df.Point(0.0, 1.0)]
interface12_vertices = [interface23_vertices[1], #
df.Point(1.0, 0.0)] # # 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
interface13_vertices = [df.Point(0.0, 0.0), # # defining the Dirichlet boundary conditions. If a domain is completely inter-
interface23_vertices[1]] # # nal, the dictionary entry should be 0: None
# subdomain1_outer_boundary_verts = {
interface15_vertices = [df.Point(0.0, 0.0), # 0: [interface12_vertices[1],
df.Point(0.0, 1.0)] # sub_domain0_vertices[2],
# df.Point(0.0, 1.0)]
interface34_vertices = [df.Point(0.0, 0.0), # }
interface23_vertices[0]] # # subdomain2
# sub_domain2_vertices = [interface24_vertices[1],
interface24_vertices = [interface23_vertices[0], # sub_domain0_vertices[1],
df.Point(0.0, -1.0)] # interface12_vertices[1],
# interface23_vertices[1],
interface45_vertices = [df.Point(-1.0, 0.0), # interface23_vertices[0]]
df.Point(0.0, 0.0)] #
# subdomain1. # subdomain2_outer_boundary_verts = {
sub_domain1_vertices = [interface23_vertices[0], # 0: [interface24_vertices[1],
interface23_vertices[1], # sub_domain0_vertices[1],
interface12_vertices[1], # interface12_vertices[1]]
sub_domain0_vertices[2], # }
df.Point(0.0, 1.0)] #
# sub_domain3_vertices = [interface23_vertices[0],
# vertex coordinates of the outer boundaries. If it can not be specified as a # interface23_vertices[1],
# polygon, use an entry per boundary polygon. This information is used for # interface13_vertices[0]]
# defining the Dirichlet boundary conditions. If a domain is completely inter- #
# nal, the dictionary entry should be 0: None # subdomain3_outer_boundary_verts = None
subdomain1_outer_boundary_verts = { #
0: [interface12_vertices[1], #
sub_domain0_vertices[2], # sub_domain4_vertices = [sub_domain0_vertices[0],
df.Point(0.0, 1.0)] # interface24_vertices[1],
} # interface34_vertices[1],
# subdomain2 # interface34_vertices[0],
sub_domain2_vertices = [interface24_vertices[1], # interface45_vertices[0]]
sub_domain0_vertices[1], #
interface12_vertices[1], # subdomain4_outer_boundary_verts = {
interface23_vertices[1], # 0: [interface45_vertices[0],
interface23_vertices[0]] # sub_domain0_vertices[0],
# interface24_vertices[1]]
subdomain2_outer_boundary_verts = { # }
0: [interface24_vertices[1], #
sub_domain0_vertices[1], # sub_domain5_vertices = [interface45_vertices[0],
interface12_vertices[1]] # interface15_vertices[0],
} # interface15_vertices[1],
# sub_domain0_vertices[3]]
sub_domain3_vertices = [interface23_vertices[0], #
interface23_vertices[1], # subdomain5_outer_boundary_verts = {
interface13_vertices[0]] # 0: [interface15_vertices[1],
# sub_domain0_vertices[3],
subdomain3_outer_boundary_verts = None # interface45_vertices[0]]
# }
#
sub_domain4_vertices = [sub_domain0_vertices[0], # # list of subdomains given by the boundary polygon vertices.
interface24_vertices[1], # # Subdomains are given as a list of dolfin points forming
interface34_vertices[1], # # a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
interface34_vertices[0], # # to create the subdomain. subdomain_def_points[0] contains the
interface45_vertices[0]] # # vertices of the global simulation domain and subdomain_def_points[i] contains
# # the vertices of the subdomain i.
subdomain4_outer_boundary_verts = { # subdomain_def_points = [sub_domain0_vertices,
0: [interface45_vertices[0], # sub_domain1_vertices,
sub_domain0_vertices[0], # sub_domain2_vertices,
interface24_vertices[1]] # sub_domain3_vertices,
} # sub_domain4_vertices,
# sub_domain5_vertices]
sub_domain5_vertices = [interface45_vertices[0], # # in the below list, index 0 corresponds to the 12 interface which has global
interface15_vertices[0], # # marker value 1
interface15_vertices[1], # interface_def_points = [interface13_vertices,
sub_domain0_vertices[3]] # interface12_vertices,
# interface23_vertices,
subdomain5_outer_boundary_verts = { # interface24_vertices,
0: [interface15_vertices[1], # interface34_vertices,
sub_domain0_vertices[3], # interface45_vertices,
interface45_vertices[0]] # interface15_vertices,]
} #
# # adjacent_subdomains[i] contains the indices of the subdomains sharing the
# list of subdomains given by the boundary polygon vertices. # # interface i (i.e. given by interface_def_points[i]).
# Subdomains are given as a list of dolfin points forming # adjacent_subdomains = [[1, 3], [1, 2], [2, 3], [2, 4], [3, 4], [4, 5], [1, 5]]
# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used #
# to create the subdomain. subdomain_def_points[0] contains the # # if a subdomain has no outer boundary write None instead, i.e.
# vertices of the global simulation domain and subdomain_def_points[i] contains # # i: None
# the vertices of the subdomain i. # # if i is the index of the inner subdomain.
subdomain_def_points = [sub_domain0_vertices, # outer_boundary_def_points = {
sub_domain1_vertices, # # subdomain number
sub_domain2_vertices, # 1: subdomain1_outer_boundary_verts,
sub_domain3_vertices, # 2: subdomain2_outer_boundary_verts,
sub_domain4_vertices, # 3: subdomain3_outer_boundary_verts,
sub_domain5_vertices] # 4: subdomain4_outer_boundary_verts,
# in the below list, index 0 corresponds to the 12 interface which has global # 5: subdomain5_outer_boundary_verts
# marker value 1 # }
interface_def_points = [interface13_vertices,
interface12_vertices,
interface23_vertices,
interface24_vertices,
interface34_vertices,
interface45_vertices,
interface15_vertices,]
# adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1, 3], [1, 2], [2, 3], [2, 4], [3, 4], [4, 5], [1, 5]]
# 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
}
isRichards = { isRichards = {
1: True, 1: True,
...@@ -650,17 +656,11 @@ p_e_sym = { ...@@ -650,17 +656,11 @@ p_e_sym = {
'nonwetting': 0*t }, 'nonwetting': 0*t },
} }
pc_e_sym = dict() pc_e_sym = hlp.generate_exact_symbolic_pc(
for subdomain, isR in isRichards.items(): isRichards=isRichards,
if isR: symbolic_pressure=p_e_sym
pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
else:
pc_e_sym.update(
{subdomain: p_e_sym[subdomain]['nonwetting']
- p_e_sym[subdomain]['wetting']}
) )
symbols = {"x": x, symbols = {"x": x,
"y": y, "y": y,
"t": t} "t": t}
...@@ -685,35 +685,18 @@ source_expression = exact_solution_example['source'] ...@@ -685,35 +685,18 @@ source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution'] exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition'] initial_condition = exact_solution_example['initial_condition']
# Dictionary of dirichlet boundary conditions. # BOUNDARY CONDITIONS #########################################################
dirichletBC = dict() # Dictionary of dirichlet boundary conditions. If an exact solution case is
# similarly to the outer boundary dictionary, if a patch has no outer boundary # used, use the hlp.generate_exact_DirichletBC() method to generate the
# None should be written instead of an expression. # Dirichlet Boundary conditions from the exact solution.
# This is a bit of a brainfuck: dirichletBC = hlp.generate_exact_DirichletBC(
# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind. isRichards=isRichards,
# Since a domain patch can have several disjoint outer boundary parts, the outer_boundary_def_points=outer_boundary_def_points,
# expressions need to get an enumaration index which starts at 0. exact_solution=exact_solution
# 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]}
) )
# 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 ############################################################# # LOG FILE OUTPUT #############################################################
# read this file and print it to std out. This way the simulation can produce a # read this file and print it to std out. This way the simulation can produce a
...@@ -723,88 +706,63 @@ print(f.read()) ...@@ -723,88 +706,63 @@ print(f.read())
f.close() 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 starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items(): for mesh_resolution, solver_tol in resolutions.items():
# initialise LDD simulation class simulation_parameter.update({"solver_tol": solver_tol})
simulation = ldd.LDDsimulation( hlp.info(simulation_parameter["use_case"])
tol=1E-14, LDDsim = mp.Process(
LDDsolver_tol=solver_tol, target=hlp.run_simulation,
debug=debugflag, args=(
max_iter_num=max_iter_num, simulation_parameter,
FEM_Lagrange_degree=FEM_Lagrange_degree, starttime,
mesh_study=mesh_study 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
# )
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment