diff --git a/LDDsimulation/domainSubstructuring.py b/LDDsimulation/domainSubstructuring.py index 9c127c9085f5715a6c1831b1f222621b1344df70..43a78d5c9418563234d077d9c4ddc1b0c312a227 100644 --- a/LDDsimulation/domainSubstructuring.py +++ b/LDDsimulation/domainSubstructuring.py @@ -529,3 +529,166 @@ class layeredSoilInnerPatch(domainSubstructuring): 5: self.__subdomain5_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, + } diff --git a/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py b/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py index 15f9e9dc9ed1fc26445d2029178bd42eedcb41a1..3f846e6361654c0b8f6af760d23a946b15655c4e 100755 --- a/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py +++ b/Two-phase-Richards/multi-patch/five_patch_domain_with_inner_patch/TP-R-multi-patch-with-inner-patch.py @@ -3,7 +3,6 @@ This program sets up an LDD simulation """ - import dolfin as df import sympy as sym import functools as ft @@ -12,8 +11,15 @@ 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() -# check if output directory exists +# 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 ") @@ -24,37 +30,38 @@ else: date = datetime.datetime.now() datestr = date.strftime("%Y-%m-%d") - -# init sympy session -sym.init_printing() -# solver_tol = 6E-7 +# Name of the usecase that will be printed during simulation. 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" + +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep max_iter_num = 700 FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### mesh_study = False resolutions = { - # 1: 2e-6, # h=2 - # 2: 2e-6, # h=1.1180 - # 4: 2e-6, # h=0.5590 - # 8: 2e-6, # h=0.2814 - # 16: 8e-6, # h=0.1412 - 32: 5e-6, + # 1: 1e-6, + # 2: 1e-6, + # 4: 1e-6, + 8: 1e-5, + # 16: 5e-6, + # 32: 5e-6, # 64: 2e-6, - # 128: 2e-6 + # 128: 1e-6, + # 256: 1e-6, } -# GRID ####################### -# mesh_resolution = 20 -timestep_size = 0.001 -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 gives a list of starttimes to run the simulation from. +# The list is looped over and a simulation is run with t_0 as initial time +# for each element t_0 in starttimes. starttimes = [0.0] +timestep_size = 0.001 +number_of_timesteps = 5 +# LDD scheme parameters ###################################################### Lw1 = 0.5 # /timestep_size Lnw1 = Lw1 @@ -94,23 +101,41 @@ lambda15_nw = 4 include_gravity = True -debugflag = False +debugflag = True analyse_condition = False -output_string = "./output/{}-{}_timesteps{}_P{}".format( - datestr, use_case, number_of_timesteps, FEM_Lagrange_degree - ) - - -# toggle what should be written to files +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 1 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 1 + +# fine grained control over data to be written to disk in the mesh study case +# as well as for a regular simuation for a fixed grid. if mesh_study: write_to_file = { + # output the relative errornorm (integration in space) w.r.t. an exact + # solution for each timestep into a csv file. 'space_errornorms': True, + # save the mesh and marker functions to disk 'meshes_and_markers': True, + # save xdmf/h5 data for each LDD iteration for timesteps determined by + # number_of_timesteps_to_analyse. I/O intensive! 'L_iterations_per_timestep': False, + # save solution to xdmf/h5. 'solutions': True, + # save absolute differences w.r.t an exact solution to xdmf/h5 file + # to monitor where on the domains errors happen 'absolute_differences': True, + # analyise condition numbers for timesteps determined by + # number_of_timesteps_to_analyse and save them over time to csv. 'condition_numbers': analyse_condition, + # output subsequent iteration errors measured in L^2 to csv for + # timesteps determined by number_of_timesteps_to_analyse. + # Usefull to monitor convergence of the acutal LDD solver. 'subsequent_errors': True } else: @@ -124,134 +149,115 @@ 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 -----------------------------------# -# ----------------------------------------------------------------------------# -# 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)] - -# interfaces - -interface23_vertices = [df.Point(0.0, -0.6), - df.Point(0.7, 0.0)] - -interface12_vertices = [interface23_vertices[1], - df.Point(1.0, 0.0)] - -interface13_vertices = [df.Point(0.0, 0.0), - interface23_vertices[1]] - -interface15_vertices = [df.Point(0.0, 0.0), - df.Point(0.0, 1.0)] - -interface34_vertices = [df.Point(0.0, 0.0), - interface23_vertices[0]] - -interface24_vertices = [interface23_vertices[0], - df.Point(0.0, -1.0)] - -interface45_vertices = [df.Point(-1.0, 0.0), - df.Point(0.0, 0.0)] -# subdomain1. -sub_domain1_vertices = [interface23_vertices[0], - interface23_vertices[1], - interface12_vertices[1], - sub_domain0_vertices[2], - df.Point(0.0, 1.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 -# defining the Dirichlet boundary conditions. If a domain is completely inter- -# nal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[1], - sub_domain0_vertices[2], - df.Point(0.0, 1.0)] -} -# subdomain2 -sub_domain2_vertices = [interface24_vertices[1], - sub_domain0_vertices[1], - interface12_vertices[1], - interface23_vertices[1], - interface23_vertices[0]] - -subdomain2_outer_boundary_verts = { - 0: [interface24_vertices[1], - sub_domain0_vertices[1], - interface12_vertices[1]] -} - -sub_domain3_vertices = [interface23_vertices[0], - interface23_vertices[1], - interface13_vertices[0]] - -subdomain3_outer_boundary_verts = None - - -sub_domain4_vertices = [sub_domain0_vertices[0], - interface24_vertices[1], - interface34_vertices[1], - interface34_vertices[0], - interface45_vertices[0]] - -subdomain4_outer_boundary_verts = { - 0: [interface45_vertices[0], - sub_domain0_vertices[0], - interface24_vertices[1]] -} - -sub_domain5_vertices = [interface45_vertices[0], - interface15_vertices[0], - interface15_vertices[1], - sub_domain0_vertices[3]] - -subdomain5_outer_boundary_verts = { - 0: [interface15_vertices[1], - sub_domain0_vertices[3], - interface45_vertices[0]] -} - -# 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, - sub_domain3_vertices, - sub_domain4_vertices, - sub_domain5_vertices] -# in the below list, index 0 corresponds to the 12 interface which has global -# 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 -} +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.chessBoardInnerPatch() +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 ######################################################### +# # subdomain1. +# sub_domain1_vertices = [interface23_vertices[0], +# interface23_vertices[1], +# interface12_vertices[1], +# sub_domain0_vertices[2], +# df.Point(0.0, 1.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 +# # defining the Dirichlet boundary conditions. If a domain is completely inter- +# # nal, the dictionary entry should be 0: None +# subdomain1_outer_boundary_verts = { +# 0: [interface12_vertices[1], +# sub_domain0_vertices[2], +# df.Point(0.0, 1.0)] +# } +# # subdomain2 +# sub_domain2_vertices = [interface24_vertices[1], +# sub_domain0_vertices[1], +# interface12_vertices[1], +# interface23_vertices[1], +# interface23_vertices[0]] +# +# subdomain2_outer_boundary_verts = { +# 0: [interface24_vertices[1], +# sub_domain0_vertices[1], +# interface12_vertices[1]] +# } +# +# sub_domain3_vertices = [interface23_vertices[0], +# interface23_vertices[1], +# interface13_vertices[0]] +# +# subdomain3_outer_boundary_verts = None +# +# +# sub_domain4_vertices = [sub_domain0_vertices[0], +# interface24_vertices[1], +# interface34_vertices[1], +# interface34_vertices[0], +# interface45_vertices[0]] +# +# subdomain4_outer_boundary_verts = { +# 0: [interface45_vertices[0], +# sub_domain0_vertices[0], +# interface24_vertices[1]] +# } +# +# sub_domain5_vertices = [interface45_vertices[0], +# interface15_vertices[0], +# interface15_vertices[1], +# sub_domain0_vertices[3]] +# +# subdomain5_outer_boundary_verts = { +# 0: [interface15_vertices[1], +# sub_domain0_vertices[3], +# interface45_vertices[0]] +# } +# +# # 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, +# sub_domain3_vertices, +# sub_domain4_vertices, +# sub_domain5_vertices] +# # in the below list, index 0 corresponds to the 12 interface which has global +# # 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 = { 1: True, @@ -650,17 +656,11 @@ p_e_sym = { 'nonwetting': 0*t }, } -pc_e_sym = dict() -for subdomain, isR in isRichards.items(): - if isR: - 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']} +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym ) - symbols = {"x": x, "y": y, "t": t} @@ -685,35 +685,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. - -# 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 + ) +# 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 @@ -723,88 +706,63 @@ print(f.read()) f.close() -# RUN ######################################################################### -for starttime in starttimes: - for mesh_resolution, solver_tol in resolutions.items(): - # initialise LDD simulation class - simulation = ldd.LDDsimulation( - tol=1E-14, - LDDsolver_tol=solver_tol, - debug=debugflag, - max_iter_num=max_iter_num, - FEM_Lagrange_degree=FEM_Lagrange_degree, - mesh_study=mesh_study - ) - - simulation.set_parameters( - use_case=use_case, - output_dir=output_string, - subdomain_def_points=subdomain_def_points, - isRichards=isRichards, - interface_def_points=interface_def_points, - outer_boundary_def_points=outer_boundary_def_points, - adjacent_subdomains=adjacent_subdomains, - mesh_resolution=mesh_resolution, - viscosity=viscosity, - porosity=porosity, - L=L, - lambda_param=lambda_param, - relative_permeability=relative_permeability, - saturation=sat_pressure_relationship, - starttime=starttime, - number_of_timesteps=number_of_timesteps, - number_of_timesteps_to_analyse=number_of_timesteps_to_analyse, - plot_timestep_every=plot_timestep_every, - timestep_size=timestep_size, - sources=source_expression, - initial_conditions=initial_condition, - dirichletBC_expression_strings=dirichletBC, - exact_solution=exact_solution, - densities=densities, - include_gravity=include_gravity, - gravity_acceleration=gravity_acceleration, - write2file=write_to_file, - ) - - simulation.initialise() - output_dir = simulation.output_dir - # simulation.write_exact_solution_to_xdmf() - output = simulation.run(analyse_condition=analyse_condition) - for subdomain_index, subdomain_output in output.items(): - mesh_h = subdomain_output['mesh_size'] - for phase, error_dict in subdomain_output['errornorm'].items(): - filename = output_dir \ - + "subdomain{}".format(subdomain_index)\ - + "-space-time-errornorm-{}-phase.csv".format(phase) - # for errortype, errornorm in error_dict.items(): - - # eocfile = open("eoc_filename", "a") - # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) - # eocfile.close() - # if subdomain.isRichards:mesh_h - data_dict = { - 'mesh_parameter': mesh_resolution, - 'mesh_h': mesh_h, - } - for norm_type, errornorm in error_dict.items(): - data_dict.update( - {norm_type: errornorm} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) is True: - with open(filename, 'a') as f: - errors.to_csv( - f, - header=False, - sep='\t', - encoding='utf-8', - index=False +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": output_string, + "subdomain_def_points": subdomain_def_points, + "isRichards": isRichards, + "interface_def_points": interface_def_points, + "outer_boundary_def_points": outer_boundary_def_points, + "adjacent_subdomains": adjacent_subdomains, + # "mesh_resolution": mesh_resolution, + "viscosity": viscosity, + "porosity": porosity, + "L": L, + "lambda_param": lambda_param, + "relative_permeability": relative_permeability, + "sat_pressure_relationship": sat_pressure_relationship, + # "starttime": starttime, + "number_of_timesteps": number_of_timesteps, + "number_of_timesteps_to_analyse": number_of_timesteps_to_analyse, + "plot_timestep_every": plot_timestep_every, + "timestep_size": timestep_size, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for starttime in starttimes: + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + starttime, + mesh_resolution ) - else: - errors.to_csv( - filename, - sep='\t', - encoding='utf-8', - index=False ) + LDDsim.start() + LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # )