diff --git a/LDDsimulation/domainSubstructuring.py b/LDDsimulation/domainSubstructuring.py index bbedaf8ef9fe465f3645fae7fc09d6a50862820a..9c127c9085f5715a6c1831b1f222621b1344df70 100644 --- a/LDDsimulation/domainSubstructuring.py +++ b/LDDsimulation/domainSubstructuring.py @@ -34,6 +34,7 @@ class domainSubstructuring(object): """Set self._outer_boundary_def_points.""" raise(NotImplementedError()) + class twoSoilLayers(domainSubstructuring): """layered soil substructuring with inner patch.""" @@ -56,7 +57,6 @@ class twoSoilLayers(domainSubstructuring): def __interface_def_points(self): """Set self._interface_def_points.""" - self.__interface12_vertices = [ df.Point(-1.0, 0.0), df.Point(1.0, 0.0) @@ -127,14 +127,20 @@ class twoSoilLayers(domainSubstructuring): } - -class layeredSoilInnerPatch(domainSubstructuring): +class layeredSoil(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() @@ -143,6 +149,150 @@ class layeredSoilInnerPatch(domainSubstructuring): def __interface_def_points(self): """Set self._interface_def_points.""" + self.__interface12_vertices = [ + df.Point(-1.0, 0.8), + df.Point(0.3, 0.8), + df.Point(0.5, 0.9), + df.Point(0.8, 0.7), + df.Point(1.0, 0.65) + ] + + # interface23 + self.__interface23_vertices = [ + df.Point(-1.0, 0.0), + # df.Point(-0.35, 0.0), + # df.Point(0.0, 0.0), + # df.Point(0.5, 0.0), + # df.Point(0.85, 0.0), + df.Point(1.0, 0.0) + ] + + self.__interface34_vertices = [ + df.Point(-1.0, -0.6), + df.Point(-0.6, -0.45), + df.Point(0.3, -0.25), + df.Point(0.65, -0.6), + df.Point(1.0, -0.7) + ] + + # interface_vertices introduces a global numbering of interfaces. + self.interface_def_points = [ + self.__interface12_vertices, + self.__interface23_vertices, + self.__interface34_vertices, + ] + + def __adjacent_subdomains(self): + """Set self._adjacent_subdomains.""" + self.adjacent_subdomains = [ + [1, 2], + [2, 3], + [3, 4], + ] + + def __subdomain_def_points(self): + """Set self._subdomain_def_points.""" + self.__subdomain1_vertices = [ + self.__interface12_vertices[0], + self.__interface12_vertices[1], + self.__interface12_vertices[2], + self.__interface12_vertices[3], + self.__interface12_vertices[4], + self.__subdomain0_vertices[2], + self.__subdomain0_vertices[3]] + + self.__subdomain2_vertices = [ + self.__interface23_vertices[0], + self.__interface23_vertices[1], + self.__subdomain1_vertices[4], + self.__subdomain1_vertices[3], + self.__subdomain1_vertices[2], + self.__subdomain1_vertices[1], + self.__subdomain1_vertices[0]] + + self.__subdomain3_vertices = [ + self.__interface34_vertices[0], + self.__interface34_vertices[1], + self.__interface34_vertices[2], + self.__interface34_vertices[3], + self.__interface34_vertices[4], + self.__subdomain2_vertices[1], + self.__subdomain2_vertices[0] + ] + + # subdomain3 + self.__subdomain4_vertices = [ + self.__subdomain0_vertices[0], + self.__subdomain0_vertices[1], + self.__subdomain3_vertices[4], + self.__subdomain3_vertices[3], + self.__subdomain3_vertices[2], + self.__subdomain3_vertices[1], + self.__subdomain3_vertices[0] + ] + + self.subdomain_def_points = [ + self.__subdomain0_vertices, + self.__subdomain1_vertices, + self.__subdomain2_vertices, + self.__subdomain3_vertices, + self.__subdomain4_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[4], + self.__subdomain0_vertices[2], + self.__subdomain0_vertices[3], + self.__interface12_vertices[0]] + } + + self.__subdomain2_outer_boundary_verts = { + 0: [self.__interface23_vertices[1], + self.__subdomain1_vertices[4]], + 1: [self.__subdomain1_vertices[0], + self.__interface23_vertices[0]] + } + + self.__subdomain3_outer_boundary_verts = { + 0: [self.__interface34_vertices[4], + self.__subdomain2_vertices[1]], + 1: [self.__subdomain2_vertices[0], + self.__interface34_vertices[0]] + } + + # if a subdomain has no outer boundary write None instead, i.e. + # i: None + self.__subdomain4_outer_boundary_verts = { + 0: [self.__subdomain4_vertices[6], + self.__subdomain4_vertices[0], + self.__subdomain4_vertices[1], + self.__subdomain4_vertices[2]] + } + + # if i is the index of the inner subdomain. + self.outer_boundary_def_points = { + # subdomain number + 1: self.__subdomain1_outer_boundary_verts, + 2: self.__subdomain2_outer_boundary_verts, + 3: self.__subdomain3_outer_boundary_verts, + 4: self.__subdomain4_outer_boundary_verts, + } + + +class layeredSoilInnerPatch(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), @@ -151,6 +301,13 @@ class layeredSoilInnerPatch(domainSubstructuring): 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.__interface12_vertices = [ df.Point(-1.0, 0.8), df.Point(0.3, 0.8), diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py index 3b054bfd98e27d85545e780d47ab1d8b48f7a16b..48766dab7d33ddb4d8aad5a46f5be49cb5a7998d 100755 --- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.py +++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-all-params-one.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,7 +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() @@ -59,7 +59,7 @@ resolutions = { # for each element t_0 in starttimes. starttimes = [0.0] timestep_size = 0.001 -number_of_timesteps = 20 +number_of_timesteps = 5 # LDD scheme parameters ###################################################### @@ -72,25 +72,25 @@ Lnw3 = Lw3 Lw4 = 0.025 # /timestep_size Lnw4 = Lw4 -lambda12_w = 40 -lambda12_nw = 40 -lambda23_w = 40 -lambda23_nw = 40 -lambda34_w = 40 -lambda34_nw = 40 +lambda12_w = 4 +lambda12_nw = 4 +lambda23_w = 4 +lambda23_nw = 4 +lambda34_w = 4 +lambda34_nw = 4 include_gravity = False debugflag = False -analyse_condition = True +analyse_condition = False # I/O CONFIG ################################################################# # when number_of_timesteps is high, it might take a long time to write all # timesteps to disk. Therefore, you can choose to only write data of every # plot_timestep_every timestep to disk. -plot_timestep_every = 4 +plot_timestep_every = 1 # Decide how many timesteps you want analysed. Analysed means, that # subsequent errors of the L-iteration within the timestep are written out. -number_of_timesteps_to_analyse = 5 +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. @@ -129,168 +129,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 - ) - +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) # DOMAIN AND INTERFACE ####################################################### -# global domain -subdomain0_vertices = [ - df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0) - ] - - -interface12_vertices = [df.Point(-1.0, 0.8), - df.Point(0.3, 0.8), - df.Point(0.5, 0.9), - df.Point(0.8, 0.7), - df.Point(1.0, 0.65)] -# subdomain1. -subdomain1_vertices = [ - interface12_vertices[0], - interface12_vertices[1], - interface12_vertices[2], - interface12_vertices[3], - interface12_vertices[4], # southern boundary, 12 interface - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3] # northern boundary, outer on_boundary - ] - - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[4], - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3], - interface12_vertices[0]] -} - - -# interface23 -interface23_vertices = [ - df.Point(-1.0, 0.0), - df.Point(-0.35, 0.0), - # df.Point(6.5, 4.5), - df.Point(0.0, 0.0), - df.Point(0.5, 0.0), - # df.Point(11.5, 3.5), - # df.Point(13.0, 3) - df.Point(0.85, 0.0), - df.Point(1.0, 0.0) - ] - -# subdomain1 -subdomain2_vertices = [ - interface23_vertices[0], - interface23_vertices[1], - interface23_vertices[2], - interface23_vertices[3], - interface23_vertices[4], - interface23_vertices[5], # southern boundary, 23 interface - subdomain1_vertices[4], # eastern boundary, outer boundary - subdomain1_vertices[3], - subdomain1_vertices[2], - subdomain1_vertices[1], - subdomain1_vertices[0] # northern boundary, 12 interface - ] - -subdomain2_outer_boundary_verts = { - 0: [interface23_vertices[5], - subdomain1_vertices[4]], - 1: [subdomain1_vertices[0], - interface23_vertices[0]] -} - - -# interface34 -interface34_vertices = [df.Point(-1.0, -0.6), - df.Point(-0.6, -0.45), - df.Point(0.3, -0.25), - df.Point(0.65, -0.6), - df.Point(1.0, -0.7)] - -# subdomain3 -subdomain3_vertices = [ - interface34_vertices[0], - interface34_vertices[1], - interface34_vertices[2], - interface34_vertices[3], - interface34_vertices[4], # southern boundary, 34 interface - subdomain2_vertices[5], # eastern boundary, outer boundary - subdomain2_vertices[4], - subdomain2_vertices[3], - subdomain2_vertices[2], - subdomain2_vertices[1], - subdomain2_vertices[0] # northern boundary, 23 interface - ] - -subdomain3_outer_boundary_verts = { - 0: [interface34_vertices[4], - subdomain2_vertices[5]], - 1: [subdomain2_vertices[0], - interface34_vertices[0]] -} - -# subdomain4 -subdomain4_vertices = [ - subdomain0_vertices[0], - subdomain0_vertices[1], # southern boundary, outer boundary - subdomain3_vertices[4], # eastern boundary, outer boundary - subdomain3_vertices[3], - subdomain3_vertices[2], - subdomain3_vertices[1], - subdomain3_vertices[0] - ] # northern boundary, 34 interface - - -subdomain4_outer_boundary_verts = { - 0: [subdomain4_vertices[6], - subdomain4_vertices[0], - subdomain4_vertices[1], - subdomain4_vertices[2]] -} - - -subdomain_def_points = [ - subdomain0_vertices, - subdomain1_vertices, - subdomain2_vertices, - subdomain3_vertices, - subdomain4_vertices - ] - -# interface_vertices introduces a global numbering of interfaces. -interface_def_points = [ - interface12_vertices, interface23_vertices, interface34_vertices - ] - -adjacent_subdomains = [[1, 2], [2, 3], [3, 4]] - -# 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 -} - +substructuring = dss.layeredSoil() +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 = { @@ -665,17 +513,11 @@ p_e_sym = { # } # } -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} @@ -701,35 +543,17 @@ exact_solution = exact_solution_example['exact_solution'] initial_condition = exact_solution_example['initial_condition'] # BOUNDARY CONDITIONS ######################################################### -# Dictionary of dirichlet boundary conditions. -dirichletBC = dict() -# similarly to the outer boundary dictionary, if a patch has no outer boundary -# None should be written instead of an expression. -# This is a bit of a brainfuck: -# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind. -# Since a domain patch can have several disjoint outer boundary parts, the -# expressions need to get an enumaration index which starts at 0. -# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of -# subdomain ind and boundary part j. -# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting'] -# return the actual expression needed for the dirichlet condition for both -# phases if present. - -# subdomain index: {outer boudary part index: {phase: expression}} -for subdomain in isRichards.keys(): - # 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]} - ) - +# 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 @@ -739,88 +563,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() + # run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py index c41f753444423526f6ce50901918730e81f7ea9a..fd9651e24c0a04f90e65a3f9141d284de79d83b3 100755 --- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py +++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.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 -# 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'): os.mkdir('./output') print("Directory ", './output', " created ") @@ -24,37 +30,40 @@ 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-layered-soil-realistic-g-same-intrinsic-perm" -# name of this very file. Needed for log output. +# The name of this very file. Needed for creating log output. thisfile = "TP-R-layered_soil-g-but-same-perm-coarse-dt-longterm.py" +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep max_iter_num = 1000 FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### mesh_study = False resolutions = { - # 1: 9e-6, # h=2 - # 2: 9e-6, # h=1.1180 - # 4: 9e-6, # h=0.5590 - # 8: 9e-6, # h=0.2814 - # 16: 5e-6, # h=0.1412 + # 1: 1e-6, + # 2: 1e-6, + # 4: 1e-6, + # 8: 1e-5, + # 16: 5e-6, 32: 5e-6, - # 64: 1e-6, - # 128: 1e-6 + # 64: 2e-6, + # 128: 1e-6, + # 256: 1e-6, } -# GRID ####################### -# mesh_resolution = 20 +# starttimes gives a list of starttimes to run the simulation from. +# The list is looped over and a simulation is run with t_0 as initial time +# for each element t_0 in starttimes. +starttimes = [0.0] timestep_size = 0.01 number_of_timesteps = 400 -plot_timestep_every = 1 -# 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 = 5 -starttimes = [0.0] + +# LDD scheme parameters ###################################################### Lw1 = 0.025 # /timestep_size Lnw1 = Lw1 Lw2 = 0.025 # /timestep_size @@ -75,20 +84,38 @@ include_gravity = True debugflag = False analyse_condition = False +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 4 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 5 -output_string = "./output/{}-{}_timesteps{}_P{}".format( - datestr, use_case, number_of_timesteps, FEM_Lagrange_degree - ) - -# toggle what should be written to files +# fine grained control over data to be written to disk in the mesh study case +# as well as for a regular simuation for a fixed grid. if mesh_study: write_to_file = { + # output the relative errornorm (integration in space) w.r.t. an exact + # solution for each timestep into a csv file. 'space_errornorms': True, + # save the mesh and marker functions to disk 'meshes_and_markers': True, + # 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: @@ -102,155 +129,19 @@ else: 'subsequent_errors': True } +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) -# global domain -subdomain0_vertices = [ - df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0) - ] - - -interface12_vertices = [df.Point(-1.0, 0.8), - df.Point(0.3, 0.8), - df.Point(0.5, 0.9), - df.Point(0.8, 0.7), - df.Point(1.0, 0.65)] -# subdomain1. -subdomain1_vertices = [ - interface12_vertices[0], - interface12_vertices[1], - interface12_vertices[2], - interface12_vertices[3], - interface12_vertices[4], # southern boundary, 12 interface - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3] # northern boundary, outer on_boundary - ] - - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[4], - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3], - interface12_vertices[0]] -} - - -# interface23 -interface23_vertices = [ - df.Point(-1.0, 0.0), - df.Point(-0.35, 0.0), - # df.Point(6.5, 4.5), - df.Point(0.0, 0.0), - df.Point(0.5, 0.0), - # df.Point(11.5, 3.5), - # df.Point(13.0, 3) - df.Point(0.85, 0.0), - df.Point(1.0, 0.0) - ] - -# subdomain1 -subdomain2_vertices = [ - interface23_vertices[0], - interface23_vertices[1], - interface23_vertices[2], - interface23_vertices[3], - interface23_vertices[4], - interface23_vertices[5], # southern boundary, 23 interface - subdomain1_vertices[4], # eastern boundary, outer boundary - subdomain1_vertices[3], - subdomain1_vertices[2], - subdomain1_vertices[1], - subdomain1_vertices[0] # northern boundary, 12 interface - ] - -subdomain2_outer_boundary_verts = { - 0: [interface23_vertices[5], - subdomain1_vertices[4]], - 1: [subdomain1_vertices[0], - interface23_vertices[0]] -} - - -# interface34 -interface34_vertices = [df.Point(-1.0, -0.6), - df.Point(-0.6, -0.45), - df.Point(0.3, -0.25), - df.Point(0.65, -0.6), - df.Point(1.0, -0.7)] - -# subdomain3 -subdomain3_vertices = [ - interface34_vertices[0], - interface34_vertices[1], - interface34_vertices[2], - interface34_vertices[3], - interface34_vertices[4], # southern boundary, 34 interface - subdomain2_vertices[5], # eastern boundary, outer boundary - subdomain2_vertices[4], - subdomain2_vertices[3], - subdomain2_vertices[2], - subdomain2_vertices[1], - subdomain2_vertices[0] # northern boundary, 23 interface - ] - -subdomain3_outer_boundary_verts = { - 0: [interface34_vertices[4], - subdomain2_vertices[5]], - 1: [subdomain2_vertices[0], - interface34_vertices[0]] -} - -# subdomain4 -subdomain4_vertices = [ - subdomain0_vertices[0], - subdomain0_vertices[1], # southern boundary, outer boundary - subdomain3_vertices[4], # eastern boundary, outer boundary - subdomain3_vertices[3], - subdomain3_vertices[2], - subdomain3_vertices[1], - subdomain3_vertices[0] - ] # northern boundary, 34 interface - - -subdomain4_outer_boundary_verts = { - 0: [subdomain4_vertices[6], - subdomain4_vertices[0], - subdomain4_vertices[1], - subdomain4_vertices[2]] -} - +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.layeredSoil() +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 -subdomain_def_points = [ - subdomain0_vertices, - subdomain1_vertices, - subdomain2_vertices, - subdomain3_vertices, - subdomain4_vertices - ] - -# interface_vertices introduces a global numbering of interfaces. -interface_def_points = [ - interface12_vertices, interface23_vertices, interface34_vertices - ] - -adjacent_subdomains = [[1, 2], [2, 3], [3, 4]] - -# 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 -} +# MODEL CONFIGURATION ######################################################### isRichards = { 1: True, @@ -638,17 +529,11 @@ p_e_sym = { # } # } -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} @@ -673,35 +558,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 @@ -711,88 +579,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() + # run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py index 8e4878af61865d60b0c0dc83f48415230e9ef5ae..1916d58fd6edf8a3ec64d61a6352fe41d0f76a8f 100755 --- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.py +++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil-g-but-same-perm.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,40 @@ else: date = datetime.datetime.now() datestr = date.strftime("%Y-%m-%d") -# init sympy session -sym.init_printing() -# solver_tol = 6E-7 -use_case = "TP-R-layered-soil-realistic-g-same-intrinsic-perm" -# name of this very file. Needed for log output. + +# Name of the usecase that will be printed during simulation. +use_case = "TP-R-layered-soil-realistic-same-perm" +# The name of this very file. Needed for creating log output. thisfile = "TP-R-layered_soil-g-but-same-perm.py" -max_iter_num = 750 +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep +max_iter_num = 1000 FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### mesh_study = False resolutions = { - # 1: 9e-6, # h=2 - # 2: 9e-6, # h=1.1180 - # 4: 9e-6, # h=0.5590 - # 8: 9e-6, # h=0.2814 - # 16: 5e-6, # h=0.1412 - 32: 4e-6, - # 64: 1e-6, - # 128: 1e-6 + # 1: 1e-6, + # 2: 1e-6, + # 4: 1e-6, + # 8: 1e-5, + # 16: 5e-6, + 32: 5e-6, + # 64: 2e-6, + # 128: 1e-6, + # 256: 1e-6, } -# GRID ####################### -# mesh_resolution = 20 -timestep_size = 0.001 -number_of_timesteps = 1000 -plot_timestep_every = 3 -# 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 = 5 +# 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 = 800 + +# LDD scheme parameters ###################################################### Lw1 = 0.025 # /timestep_size Lnw1 = Lw1 Lw2 = 0.025 # /timestep_size @@ -75,20 +84,38 @@ include_gravity = True debugflag = False analyse_condition = False +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 4 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 5 -output_string = "./output/{}-{}_timesteps{}_P{}".format( - datestr, use_case, number_of_timesteps, FEM_Lagrange_degree - ) - -# toggle what should be written to files +# fine grained control over data to be written to disk in the mesh study case +# as well as for a regular simuation for a fixed grid. if mesh_study: write_to_file = { + # output the relative errornorm (integration in space) w.r.t. an exact + # solution for each timestep into a csv file. 'space_errornorms': True, + # save the mesh and marker functions to disk 'meshes_and_markers': True, + # 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: @@ -102,156 +129,19 @@ else: 'subsequent_errors': True } +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) -# global domain -subdomain0_vertices = [ - df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0) - ] - - -interface12_vertices = [df.Point(-1.0, 0.8), - df.Point(0.3, 0.8), - df.Point(0.5, 0.9), - df.Point(0.8, 0.7), - df.Point(1.0, 0.65)] -# subdomain1. -subdomain1_vertices = [ - interface12_vertices[0], - interface12_vertices[1], - interface12_vertices[2], - interface12_vertices[3], - interface12_vertices[4], # southern boundary, 12 interface - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3] # northern boundary, outer on_boundary - ] - - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[4], - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3], - interface12_vertices[0]] -} - - -# interface23 -interface23_vertices = [ - df.Point(-1.0, 0.0), - df.Point(-0.35, 0.0), - # df.Point(6.5, 4.5), - df.Point(0.0, 0.0), - df.Point(0.5, 0.0), - # df.Point(11.5, 3.5), - # df.Point(13.0, 3) - df.Point(0.85, 0.0), - df.Point(1.0, 0.0) - ] - -# subdomain1 -subdomain2_vertices = [ - interface23_vertices[0], - interface23_vertices[1], - interface23_vertices[2], - interface23_vertices[3], - interface23_vertices[4], - interface23_vertices[5], # southern boundary, 23 interface - subdomain1_vertices[4], # eastern boundary, outer boundary - subdomain1_vertices[3], - subdomain1_vertices[2], - subdomain1_vertices[1], - subdomain1_vertices[0] # northern boundary, 12 interface - ] - -subdomain2_outer_boundary_verts = { - 0: [interface23_vertices[5], - subdomain1_vertices[4]], - 1: [subdomain1_vertices[0], - interface23_vertices[0]] -} - - -# interface34 -interface34_vertices = [df.Point(-1.0, -0.6), - df.Point(-0.6, -0.45), - df.Point(0.3, -0.25), - df.Point(0.65, -0.6), - df.Point(1.0, -0.7)] - -# subdomain3 -subdomain3_vertices = [ - interface34_vertices[0], - interface34_vertices[1], - interface34_vertices[2], - interface34_vertices[3], - interface34_vertices[4], # southern boundary, 34 interface - subdomain2_vertices[5], # eastern boundary, outer boundary - subdomain2_vertices[4], - subdomain2_vertices[3], - subdomain2_vertices[2], - subdomain2_vertices[1], - subdomain2_vertices[0] # northern boundary, 23 interface - ] - -subdomain3_outer_boundary_verts = { - 0: [interface34_vertices[4], - subdomain2_vertices[5]], - 1: [subdomain2_vertices[0], - interface34_vertices[0]] -} - -# subdomain4 -subdomain4_vertices = [ - subdomain0_vertices[0], - subdomain0_vertices[1], # southern boundary, outer boundary - subdomain3_vertices[4], # eastern boundary, outer boundary - subdomain3_vertices[3], - subdomain3_vertices[2], - subdomain3_vertices[1], - subdomain3_vertices[0] - ] # northern boundary, 34 interface - - -subdomain4_outer_boundary_verts = { - 0: [subdomain4_vertices[6], - subdomain4_vertices[0], - subdomain4_vertices[1], - subdomain4_vertices[2]] -} - - -subdomain_def_points = [ - subdomain0_vertices, - subdomain1_vertices, - subdomain2_vertices, - subdomain3_vertices, - subdomain4_vertices - ] - -# interface_vertices introduces a global numbering of interfaces. -interface_def_points = [ - interface12_vertices, interface23_vertices, interface34_vertices - ] - -adjacent_subdomains = [[1, 2], [2, 3], [3, 4]] - -# 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 -} +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.layeredSoil() +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 = { 1: True, 2: True, @@ -638,17 +528,11 @@ p_e_sym = { # } # } -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} @@ -673,35 +557,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 @@ -711,88 +578,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() + # run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py index 49141a97ba720dc323f76c9f28e4a35df115f8c2..c115fa15ed2638fa20b14bf6e466e8008c68ded7 100755 --- a/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py +++ b/Two-phase-Richards/multi-patch/layered_soil/TP-R-layered_soil.py @@ -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() @@ -37,7 +38,7 @@ thisfile = "TP-R-layered_soil.py" # GENERAL SOLVER CONFIG ###################################################### # maximal iteration per timestep -max_iter_num = 300 +max_iter_num = 10 FEM_Lagrange_degree = 1 # GRID AND MESH STUDY SPECIFICATIONS ######################################### @@ -46,10 +47,10 @@ resolutions = { # 1: 1e-6, # 2: 1e-6, # 4: 1e-6, - # 8: 1e-6, + # 8: 1e-5, # 16: 5e-6, - # 32: 5e-6, - 64: 2e-6, + 32: 5e-6, + # 64: 2e-6, # 128: 1e-6, # 256: 1e-6, } @@ -59,7 +60,7 @@ resolutions = { # for each element t_0 in starttimes. starttimes = [0.0] timestep_size = 0.001 -number_of_timesteps = 20 +number_of_timesteps = 5 # LDD scheme parameters ###################################################### @@ -80,17 +81,17 @@ lambda34_w = 40 lambda34_nw = 40 include_gravity = False -debugflag = False -analyse_condition = True +debugflag = True +analyse_condition = False # I/O CONFIG ################################################################# # when number_of_timesteps is high, it might take a long time to write all # timesteps to disk. Therefore, you can choose to only write data of every # plot_timestep_every timestep to disk. -plot_timestep_every = 4 +plot_timestep_every = 1 # Decide how many timesteps you want analysed. Analysed means, that # subsequent errors of the L-iteration within the timestep are written out. -number_of_timesteps_to_analyse = 5 +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. @@ -129,168 +130,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 - ) - +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) # DOMAIN AND INTERFACE ####################################################### -# global domain -subdomain0_vertices = [ - df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0) - ] - - -interface12_vertices = [df.Point(-1.0, 0.8), - df.Point(0.3, 0.8), - df.Point(0.5, 0.9), - df.Point(0.8, 0.7), - df.Point(1.0, 0.65)] -# subdomain1. -subdomain1_vertices = [ - interface12_vertices[0], - interface12_vertices[1], - interface12_vertices[2], - interface12_vertices[3], - interface12_vertices[4], # southern boundary, 12 interface - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3] # northern boundary, outer on_boundary - ] - - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[4], - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3], - interface12_vertices[0]] -} - - -# interface23 -interface23_vertices = [ - df.Point(-1.0, 0.0), - df.Point(-0.35, 0.0), - # df.Point(6.5, 4.5), - df.Point(0.0, 0.0), - df.Point(0.5, 0.0), - # df.Point(11.5, 3.5), - # df.Point(13.0, 3) - df.Point(0.85, 0.0), - df.Point(1.0, 0.0) - ] - -# subdomain1 -subdomain2_vertices = [ - interface23_vertices[0], - interface23_vertices[1], - interface23_vertices[2], - interface23_vertices[3], - interface23_vertices[4], - interface23_vertices[5], # southern boundary, 23 interface - subdomain1_vertices[4], # eastern boundary, outer boundary - subdomain1_vertices[3], - subdomain1_vertices[2], - subdomain1_vertices[1], - subdomain1_vertices[0] # northern boundary, 12 interface - ] - -subdomain2_outer_boundary_verts = { - 0: [interface23_vertices[5], - subdomain1_vertices[4]], - 1: [subdomain1_vertices[0], - interface23_vertices[0]] -} - - -# interface34 -interface34_vertices = [df.Point(-1.0, -0.6), - df.Point(-0.6, -0.45), - df.Point(0.3, -0.25), - df.Point(0.65, -0.6), - df.Point(1.0, -0.7)] - -# subdomain3 -subdomain3_vertices = [ - interface34_vertices[0], - interface34_vertices[1], - interface34_vertices[2], - interface34_vertices[3], - interface34_vertices[4], # southern boundary, 34 interface - subdomain2_vertices[5], # eastern boundary, outer boundary - subdomain2_vertices[4], - subdomain2_vertices[3], - subdomain2_vertices[2], - subdomain2_vertices[1], - subdomain2_vertices[0] # northern boundary, 23 interface - ] - -subdomain3_outer_boundary_verts = { - 0: [interface34_vertices[4], - subdomain2_vertices[5]], - 1: [subdomain2_vertices[0], - interface34_vertices[0]] -} - -# subdomain4 -subdomain4_vertices = [ - subdomain0_vertices[0], - subdomain0_vertices[1], # southern boundary, outer boundary - subdomain3_vertices[4], # eastern boundary, outer boundary - subdomain3_vertices[3], - subdomain3_vertices[2], - subdomain3_vertices[1], - subdomain3_vertices[0] - ] # northern boundary, 34 interface - - -subdomain4_outer_boundary_verts = { - 0: [subdomain4_vertices[6], - subdomain4_vertices[0], - subdomain4_vertices[1], - subdomain4_vertices[2]] -} - - -subdomain_def_points = [ - subdomain0_vertices, - subdomain1_vertices, - subdomain2_vertices, - subdomain3_vertices, - subdomain4_vertices - ] - -# interface_vertices introduces a global numbering of interfaces. -interface_def_points = [ - interface12_vertices, interface23_vertices, interface34_vertices - ] - -adjacent_subdomains = [[1, 2], [2, 3], [3, 4]] - -# 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 -} - +substructuring = dss.layeredSoil() +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 = { @@ -676,17 +525,11 @@ p_e_sym = { # } # } -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} @@ -712,35 +555,17 @@ exact_solution = exact_solution_example['exact_solution'] initial_condition = exact_solution_example['initial_condition'] # BOUNDARY CONDITIONS ######################################################### -# Dictionary of dirichlet boundary conditions. -dirichletBC = dict() -# similarly to the outer boundary dictionary, if a patch has no outer boundary -# None should be written instead of an expression. -# This is a bit of a brainfuck: -# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind. -# Since a domain patch can have several disjoint outer boundary parts, the -# expressions need to get an enumaration index which starts at 0. -# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of -# subdomain ind and boundary part j. -# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting'] -# return the actual expression needed for the dirichlet condition for both -# phases if present. - -# subdomain index: {outer boudary part index: {phase: expression}} -for subdomain in isRichards.keys(): - # 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]} - ) - +# 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 @@ -750,88 +575,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() + # run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.py index 1e2698aa49f173e801c891a5d9f33b3a9584f4de..3c14df1e9159379de770ec5f8c9813516acc9aca 100755 --- a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.py +++ b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-g-but-same-perm-mesh-study.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,15 +30,17 @@ 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-layered-soil-realistic-g-same-intrinsic-perm" -# name of this very file. Needed for log output. +# The name of this very file. Needed for creating log output. thisfile = "TP-R-layered_soil-g-but-same-perm-mesh-study.py" +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep max_iter_num = 750 FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### mesh_study = True resolutions = { 1: 9e-6, # h=2 @@ -45,16 +53,15 @@ resolutions = { 128: 1e-6 } -# GRID ####################### -# mesh_resolution = 20 +# starttimes gives a list of starttimes to run the simulation from. +# The list is looped over and a simulation is run with t_0 as initial time +# for each element t_0 in starttimes. +starttimes = [0.0] timestep_size = 0.001 number_of_timesteps = 1000 -plot_timestep_every = 4 -# 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 = 5 -starttimes = [0.0] + +# LDD scheme parameters ###################################################### Lw1 = 0.025 # /timestep_size Lnw1 = Lw1 Lw2 = 0.025 # /timestep_size @@ -75,20 +82,38 @@ include_gravity = True debugflag = False analyse_condition = False +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 4 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 5 -output_string = "./output/{}-{}_timesteps{}_P{}".format( - datestr, use_case, number_of_timesteps, FEM_Lagrange_degree - ) - -# toggle what should be written to files +# fine grained control over data to be written to disk in the mesh study case +# as well as for a regular simuation for a fixed grid. if mesh_study: write_to_file = { + # output the relative errornorm (integration in space) w.r.t. an exact + # solution for each timestep into a csv file. 'space_errornorms': True, + # save the mesh and marker functions to disk 'meshes_and_markers': True, + # 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: @@ -102,156 +127,19 @@ else: 'subsequent_errors': True } +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) -# global domain -subdomain0_vertices = [ - df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0) - ] - - -interface12_vertices = [df.Point(-1.0, 0.8), - df.Point(0.3, 0.8), - df.Point(0.5, 0.9), - df.Point(0.8, 0.7), - df.Point(1.0, 0.65)] -# subdomain1. -subdomain1_vertices = [ - interface12_vertices[0], - interface12_vertices[1], - interface12_vertices[2], - interface12_vertices[3], - interface12_vertices[4], # southern boundary, 12 interface - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3] # northern boundary, outer on_boundary - ] - - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[4], - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3], - interface12_vertices[0]] -} - - -# interface23 -interface23_vertices = [ - df.Point(-1.0, 0.0), - df.Point(-0.35, 0.0), - # df.Point(6.5, 4.5), - df.Point(0.0, 0.0), - df.Point(0.5, 0.0), - # df.Point(11.5, 3.5), - # df.Point(13.0, 3) - df.Point(0.85, 0.0), - df.Point(1.0, 0.0) - ] - -# subdomain1 -subdomain2_vertices = [ - interface23_vertices[0], - interface23_vertices[1], - interface23_vertices[2], - interface23_vertices[3], - interface23_vertices[4], - interface23_vertices[5], # southern boundary, 23 interface - subdomain1_vertices[4], # eastern boundary, outer boundary - subdomain1_vertices[3], - subdomain1_vertices[2], - subdomain1_vertices[1], - subdomain1_vertices[0] # northern boundary, 12 interface - ] - -subdomain2_outer_boundary_verts = { - 0: [interface23_vertices[5], - subdomain1_vertices[4]], - 1: [subdomain1_vertices[0], - interface23_vertices[0]] -} - - -# interface34 -interface34_vertices = [df.Point(-1.0, -0.6), - df.Point(-0.6, -0.45), - df.Point(0.3, -0.25), - df.Point(0.65, -0.6), - df.Point(1.0, -0.7)] - -# subdomain3 -subdomain3_vertices = [ - interface34_vertices[0], - interface34_vertices[1], - interface34_vertices[2], - interface34_vertices[3], - interface34_vertices[4], # southern boundary, 34 interface - subdomain2_vertices[5], # eastern boundary, outer boundary - subdomain2_vertices[4], - subdomain2_vertices[3], - subdomain2_vertices[2], - subdomain2_vertices[1], - subdomain2_vertices[0] # northern boundary, 23 interface - ] - -subdomain3_outer_boundary_verts = { - 0: [interface34_vertices[4], - subdomain2_vertices[5]], - 1: [subdomain2_vertices[0], - interface34_vertices[0]] -} - -# subdomain4 -subdomain4_vertices = [ - subdomain0_vertices[0], - subdomain0_vertices[1], # southern boundary, outer boundary - subdomain3_vertices[4], # eastern boundary, outer boundary - subdomain3_vertices[3], - subdomain3_vertices[2], - subdomain3_vertices[1], - subdomain3_vertices[0] - ] # northern boundary, 34 interface - - -subdomain4_outer_boundary_verts = { - 0: [subdomain4_vertices[6], - subdomain4_vertices[0], - subdomain4_vertices[1], - subdomain4_vertices[2]] -} - - -subdomain_def_points = [ - subdomain0_vertices, - subdomain1_vertices, - subdomain2_vertices, - subdomain3_vertices, - subdomain4_vertices - ] - -# interface_vertices introduces a global numbering of interfaces. -interface_def_points = [ - interface12_vertices, interface23_vertices, interface34_vertices - ] - -adjacent_subdomains = [[1, 2], [2, 3], [3, 4]] - -# 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 -} +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.layeredSoil() +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 = { 1: True, 2: True, @@ -638,17 +526,11 @@ p_e_sym = { # } # } -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} @@ -673,35 +555,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 @@ -711,88 +576,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() + # run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.py b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.py index 3269d5cc37412c5b68b768a6af14d04fac051740..bee65c9a8806119abdf2e3ec4c193bf51b3d52c8 100755 --- a/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.py +++ b/Two-phase-Richards/multi-patch/layered_soil/mesh_study/TP-R-layered_soil-mesh-study.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 -# 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'): os.mkdir('./output') print("Directory ", './output', " created ") @@ -24,45 +30,46 @@ 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-layered-soil-realistic-same-intrinsic-perm" -# name of this very file. Needed for log output. +# The name of this very file. Needed for creating log output. thisfile = "TP-R-layered_soil-mesh-study.py" +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep max_iter_num = 750 FEM_Lagrange_degree = 1 -mesh_study = False + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### +mesh_study = True resolutions = { - # 1: 9e-6, # h=2 - # 2: 9e-6, # h=1.1180 - # 4: 9e-6, # h=0.5590 - # 8: 9e-6, # h=0.2814 - 16: 2e-6, # h=0.1412 - # 32: 2e-6, - # 64: 2e-6, - # 128: 2e-6 + 1: 9e-6, # h=2 + 2: 9e-6, # h=1.1180 + 4: 9e-6, # h=0.5590 + 8: 9e-6, # h=0.2814 + 16: 5e-6, # h=0.1412 + 32: 4e-6, + 64: 1e-6, + 128: 1e-6 } -# GRID ####################### -# mesh_resolution = 20 -timestep_size = 0.001 -number_of_timesteps = 5 -plot_timestep_every = 1 -# 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 = 5 +# 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 = 1000 + +# LDD scheme parameters ###################################################### Lw1 = 0.025 # /timestep_size -Lnw1 = 0.025 +Lnw1 = Lw1 Lw2 = 0.025 # /timestep_size -Lnw2 = 0.025 +Lnw2 = Lw2 Lw3 = 0.025 # /timestep_size -Lnw3 = 0.025 +Lnw3 = Lw3 Lw4 = 0.025 # /timestep_size -Lnw4 = 0.025 +Lnw4 = Lw4 lambda12_w = 4 lambda12_nw = 4 @@ -72,23 +79,41 @@ lambda34_w = 4 lambda34_nw = 4 include_gravity = False -debugflag = True +debugflag = False analyse_condition = False +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 4 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 5 -output_string = "./output/{}-{}_timesteps{}_P{}".format( - datestr, use_case, number_of_timesteps, FEM_Lagrange_degree - ) - -# toggle what should be written to files +# fine grained control over data to be written to disk in the mesh study case +# as well as for a regular simuation for a fixed grid. if mesh_study: write_to_file = { + # output the relative errornorm (integration in space) w.r.t. an exact + # solution for each timestep into a csv file. 'space_errornorms': True, + # save the mesh and marker functions to disk 'meshes_and_markers': True, + # 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: @@ -102,156 +127,19 @@ else: 'subsequent_errors': True } +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) -# global domain -subdomain0_vertices = [ - df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0) - ] - - -interface12_vertices = [df.Point(-1.0, 0.8), - df.Point(0.3, 0.8), - df.Point(0.5, 0.9), - df.Point(0.8, 0.7), - df.Point(1.0, 0.65)] -# subdomain1. -subdomain1_vertices = [ - interface12_vertices[0], - interface12_vertices[1], - interface12_vertices[2], - interface12_vertices[3], - interface12_vertices[4], # southern boundary, 12 interface - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3] # northern boundary, outer on_boundary - ] - - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[4], - subdomain0_vertices[2], # eastern boundary, outer boundary - subdomain0_vertices[3], - interface12_vertices[0]] -} - - -# interface23 -interface23_vertices = [ - df.Point(-1.0, 0.0), - df.Point(-0.35, 0.0), - # df.Point(6.5, 4.5), - df.Point(0.0, 0.0), - df.Point(0.5, 0.0), - # df.Point(11.5, 3.5), - # df.Point(13.0, 3) - df.Point(0.85, 0.0), - df.Point(1.0, 0.0) - ] - -# subdomain1 -subdomain2_vertices = [ - interface23_vertices[0], - interface23_vertices[1], - interface23_vertices[2], - interface23_vertices[3], - interface23_vertices[4], - interface23_vertices[5], # southern boundary, 23 interface - subdomain1_vertices[4], # eastern boundary, outer boundary - subdomain1_vertices[3], - subdomain1_vertices[2], - subdomain1_vertices[1], - subdomain1_vertices[0] # northern boundary, 12 interface - ] - -subdomain2_outer_boundary_verts = { - 0: [interface23_vertices[5], - subdomain1_vertices[4]], - 1: [subdomain1_vertices[0], - interface23_vertices[0]] -} - - -# interface34 -interface34_vertices = [df.Point(-1.0, -0.6), - df.Point(-0.6, -0.45), - df.Point(0.3, -0.25), - df.Point(0.65, -0.6), - df.Point(1.0, -0.7)] - -# subdomain3 -subdomain3_vertices = [ - interface34_vertices[0], - interface34_vertices[1], - interface34_vertices[2], - interface34_vertices[3], - interface34_vertices[4], # southern boundary, 34 interface - subdomain2_vertices[5], # eastern boundary, outer boundary - subdomain2_vertices[4], - subdomain2_vertices[3], - subdomain2_vertices[2], - subdomain2_vertices[1], - subdomain2_vertices[0] # northern boundary, 23 interface - ] - -subdomain3_outer_boundary_verts = { - 0: [interface34_vertices[4], - subdomain2_vertices[5]], - 1: [subdomain2_vertices[0], - interface34_vertices[0]] -} - -# subdomain4 -subdomain4_vertices = [ - subdomain0_vertices[0], - subdomain0_vertices[1], # southern boundary, outer boundary - subdomain3_vertices[4], # eastern boundary, outer boundary - subdomain3_vertices[3], - subdomain3_vertices[2], - subdomain3_vertices[1], - subdomain3_vertices[0] - ] # northern boundary, 34 interface - - -subdomain4_outer_boundary_verts = { - 0: [subdomain4_vertices[6], - subdomain4_vertices[0], - subdomain4_vertices[1], - subdomain4_vertices[2]] -} - - -subdomain_def_points = [ - subdomain0_vertices, - subdomain1_vertices, - subdomain2_vertices, - subdomain3_vertices, - subdomain4_vertices - ] - -# interface_vertices introduces a global numbering of interfaces. -interface_def_points = [ - interface12_vertices, interface23_vertices, interface34_vertices - ] - -adjacent_subdomains = [[1, 2], [2, 3], [3, 4]] - -# 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 -} +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.layeredSoil() +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 = { 1: True, 2: True, @@ -638,17 +526,11 @@ p_e_sym = { # } # } -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} @@ -673,35 +555,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 @@ -711,88 +576,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() + # run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # )