diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-realistic-parameters-densities-scaled.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/Archive/TP-R-2-realistic-parameters-densities-scaled.py similarity index 100% rename from Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-realistic-parameters-densities-scaled.py rename to Two-phase-Richards/two-patch/TP-R-two-patch-test-case/Archive/TP-R-2-realistic-parameters-densities-scaled.py diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py index 9ae8b6b42002422144213fd485e605f8d890bd34..d0d78747e81577d7178e04fce5e472cf1a19b1a2 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py @@ -1,147 +1,138 @@ -#!/usr/bin/python3 +"""TPR 2 patch soil simulation. + +This program sets up an LDD simulation +""" + import dolfin as df -import mshr -import numpy as np import sympy as sym -import typing as tp -import domainPatch as dp -import LDDsimulation as ldd import functools as ft +import LDDsimulation as ldd import helpers as hlp import datetime import os import pandas as pd +import multiprocessing as mp +import domainSubstructuring as dss + +# init sympy session +sym.init_printing() + +# PREREQUISITS ############################################################### +# check if output directory "./output" exists. This will be used in +# the generation of the output string. +if not os.path.exists('./output'): + os.mkdir('./output') + print("Directory ", './output', " created ") +else: + print("Directory ", './output', " already exists. Will use as output \ + directory") date = datetime.datetime.now() datestr = date.strftime("%Y-%m-%d") -#import ufl as ufl -# init sympy session -sym.init_printing() +# Name of the usecase that will be printed during simulation. +use_case = "TPR-2-patch-realistic-pure-dd" +# The name of this very file. Needed for creating log output. +thisfile = "TP-R-2-patch-pure-dd-realistic.py" -use_case = "TP-R-2-patch-realistic-pure-dd" -# solver_tol = 6E-7 -max_iter_num = 1000 +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep +max_iter_num = 500 FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### mesh_study = False resolutions = { - # 1: 7e-7, - # 2: 7e-7, - # 4: 7e-7, - # 8: 7e-7, - # 16: 7e-7, - 32: 1e-6, - # 64: 7e-7, - # 128: 7e-7, - # 256: 7e-7 + # 1: 1e-5, + # 2: 1e-5, + # 4: 1e-5, + # 8: 1e-5, + 16: 5e-6, + # 32: 5e-6, + # 64: 2e-6, + # 128: 2e-6, + # 256: 1e-6, } -############ GRID ####################### -timestep_size = 0.00001 -number_of_timesteps = 20 -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 = 4 +# 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.00001 +number_of_timesteps = 5 -Lw = 0.5 #/timestep_size -Lnw= 0.5 +# LDD scheme parameters ###################################################### +Lw1 = 0.5 #/timestep_size +Lnw1= 0.5 + +Lw2 = 0.5 #/timestep_size +Lnw2= 0.5 lambda_w = 1 lambda_nw = 1 include_gravity = True -debugflag = True +debugflag = False analyse_condition = False -output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree) +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 1 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 1 -# 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, - 'solutions': False, - 'absolute_differences': 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, - 'subsequent_errors': False + # 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: write_to_file = { 'space_errornorms': True, 'meshes_and_markers': True, - 'L_iterations_per_timestep': True, + 'L_iterations_per_timestep': False, 'solutions': True, 'absolute_differences': True, 'condition_numbers': analyse_condition, 'subsequent_errors': True } +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) -##### Domain and Interface #### -# global simulation domain domain -sub_domain0_vertices = [df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0)] -# interface between subdomain1 and subdomain2 -interface12_vertices = [df.Point(-1.0, 0.0), - df.Point(1.0, 0.0) ] -# subdomain1. -sub_domain1_vertices = [interface12_vertices[0], - interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3]] - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for defining -# the Dirichlet boundary conditions. If a domain is completely internal, the -# dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[1], # - sub_domain0_vertices[2], - sub_domain0_vertices[3], # - interface12_vertices[0]] -} -# subdomain2 -sub_domain2_vertices = [sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1], - interface12_vertices[0] ] - -subdomain2_outer_boundary_verts = { - 0: [interface12_vertices[0], # - sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1]] -} - -# list of subdomains given by the boundary polygon vertices. -# Subdomains are given as a list of dolfin points forming -# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used -# to create the subdomain. subdomain_def_points[0] contains the -# vertices of the global simulation domain and subdomain_def_points[i] contains the -# vertices of the subdomain i. -subdomain_def_points = [sub_domain0_vertices,# - sub_domain1_vertices,# - sub_domain2_vertices] -# in the below list, index 0 corresponds to the 12 interface which has index 1 -interface_def_points = [interface12_vertices] - -# if a subdomain has no outer boundary write None instead, i.e. -# i: None -# if i is the index of the inner subdomain. -outer_boundary_def_points = { - # subdomain number - 1 : subdomain1_outer_boundary_verts, - 2 : subdomain2_outer_boundary_verts -} +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points -# adjacent_subdomains[i] contains the indices of the subdomains sharing the -# interface i (i.e. given by interface_def_points[i]). -adjacent_subdomains = [[1,2]] +# MODEL CONFIGURATION ######################################################### isRichards = { 1: True, # 2: False @@ -387,14 +378,10 @@ p_e_sym = { 'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, } -pc_e_sym = dict() -for subdomain, isR in isRichards.items(): - if isR: - pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) - else: - pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - - p_e_sym[subdomain]['wetting'].copy()}) - +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) symbols = {"x": x, "y": y, @@ -420,34 +407,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 @@ -457,88 +428,63 @@ print(f.read()) f.close() -# RUN ######################################################################### -for starttime in starttimes: - for mesh_resolution, solver_tol in resolutions.items(): - # initialise LDD simulation class - simulation = ldd.LDDsimulation( - tol=1E-14, - LDDsolver_tol=solver_tol, - debug=debugflag, - max_iter_num=max_iter_num, - FEM_Lagrange_degree=FEM_Lagrange_degree, - mesh_study=mesh_study - ) - - simulation.set_parameters( - use_case=use_case, - output_dir=output_string, - subdomain_def_points=subdomain_def_points, - isRichards=isRichards, - interface_def_points=interface_def_points, - outer_boundary_def_points=outer_boundary_def_points, - adjacent_subdomains=adjacent_subdomains, - mesh_resolution=mesh_resolution, - viscosity=viscosity, - porosity=porosity, - L=L, - lambda_param=lambda_param, - relative_permeability=relative_permeability, - saturation=sat_pressure_relationship, - starttime=starttime, - number_of_timesteps=number_of_timesteps, - number_of_timesteps_to_analyse=number_of_timesteps_to_analyse, - plot_timestep_every=plot_timestep_every, - timestep_size=timestep_size, - sources=source_expression, - initial_conditions=initial_condition, - dirichletBC_expression_strings=dirichletBC, - exact_solution=exact_solution, - densities=densities, - include_gravity=include_gravity, - gravity_acceleration=gravity_acceleration, - write2file=write_to_file, - ) - - simulation.initialise() - output_dir = simulation.output_dir - # simulation.write_exact_solution_to_xdmf() - output = simulation.run(analyse_condition=analyse_condition) - for subdomain_index, subdomain_output in output.items(): - mesh_h = subdomain_output['mesh_size'] - for phase, error_dict in subdomain_output['errornorm'].items(): - filename = output_dir \ - + "subdomain{}".format(subdomain_index)\ - + "-space-time-errornorm-{}-phase.csv".format(phase) - # for errortype, errornorm in error_dict.items(): - - # eocfile = open("eoc_filename", "a") - # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) - # eocfile.close() - # if subdomain.isRichards:mesh_h - data_dict = { - 'mesh_parameter': mesh_resolution, - 'mesh_h': mesh_h, - } - for norm_type, errornorm in error_dict.items(): - data_dict.update( - {norm_type: errornorm} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) is True: - with open(filename, 'a') as f: - errors.to_csv( - f, - header=False, - sep='\t', - encoding='utf-8', - index=False +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": output_string, + "subdomain_def_points": subdomain_def_points, + "isRichards": isRichards, + "interface_def_points": interface_def_points, + "outer_boundary_def_points": outer_boundary_def_points, + "adjacent_subdomains": adjacent_subdomains, + # "mesh_resolution": mesh_resolution, + "viscosity": viscosity, + "porosity": porosity, + "L": L, + "lambda_param": lambda_param, + "relative_permeability": relative_permeability, + "sat_pressure_relationship": sat_pressure_relationship, + # "starttime": starttime, + "number_of_timesteps": number_of_timesteps, + "number_of_timesteps_to_analyse": number_of_timesteps_to_analyse, + "plot_timestep_every": plot_timestep_every, + "timestep_size": timestep_size, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for starttime in starttimes: + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + starttime, + mesh_resolution ) - else: - errors.to_csv( - filename, - sep='\t', - encoding='utf-8', - index=False ) + LDDsim.start() + LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py index 2a42078bdd12a33967010e246814fd2660292fbb..7e26fbfc5c077d05f0d79255d3291cb91764a39d 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py @@ -1,73 +1,113 @@ -#!/usr/bin/python3 +"""TPR 2 patch soil simulation. + +This program sets up an LDD simulation +""" + import dolfin as df -import mshr -import numpy as np import sympy as sym -import typing as tp -import domainPatch as dp -import LDDsimulation as ldd import functools as ft +import LDDsimulation as ldd import helpers as hlp import datetime import os import pandas as pd +import multiprocessing as mp +import domainSubstructuring as dss + +# init sympy session +sym.init_printing() + +# PREREQUISITS ############################################################### +# check if output directory "./output" exists. This will be used in +# the generation of the output string. +if not os.path.exists('./output'): + os.mkdir('./output') + print("Directory ", './output', " created ") +else: + print("Directory ", './output', " already exists. Will use as output \ + directory") date = datetime.datetime.now() datestr = date.strftime("%Y-%m-%d") -#import ufl as ufl -# init sympy session -sym.init_printing() +# Name of the usecase that will be printed during simulation. +use_case = "TPR-2-patch-params-one-pure-dd" +# The name of this very file. Needed for creating log output. +thisfile = "TP-R-2-patch-pure-dd.py" -use_case = "TP-R-2-patch-pure-dd" -# solver_tol = 6E-7 -max_iter_num = 1000 +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep +max_iter_num = 500 FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### mesh_study = False resolutions = { - # 1: 7e-7, - # 2: 7e-7, - # 4: 7e-7, - # 8: 7e-7, - # 16: 5e-7, - 32: 1e-7, - # 64: 7e-7, - # 128: 7e-7, - # 256: 7e-7 + # 1: 1e-5, + # 2: 1e-5, + # 4: 1e-5, + # 8: 1e-5, + 16: 5e-6, + # 32: 5e-6, + # 64: 2e-6, + # 128: 2e-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.005 -number_of_timesteps = 200 +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 = 0 -starttimes = [0.0] +# LDD scheme parameters ###################################################### +Lw1 = 0.05 #/timestep_size +Lnw1= 0.05 -Lw = 0.05 #/timestep_size -Lnw= 0.05 +Lw2 = 0.05 #/timestep_size +Lnw2= 0.05 lambda_w = 20 lambda_nw = 20 + include_gravity = True debugflag = False -analyse_condition = True +analyse_condition = False -output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree) +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 1 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 1 -# 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, - 'solutions': False, - 'absolute_differences': 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, - 'subsequent_errors': False + # 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: write_to_file = { @@ -80,69 +120,154 @@ else: 'subsequent_errors': True } +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) -##### Domain and Interface #### -# global simulation domain domain -sub_domain0_vertices = [df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0)] -# interface between subdomain1 and subdomain2 -interface12_vertices = [df.Point(-1.0, 0.0), - df.Point(1.0, 0.0) ] -# subdomain1. -sub_domain1_vertices = [interface12_vertices[0], - interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3]] - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for defining -# the Dirichlet boundary conditions. If a domain is completely internal, the -# dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[1], # - sub_domain0_vertices[2], - sub_domain0_vertices[3], # - interface12_vertices[0]] -} -# subdomain2 -sub_domain2_vertices = [sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1], - interface12_vertices[0] ] - -subdomain2_outer_boundary_verts = { - 0: [interface12_vertices[0], # - sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1]] -} +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points -# list of subdomains given by the boundary polygon vertices. -# Subdomains are given as a list of dolfin points forming -# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used -# to create the subdomain. subdomain_def_points[0] contains the -# vertices of the global simulation domain and subdomain_def_points[i] contains the -# vertices of the subdomain i. -subdomain_def_points = [sub_domain0_vertices,# - sub_domain1_vertices,# - sub_domain2_vertices] -# in the below list, index 0 corresponds to the 12 interface which has index 1 -interface_def_points = [interface12_vertices] - -# if a subdomain has no outer boundary write None instead, i.e. -# i: None -# if i is the index of the inner subdomain. -outer_boundary_def_points = { - # subdomain number - 1 : subdomain1_outer_boundary_verts, - 2 : subdomain2_outer_boundary_verts -} +# MODEL CONFIGURATION ######################################################### +"""TPR 2 patch soil simulation. -# adjacent_subdomains[i] contains the indices of the subdomains sharing the -# interface i (i.e. given by interface_def_points[i]). -adjacent_subdomains = [[1,2]] +This program sets up an LDD simulation +""" + +import dolfin as df +import sympy as sym +import functools as ft +import LDDsimulation as ldd +import helpers as hlp +import datetime +import os +import pandas as pd +import multiprocessing as mp +import domainSubstructuring as dss + +# init sympy session +sym.init_printing() + +# PREREQUISITS ############################################################### +# check if output directory "./output" exists. This will be used in +# the generation of the output string. +if not os.path.exists('./output'): + os.mkdir('./output') + print("Directory ", './output', " created ") +else: + print("Directory ", './output', " already exists. Will use as output \ + directory") + +date = datetime.datetime.now() +datestr = date.strftime("%Y-%m-%d") + +# Name of the usecase that will be printed during simulation. +use_case = "TPR-2-patch-realistic" +# The name of this very file. Needed for creating log output. +thisfile = "TP-R-2-patch-realistic.py" + +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep +max_iter_num = 500 +FEM_Lagrange_degree = 1 + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### +mesh_study = False +resolutions = { + # 1: 1e-5, + # 2: 1e-5, + # 4: 1e-5, + # 8: 1e-5, + 16: 5e-6, + # 32: 5e-6, + # 64: 2e-6, + # 128: 2e-6, + # 256: 1e-6, + } + +# starttimes gives a list of starttimes to run the simulation from. +# The list is looped over and a simulation is run with t_0 as initial time +# for each element t_0 in starttimes. +starttimes = [0.0] +timestep_size = 0.001 +number_of_timesteps = 5 + +# LDD scheme parameters ###################################################### +Lw1 = 0.025 #/timestep_size +Lnw1= 0.025 + +Lw2 = 0.025 #/timestep_size +Lnw2= 0.025 + +lambda_w = 4 +lambda_nw = 4 + +include_gravity = False +debugflag = False +analyse_condition = False + +# I/O CONFIG ################################################################# +# when number_of_timesteps is high, it might take a long time to write all +# timesteps to disk. Therefore, you can choose to only write data of every +# plot_timestep_every timestep to disk. +plot_timestep_every = 1 +# Decide how many timesteps you want analysed. Analysed means, that +# subsequent errors of the L-iteration within the timestep are written out. +number_of_timesteps_to_analyse = 1 + +# fine grained control over data to be written to disk in the mesh study case +# as well as for a regular simuation for a fixed grid. +if mesh_study: + write_to_file = { + # output the relative errornorm (integration in space) w.r.t. an exact + # solution for each timestep into a csv file. + 'space_errornorms': True, + # save the mesh and marker functions to disk + 'meshes_and_markers': True, + # save xdmf/h5 data for each LDD iteration for timesteps determined by + # number_of_timesteps_to_analyse. I/O intensive! + 'L_iterations_per_timestep': False, + # save solution to xdmf/h5. + 'solutions': True, + # save absolute differences w.r.t an exact solution to xdmf/h5 file + # to monitor where on the domains errors happen + 'absolute_differences': True, + # analyise condition numbers for timesteps determined by + # number_of_timesteps_to_analyse and save them over time to csv. + 'condition_numbers': analyse_condition, + # output subsequent iteration errors measured in L^2 to csv for + # timesteps determined by number_of_timesteps_to_analyse. + # Usefull to monitor convergence of the acutal LDD solver. + 'subsequent_errors': True + } +else: + write_to_file = { + 'space_errornorms': True, + 'meshes_and_markers': True, + 'L_iterations_per_timestep': False, + 'solutions': True, + 'absolute_differences': True, + 'condition_numbers': analyse_condition, + 'subsequent_errors': True + } + +# OUTPUT FILE STRING ######################################################### +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) + +# DOMAIN AND INTERFACE ####################################################### +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points + +# MODEL CONFIGURATION ######################################################### isRichards = { 1: True, # 2: False @@ -380,14 +505,10 @@ p_e_sym = { 'nonwetting': -t**2*(1.1 + y*x)*y**4},#*(sym.sin((x*y-0.5*t)*y**3))**4}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, } -pc_e_sym = dict() -for subdomain, isR in isRichards.items(): - if isR: - pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) - else: - pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - - p_e_sym[subdomain]['wetting'].copy()}) - +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) symbols = {"x": x, "y": y, @@ -413,34 +534,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 @@ -450,88 +555,63 @@ print(f.read()) f.close() -# RUN ######################################################################### -for starttime in starttimes: - for mesh_resolution, solver_tol in resolutions.items(): - # initialise LDD simulation class - simulation = ldd.LDDsimulation( - tol=1E-14, - LDDsolver_tol=solver_tol, - debug=debugflag, - max_iter_num=max_iter_num, - FEM_Lagrange_degree=FEM_Lagrange_degree, - mesh_study=mesh_study - ) - - simulation.set_parameters( - use_case=use_case, - output_dir=output_string, - subdomain_def_points=subdomain_def_points, - isRichards=isRichards, - interface_def_points=interface_def_points, - outer_boundary_def_points=outer_boundary_def_points, - adjacent_subdomains=adjacent_subdomains, - mesh_resolution=mesh_resolution, - viscosity=viscosity, - porosity=porosity, - L=L, - lambda_param=lambda_param, - relative_permeability=relative_permeability, - saturation=sat_pressure_relationship, - starttime=starttime, - number_of_timesteps=number_of_timesteps, - number_of_timesteps_to_analyse=number_of_timesteps_to_analyse, - plot_timestep_every=plot_timestep_every, - timestep_size=timestep_size, - sources=source_expression, - initial_conditions=initial_condition, - dirichletBC_expression_strings=dirichletBC, - exact_solution=exact_solution, - densities=densities, - include_gravity=include_gravity, - gravity_acceleration=gravity_acceleration, - write2file=write_to_file, - ) - - simulation.initialise() - output_dir = simulation.output_dir - # simulation.write_exact_solution_to_xdmf() - output = simulation.run(analyse_condition=analyse_condition) - for subdomain_index, subdomain_output in output.items(): - mesh_h = subdomain_output['mesh_size'] - for phase, error_dict in subdomain_output['errornorm'].items(): - filename = output_dir \ - + "subdomain{}".format(subdomain_index)\ - + "-space-time-errornorm-{}-phase.csv".format(phase) - # for errortype, errornorm in error_dict.items(): - - # eocfile = open("eoc_filename", "a") - # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) - # eocfile.close() - # if subdomain.isRichards:mesh_h - data_dict = { - 'mesh_parameter': mesh_resolution, - 'mesh_h': mesh_h, - } - for norm_type, errornorm in error_dict.items(): - data_dict.update( - {norm_type: errornorm} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) is True: - with open(filename, 'a') as f: - errors.to_csv( - f, - header=False, - sep='\t', - encoding='utf-8', - index=False +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": output_string, + "subdomain_def_points": subdomain_def_points, + "isRichards": isRichards, + "interface_def_points": interface_def_points, + "outer_boundary_def_points": outer_boundary_def_points, + "adjacent_subdomains": adjacent_subdomains, + # "mesh_resolution": mesh_resolution, + "viscosity": viscosity, + "porosity": porosity, + "L": L, + "lambda_param": lambda_param, + "relative_permeability": relative_permeability, + "sat_pressure_relationship": sat_pressure_relationship, + # "starttime": starttime, + "number_of_timesteps": number_of_timesteps, + "number_of_timesteps_to_analyse": number_of_timesteps_to_analyse, + "plot_timestep_every": plot_timestep_every, + "timestep_size": timestep_size, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for starttime in starttimes: + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + starttime, + mesh_resolution ) - else: - errors.to_csv( - filename, - sep='\t', - encoding='utf-8', - index=False ) + LDDsim.start() + LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py index 3f9f35aca65c43aa538f6f0431264a07905c34df..10586745f8a82d2ad21f22b857751715f9b70352 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-gravity-but-same-intrinsic-perm.py @@ -12,6 +12,9 @@ import helpers as hlp import datetime import os import pandas as pd +import multiprocessing as mp +import domainSubstructuring as dss + # init sympy session sym.init_printing() @@ -120,82 +123,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 simulation domain domain -sub_domain0_vertices = [df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0)] -# interface between subdomain1 and subdomain2 -interface12_vertices = [df.Point(-1.0, 0.0), - df.Point(1.0, 0.0)] -# subdomain1. -sub_domain1_vertices = [interface12_vertices[0], - interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3]] - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3], - interface12_vertices[0]] -} -# subdomain2 -sub_domain2_vertices = [sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1], - interface12_vertices[0]] - -subdomain2_outer_boundary_verts = { - 0: [interface12_vertices[0], - sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1]] -} - -# list of subdomains given by the boundary polygon vertices. -# Subdomains are given as a list of dolfin points forming -# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used -# to create the subdomain. subdomain_def_points[0] contains the -# vertices of the global simulation domain and subdomain_def_points[i] contains -# the vertices of the subdomain i. -subdomain_def_points = [ - sub_domain0_vertices, - sub_domain1_vertices, - sub_domain2_vertices - ] -# in the below list, index 0 corresponds to the 12 interface which has index 1 -interface_def_points = [interface12_vertices] - -# if a subdomain has no outer boundary write None instead, i.e. -# i: None -# if i is the index of the inner subdomain. -outer_boundary_def_points = { - # subdomain number - 1: subdomain1_outer_boundary_verts, - 2: subdomain2_outer_boundary_verts -} - -# adjacent_subdomains[i] contains the indices of the subdomains sharing the -# interface i (i.e. given by interface_def_points[i]). -adjacent_subdomains = [[1, 2]] +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points # MODEL CONFIGURATION ######################################################### isRichards = { @@ -445,14 +382,10 @@ p_e_sym = { } #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x) -pc_e_sym = dict() -for subdomain, isR in isRichards.items(): - if isR: - pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) - else: - pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - - p_e_sym[subdomain]['wetting'].copy()}) - +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) symbols = {"x": x, "y": y, @@ -478,35 +411,18 @@ source_expression = exact_solution_example['source'] exact_solution = exact_solution_example['exact_solution'] initial_condition = exact_solution_example['initial_condition'] -# Dictionary of dirichlet boundary conditions. -dirichletBC = dict() -# similarly to the outer boundary dictionary, if a patch has no outer boundary -# None should be written instead of an expression. -# This is a bit of a brainfuck: -# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind. -# Since a domain patch can have several disjoint outer boundary parts, the -# expressions need to get an enumaration index which starts at 0. -# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of -# subdomain ind and boundary part j. -# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting'] -# return the actual expression needed for the dirichlet condition for both -# phases if present. - # BOUNDARY CONDITIONS ######################################################### -# subdomain index: {outer boudary part index: {phase: expression}} -for subdomain in isRichards.keys(): - # subdomain can have no outer boundary - if outer_boundary_def_points[subdomain] is None: - dirichletBC.update({subdomain: None}) - else: - dirichletBC.update({subdomain: dict()}) - # set the dirichlet conditions to be the same code as exact solution on - # the subdomain. - for outer_boundary_ind in outer_boundary_def_points[subdomain].keys(): - dirichletBC[subdomain].update( - {outer_boundary_ind: exact_solution[subdomain]} - ) - +# Dictionary of dirichlet boundary conditions. If an exact solution case is +# used, use the hlp.generate_exact_DirichletBC() method to generate the +# Dirichlet Boundary conditions from the exact solution. +dirichletBC = hlp.generate_exact_DirichletBC( + isRichards=isRichards, + outer_boundary_def_points=outer_boundary_def_points, + exact_solution=exact_solution + ) +# If no exact solution is provided you need to provide a dictionary of boundary +# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see +# the structure. # LOG FILE OUTPUT ############################################################# # read this file and print it to std out. This way the simulation can produce a @@ -516,88 +432,63 @@ print(f.read()) f.close() -# RUN ######################################################################### -for starttime in starttimes: - for mesh_resolution, solver_tol in resolutions.items(): - # initialise LDD simulation class - simulation = ldd.LDDsimulation( - tol=1E-14, - LDDsolver_tol=solver_tol, - debug=debugflag, - max_iter_num=max_iter_num, - FEM_Lagrange_degree=FEM_Lagrange_degree, - mesh_study=mesh_study - ) - - simulation.set_parameters( - use_case=use_case, - output_dir=output_string, - subdomain_def_points=subdomain_def_points, - isRichards=isRichards, - interface_def_points=interface_def_points, - outer_boundary_def_points=outer_boundary_def_points, - adjacent_subdomains=adjacent_subdomains, - mesh_resolution=mesh_resolution, - viscosity=viscosity, - porosity=porosity, - L=L, - lambda_param=lambda_param, - relative_permeability=relative_permeability, - saturation=sat_pressure_relationship, - starttime=starttime, - number_of_timesteps=number_of_timesteps, - number_of_timesteps_to_analyse=number_of_timesteps_to_analyse, - plot_timestep_every=plot_timestep_every, - timestep_size=timestep_size, - sources=source_expression, - initial_conditions=initial_condition, - dirichletBC_expression_strings=dirichletBC, - exact_solution=exact_solution, - densities=densities, - include_gravity=include_gravity, - gravity_acceleration=gravity_acceleration, - write2file=write_to_file, - ) - - simulation.initialise() - output_dir = simulation.output_dir - # simulation.write_exact_solution_to_xdmf() - output = simulation.run(analyse_condition=analyse_condition) - for subdomain_index, subdomain_output in output.items(): - mesh_h = subdomain_output['mesh_size'] - for phase, error_dict in subdomain_output['errornorm'].items(): - filename = output_dir \ - + "subdomain{}".format(subdomain_index)\ - + "-space-time-errornorm-{}-phase.csv".format(phase) - # for errortype, errornorm in error_dict.items(): - - # eocfile = open("eoc_filename", "a") - # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) - # eocfile.close() - # if subdomain.isRichards:mesh_h - data_dict = { - 'mesh_parameter': mesh_resolution, - 'mesh_h': mesh_h, - } - for norm_type, errornorm in error_dict.items(): - data_dict.update( - {norm_type: errornorm} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) is True: - with open(filename, 'a') as f: - errors.to_csv( - f, - header=False, - sep='\t', - encoding='utf-8', - index=False +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": output_string, + "subdomain_def_points": subdomain_def_points, + "isRichards": isRichards, + "interface_def_points": interface_def_points, + "outer_boundary_def_points": outer_boundary_def_points, + "adjacent_subdomains": adjacent_subdomains, + # "mesh_resolution": mesh_resolution, + "viscosity": viscosity, + "porosity": porosity, + "L": L, + "lambda_param": lambda_param, + "relative_permeability": relative_permeability, + "sat_pressure_relationship": sat_pressure_relationship, + # "starttime": starttime, + "number_of_timesteps": number_of_timesteps, + "number_of_timesteps_to_analyse": number_of_timesteps_to_analyse, + "plot_timestep_every": plot_timestep_every, + "timestep_size": timestep_size, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for starttime in starttimes: + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + starttime, + mesh_resolution ) - else: - errors.to_csv( - filename, - sep='\t', - encoding='utf-8', - index=False ) + LDDsim.start() + LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py index d51852eb0b811ee78b3aad39c593a54f15cd1508..f2f48ffab2ff83c1f1a1cbc3c9fd1223302d289e 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-perm.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-no-gravity-but-varying-intrinsic-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,6 +11,8 @@ import helpers as hlp import datetime import os import pandas as pd +import multiprocessing as mp +import domainSubstructuring as dss # init sympy session sym.init_printing() @@ -120,82 +121,16 @@ else: } # OUTPUT FILE STRING ######################################################### -if mesh_study: - output_string = "./output/{}-{}_timesteps{}_P{}".format( - datestr, use_case, number_of_timesteps, FEM_Lagrange_degree - ) -else: - for tol in resolutions.values(): - solver_tol = tol - output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format( - datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol - ) - +output_string = "./output/{}-{}_timesteps{}_P{}".format( + datestr, use_case, number_of_timesteps, FEM_Lagrange_degree + ) # DOMAIN AND INTERFACE ####################################################### -# global simulation domain domain -sub_domain0_vertices = [df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0)] -# interface between subdomain1 and subdomain2 -interface12_vertices = [df.Point(-1.0, 0.0), - df.Point(1.0, 0.0)] -# subdomain1. -sub_domain1_vertices = [interface12_vertices[0], - interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3]] - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3], - interface12_vertices[0]] -} -# subdomain2 -sub_domain2_vertices = [sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1], - interface12_vertices[0]] - -subdomain2_outer_boundary_verts = { - 0: [interface12_vertices[0], - sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1]] -} - -# list of subdomains given by the boundary polygon vertices. -# Subdomains are given as a list of dolfin points forming -# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used -# to create the subdomain. subdomain_def_points[0] contains the -# vertices of the global simulation domain and subdomain_def_points[i] contains -# the vertices of the subdomain i. -subdomain_def_points = [ - sub_domain0_vertices, - sub_domain1_vertices, - sub_domain2_vertices - ] -# in the below list, index 0 corresponds to the 12 interface which has index 1 -interface_def_points = [interface12_vertices] - -# if a subdomain has no outer boundary write None instead, i.e. -# i: None -# if i is the index of the inner subdomain. -outer_boundary_def_points = { - # subdomain number - 1: subdomain1_outer_boundary_verts, - 2: subdomain2_outer_boundary_verts -} - -# adjacent_subdomains[i] contains the indices of the subdomains sharing the -# interface i (i.e. given by interface_def_points[i]). -adjacent_subdomains = [[1, 2]] +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points # MODEL CONFIGURATION ######################################################### isRichards = { @@ -444,15 +379,10 @@ p_e_sym = { 'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2}, } #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x) - -pc_e_sym = dict() -for subdomain, isR in isRichards.items(): - if isR: - pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) - else: - pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - - p_e_sym[subdomain]['wetting'].copy()}) - +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) symbols = {"x": x, "y": y, @@ -478,35 +408,18 @@ source_expression = exact_solution_example['source'] exact_solution = exact_solution_example['exact_solution'] initial_condition = exact_solution_example['initial_condition'] -# Dictionary of dirichlet boundary conditions. -dirichletBC = dict() -# similarly to the outer boundary dictionary, if a patch has no outer boundary -# None should be written instead of an expression. -# This is a bit of a brainfuck: -# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind. -# Since a domain patch can have several disjoint outer boundary parts, the -# expressions need to get an enumaration index which starts at 0. -# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of -# subdomain ind and boundary part j. -# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting'] -# return the actual expression needed for the dirichlet condition for both -# phases if present. - # BOUNDARY CONDITIONS ######################################################### -# subdomain index: {outer boudary part index: {phase: expression}} -for subdomain in isRichards.keys(): - # subdomain can have no outer boundary - if outer_boundary_def_points[subdomain] is None: - dirichletBC.update({subdomain: None}) - else: - dirichletBC.update({subdomain: dict()}) - # set the dirichlet conditions to be the same code as exact solution on - # the subdomain. - for outer_boundary_ind in outer_boundary_def_points[subdomain].keys(): - dirichletBC[subdomain].update( - {outer_boundary_ind: exact_solution[subdomain]} - ) - +# Dictionary of dirichlet boundary conditions. If an exact solution case is +# used, use the hlp.generate_exact_DirichletBC() method to generate the +# Dirichlet Boundary conditions from the exact solution. +dirichletBC = hlp.generate_exact_DirichletBC( + isRichards=isRichards, + outer_boundary_def_points=outer_boundary_def_points, + exact_solution=exact_solution + ) +# If no exact solution is provided you need to provide a dictionary of boundary +# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see +# the structure. # LOG FILE OUTPUT ############################################################# # read this file and print it to std out. This way the simulation can produce a @@ -516,88 +429,63 @@ print(f.read()) f.close() -# RUN ######################################################################### -for starttime in starttimes: - for mesh_resolution, solver_tol in resolutions.items(): - # initialise LDD simulation class - simulation = ldd.LDDsimulation( - tol=1E-14, - LDDsolver_tol=solver_tol, - debug=debugflag, - max_iter_num=max_iter_num, - FEM_Lagrange_degree=FEM_Lagrange_degree, - mesh_study=mesh_study - ) - - simulation.set_parameters( - use_case=use_case, - output_dir=output_string, - subdomain_def_points=subdomain_def_points, - isRichards=isRichards, - interface_def_points=interface_def_points, - outer_boundary_def_points=outer_boundary_def_points, - adjacent_subdomains=adjacent_subdomains, - mesh_resolution=mesh_resolution, - viscosity=viscosity, - porosity=porosity, - L=L, - lambda_param=lambda_param, - relative_permeability=relative_permeability, - saturation=sat_pressure_relationship, - starttime=starttime, - number_of_timesteps=number_of_timesteps, - number_of_timesteps_to_analyse=number_of_timesteps_to_analyse, - plot_timestep_every=plot_timestep_every, - timestep_size=timestep_size, - sources=source_expression, - initial_conditions=initial_condition, - dirichletBC_expression_strings=dirichletBC, - exact_solution=exact_solution, - densities=densities, - include_gravity=include_gravity, - gravity_acceleration=gravity_acceleration, - write2file=write_to_file, - ) - - simulation.initialise() - output_dir = simulation.output_dir - # simulation.write_exact_solution_to_xdmf() - output = simulation.run(analyse_condition=analyse_condition) - for subdomain_index, subdomain_output in output.items(): - mesh_h = subdomain_output['mesh_size'] - for phase, error_dict in subdomain_output['errornorm'].items(): - filename = output_dir \ - + "subdomain{}".format(subdomain_index)\ - + "-space-time-errornorm-{}-phase.csv".format(phase) - # for errortype, errornorm in error_dict.items(): - - # eocfile = open("eoc_filename", "a") - # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) - # eocfile.close() - # if subdomain.isRichards:mesh_h - data_dict = { - 'mesh_parameter': mesh_resolution, - 'mesh_h': mesh_h, - } - for norm_type, errornorm in error_dict.items(): - data_dict.update( - {norm_type: errornorm} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) is True: - with open(filename, 'a') as f: - errors.to_csv( - f, - header=False, - sep='\t', - encoding='utf-8', - index=False +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": output_string, + "subdomain_def_points": subdomain_def_points, + "isRichards": isRichards, + "interface_def_points": interface_def_points, + "outer_boundary_def_points": outer_boundary_def_points, + "adjacent_subdomains": adjacent_subdomains, + # "mesh_resolution": mesh_resolution, + "viscosity": viscosity, + "porosity": porosity, + "L": L, + "lambda_param": lambda_param, + "relative_permeability": relative_permeability, + "sat_pressure_relationship": sat_pressure_relationship, + # "starttime": starttime, + "number_of_timesteps": number_of_timesteps, + "number_of_timesteps_to_analyse": number_of_timesteps_to_analyse, + "plot_timestep_every": plot_timestep_every, + "timestep_size": timestep_size, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for starttime in starttimes: + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + starttime, + mesh_resolution ) - else: - errors.to_csv( - filename, - sep='\t', - encoding='utf-8', - index=False ) + LDDsim.start() + LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py index 163ff3d323e60e006beaeb4f68da10c1ace2c412..92b123ffc013ce529113510ea5c37adcf068c34d 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-perm.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic-same-intrinsic-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,6 +11,8 @@ import helpers as hlp import datetime import os import pandas as pd +import multiprocessing as mp +import domainSubstructuring as dss # init sympy session sym.init_printing() @@ -125,69 +126,11 @@ output_string = "./output/{}-{}_timesteps{}_P{}".format( ) # DOMAIN AND INTERFACE ####################################################### -# global simulation domain domain -sub_domain0_vertices = [df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0)] -# interface between subdomain1 and subdomain2 -interface12_vertices = [df.Point(-1.0, 0.0), - df.Point(1.0, 0.0)] -# subdomain1. -sub_domain1_vertices = [interface12_vertices[0], - interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3]] - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for -# defining the Dirichlet boundary conditions. If a domain is completely -# internal, the dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3], - interface12_vertices[0]] -} -# subdomain2 -sub_domain2_vertices = [sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1], - interface12_vertices[0]] - -subdomain2_outer_boundary_verts = { - 0: [interface12_vertices[0], - sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1]] -} - -# list of subdomains given by the boundary polygon vertices. -# Subdomains are given as a list of dolfin points forming -# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used -# to create the subdomain. subdomain_def_points[0] contains the -# vertices of the global simulation domain and subdomain_def_points[i] contains -# the vertices of the subdomain i. -subdomain_def_points = [ - sub_domain0_vertices, - sub_domain1_vertices, - sub_domain2_vertices - ] -# in the below list, index 0 corresponds to the 12 interface which has index 1 -interface_def_points = [interface12_vertices] - -# if a subdomain has no outer boundary write None instead, i.e. -# i: None -# if i is the index of the inner subdomain. -outer_boundary_def_points = { - # subdomain number - 1: subdomain1_outer_boundary_verts, - 2: subdomain2_outer_boundary_verts -} - -# adjacent_subdomains[i] contains the indices of the subdomains sharing the -# interface i (i.e. given by interface_def_points[i]). -adjacent_subdomains = [[1, 2]] +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points # MODEL CONFIGURATION ######################################################### isRichards = { @@ -436,15 +379,10 @@ p_e_sym = { 'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2}, } #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x) - -pc_e_sym = dict() -for subdomain, isR in isRichards.items(): - if isR: - pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) - else: - pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - - p_e_sym[subdomain]['wetting'].copy()}) - +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) symbols = {"x": x, "y": y, @@ -470,35 +408,18 @@ source_expression = exact_solution_example['source'] exact_solution = exact_solution_example['exact_solution'] initial_condition = exact_solution_example['initial_condition'] -# Dictionary of dirichlet boundary conditions. -dirichletBC = dict() -# similarly to the outer boundary dictionary, if a patch has no outer boundary -# None should be written instead of an expression. -# This is a bit of a brainfuck: -# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind. -# Since a domain patch can have several disjoint outer boundary parts, the -# expressions need to get an enumaration index which starts at 0. -# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of -# subdomain ind and boundary part j. -# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting'] -# return the actual expression needed for the dirichlet condition for both -# phases if present. - # BOUNDARY CONDITIONS ######################################################### -# subdomain index: {outer boudary part index: {phase: expression}} -for subdomain in isRichards.keys(): - # subdomain can have no outer boundary - if outer_boundary_def_points[subdomain] is None: - dirichletBC.update({subdomain: None}) - else: - dirichletBC.update({subdomain: dict()}) - # set the dirichlet conditions to be the same code as exact solution on - # the subdomain. - for outer_boundary_ind in outer_boundary_def_points[subdomain].keys(): - dirichletBC[subdomain].update( - {outer_boundary_ind: exact_solution[subdomain]} - ) - +# Dictionary of dirichlet boundary conditions. If an exact solution case is +# used, use the hlp.generate_exact_DirichletBC() method to generate the +# Dirichlet Boundary conditions from the exact solution. +dirichletBC = hlp.generate_exact_DirichletBC( + isRichards=isRichards, + outer_boundary_def_points=outer_boundary_def_points, + exact_solution=exact_solution + ) +# If no exact solution is provided you need to provide a dictionary of boundary +# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see +# the structure. # LOG FILE OUTPUT ############################################################# # read this file and print it to std out. This way the simulation can produce a @@ -508,88 +429,63 @@ print(f.read()) f.close() -# RUN ######################################################################### -for starttime in starttimes: - for mesh_resolution, solver_tol in resolutions.items(): - # initialise LDD simulation class - simulation = ldd.LDDsimulation( - tol=1E-14, - LDDsolver_tol=solver_tol, - debug=debugflag, - max_iter_num=max_iter_num, - FEM_Lagrange_degree=FEM_Lagrange_degree, - mesh_study=mesh_study - ) - - simulation.set_parameters( - use_case=use_case, - output_dir=output_string, - subdomain_def_points=subdomain_def_points, - isRichards=isRichards, - interface_def_points=interface_def_points, - outer_boundary_def_points=outer_boundary_def_points, - adjacent_subdomains=adjacent_subdomains, - mesh_resolution=mesh_resolution, - viscosity=viscosity, - porosity=porosity, - L=L, - lambda_param=lambda_param, - relative_permeability=relative_permeability, - saturation=sat_pressure_relationship, - starttime=starttime, - number_of_timesteps=number_of_timesteps, - number_of_timesteps_to_analyse=number_of_timesteps_to_analyse, - plot_timestep_every=plot_timestep_every, - timestep_size=timestep_size, - sources=source_expression, - initial_conditions=initial_condition, - dirichletBC_expression_strings=dirichletBC, - exact_solution=exact_solution, - densities=densities, - include_gravity=include_gravity, - gravity_acceleration=gravity_acceleration, - write2file=write_to_file, - ) - - simulation.initialise() - output_dir = simulation.output_dir - # simulation.write_exact_solution_to_xdmf() - output = simulation.run(analyse_condition=analyse_condition) - for subdomain_index, subdomain_output in output.items(): - mesh_h = subdomain_output['mesh_size'] - for phase, error_dict in subdomain_output['errornorm'].items(): - filename = output_dir \ - + "subdomain{}".format(subdomain_index)\ - + "-space-time-errornorm-{}-phase.csv".format(phase) - # for errortype, errornorm in error_dict.items(): - - # eocfile = open("eoc_filename", "a") - # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) - # eocfile.close() - # if subdomain.isRichards:mesh_h - data_dict = { - 'mesh_parameter': mesh_resolution, - 'mesh_h': mesh_h, - } - for norm_type, errornorm in error_dict.items(): - data_dict.update( - {norm_type: errornorm} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) is True: - with open(filename, 'a') as f: - errors.to_csv( - f, - header=False, - sep='\t', - encoding='utf-8', - index=False +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": output_string, + "subdomain_def_points": subdomain_def_points, + "isRichards": isRichards, + "interface_def_points": interface_def_points, + "outer_boundary_def_points": outer_boundary_def_points, + "adjacent_subdomains": adjacent_subdomains, + # "mesh_resolution": mesh_resolution, + "viscosity": viscosity, + "porosity": porosity, + "L": L, + "lambda_param": lambda_param, + "relative_permeability": relative_permeability, + "sat_pressure_relationship": sat_pressure_relationship, + # "starttime": starttime, + "number_of_timesteps": number_of_timesteps, + "number_of_timesteps_to_analyse": number_of_timesteps_to_analyse, + "plot_timestep_every": plot_timestep_every, + "timestep_size": timestep_size, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for starttime in starttimes: + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + starttime, + mesh_resolution ) - else: - errors.to_csv( - filename, - sep='\t', - encoding='utf-8', - index=False ) + LDDsim.start() + LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py index 5cad90c4e21b5b79b13080d2a58663cb14736af8..c95a5d15f6709b4617b40d15fd8f8e6bd2c287ef 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-realistic.py @@ -122,17 +122,9 @@ 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 ####################################################### substructuring = dss.twoSoilLayers() diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py index a3e6306eaa8d6fed2c84d2344d74ea4074fe7fb5..09106bdcd97cc86a865d580128ff3ae2bdc0ffd1 100755 --- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.py +++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-test.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() @@ -124,69 +125,12 @@ output_string = "./output/{}-{}_timesteps{}_P{}".format( datestr, use_case, number_of_timesteps, FEM_Lagrange_degree ) - # DOMAIN AND INTERFACE ####################################################### -# global simulation domain domain -sub_domain0_vertices = [df.Point(-1.0, -1.0), - df.Point(1.0, -1.0), - df.Point(1.0, 1.0), - df.Point(-1.0, 1.0)] -# interface between subdomain1 and subdomain2 -interface12_vertices = [df.Point(-1.0, 0.0), - df.Point(1.0, 0.0) ] -# subdomain1. -sub_domain1_vertices = [interface12_vertices[0], - interface12_vertices[1], - sub_domain0_vertices[2], - sub_domain0_vertices[3]] - -# vertex coordinates of the outer boundaries. If it can not be specified as a -# polygon, use an entry per boundary polygon. This information is used for defining -# the Dirichlet boundary conditions. If a domain is completely internal, the -# dictionary entry should be 0: None -subdomain1_outer_boundary_verts = { - 0: [interface12_vertices[1], # - sub_domain0_vertices[2], - sub_domain0_vertices[3], # - interface12_vertices[0]] -} -# subdomain2 -sub_domain2_vertices = [sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1], - interface12_vertices[0] ] - -subdomain2_outer_boundary_verts = { - 0: [interface12_vertices[0], # - sub_domain0_vertices[0], - sub_domain0_vertices[1], - interface12_vertices[1]] -} - -# list of subdomains given by the boundary polygon vertices. -# Subdomains are given as a list of dolfin points forming -# a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used -# to create the subdomain. subdomain_def_points[0] contains the -# vertices of the global simulation domain and subdomain_def_points[i] contains the -# vertices of the subdomain i. -subdomain_def_points = [sub_domain0_vertices,# - sub_domain1_vertices,# - sub_domain2_vertices] -# in the below list, index 0 corresponds to the 12 interface which has index 1 -interface_def_points = [interface12_vertices] - -# if a subdomain has no outer boundary write None instead, i.e. -# i: None -# if i is the index of the inner subdomain. -outer_boundary_def_points = { - # subdomain number - 1 : subdomain1_outer_boundary_verts, - 2 : subdomain2_outer_boundary_verts -} - -# adjacent_subdomains[i] contains the indices of the subdomains sharing the -# interface i (i.e. given by interface_def_points[i]). -adjacent_subdomains = [[1,2]] +substructuring = dss.twoSoilLayers() +interface_def_points = substructuring.interface_def_points +adjacent_subdomains = substructuring.adjacent_subdomains +subdomain_def_points = substructuring.subdomain_def_points +outer_boundary_def_points = substructuring.outer_boundary_def_points # MODEL CONFIGURATION ######################################################### isRichards = { @@ -423,15 +367,10 @@ p_e_sym = { 'nonwetting': (-2-t*(1.1+y + x**2))*y**2}, #*(1-x)**2*(1+x)**2*(1+y)**2}, } #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x) - -pc_e_sym = dict() -for subdomain, isR in isRichards.items(): - if isR: - pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting'].copy()}) - else: - pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting'].copy() - - p_e_sym[subdomain]['wetting'].copy()}) - +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) symbols = {"x": x, "y": y, @@ -458,34 +397,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(): - # subdomain can have no outer boundary - if outer_boundary_def_points[subdomain] is None: - dirichletBC.update({subdomain: None}) - else: - dirichletBC.update({subdomain: dict()}) - # set the dirichlet conditions to be the same code as exact solution on - # the subdomain. - for outer_boundary_ind in outer_boundary_def_points[subdomain].keys(): - dirichletBC[subdomain].update( - {outer_boundary_ind: exact_solution[subdomain]} - ) - +# Dictionary of dirichlet boundary conditions. If an exact solution case is +# used, use the hlp.generate_exact_DirichletBC() method to generate the +# Dirichlet Boundary conditions from the exact solution. +dirichletBC = hlp.generate_exact_DirichletBC( + isRichards=isRichards, + outer_boundary_def_points=outer_boundary_def_points, + exact_solution=exact_solution + ) +# If no exact solution is provided you need to provide a dictionary of boundary +# conditions. See the definiton of hlp.generate_exact_DirichletBC() to see +# the structure. # LOG FILE OUTPUT ############################################################# # read this file and print it to std out. This way the simulation can produce a @@ -495,88 +417,63 @@ print(f.read()) f.close() -# RUN ######################################################################### -for starttime in starttimes: - for mesh_resolution, solver_tol in resolutions.items(): - # initialise LDD simulation class - simulation = ldd.LDDsimulation( - tol=1E-14, - LDDsolver_tol=solver_tol, - debug=debugflag, - max_iter_num=max_iter_num, - FEM_Lagrange_degree=FEM_Lagrange_degree, - mesh_study=mesh_study - ) - - simulation.set_parameters( - use_case=use_case, - output_dir=output_string, - subdomain_def_points=subdomain_def_points, - isRichards=isRichards, - interface_def_points=interface_def_points, - outer_boundary_def_points=outer_boundary_def_points, - adjacent_subdomains=adjacent_subdomains, - mesh_resolution=mesh_resolution, - viscosity=viscosity, - porosity=porosity, - L=L, - lambda_param=lambda_param, - relative_permeability=relative_permeability, - saturation=sat_pressure_relationship, - starttime=starttime, - number_of_timesteps=number_of_timesteps, - number_of_timesteps_to_analyse=number_of_timesteps_to_analyse, - plot_timestep_every=plot_timestep_every, - timestep_size=timestep_size, - sources=source_expression, - initial_conditions=initial_condition, - dirichletBC_expression_strings=dirichletBC, - exact_solution=exact_solution, - densities=densities, - include_gravity=include_gravity, - gravity_acceleration=gravity_acceleration, - write2file=write_to_file, - ) - - simulation.initialise() - output_dir = simulation.output_dir - # simulation.write_exact_solution_to_xdmf() - output = simulation.run(analyse_condition=analyse_condition) - for subdomain_index, subdomain_output in output.items(): - mesh_h = subdomain_output['mesh_size'] - for phase, error_dict in subdomain_output['errornorm'].items(): - filename = output_dir \ - + "subdomain{}".format(subdomain_index)\ - + "-space-time-errornorm-{}-phase.csv".format(phase) - # for errortype, errornorm in error_dict.items(): - - # eocfile = open("eoc_filename", "a") - # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" ) - # eocfile.close() - # if subdomain.isRichards:mesh_h - data_dict = { - 'mesh_parameter': mesh_resolution, - 'mesh_h': mesh_h, - } - for norm_type, errornorm in error_dict.items(): - data_dict.update( - {norm_type: errornorm} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) is True: - with open(filename, 'a') as f: - errors.to_csv( - f, - header=False, - sep='\t', - encoding='utf-8', - index=False +# MAIN ######################################################################## +if __name__ == '__main__': + # dictionary of simualation parameters to pass to the run function. + # mesh_resolution and starttime are excluded, as they get passed explicitly + # to achieve parallelisation in these parameters in these parameters for + # mesh studies etc. + simulation_parameter = { + "tol": 1E-14, + "debugflag": debugflag, + "max_iter_num": max_iter_num, + "FEM_Lagrange_degree": FEM_Lagrange_degree, + "mesh_study": mesh_study, + "use_case": use_case, + "output_string": output_string, + "subdomain_def_points": subdomain_def_points, + "isRichards": isRichards, + "interface_def_points": interface_def_points, + "outer_boundary_def_points": outer_boundary_def_points, + "adjacent_subdomains": adjacent_subdomains, + # "mesh_resolution": mesh_resolution, + "viscosity": viscosity, + "porosity": porosity, + "L": L, + "lambda_param": lambda_param, + "relative_permeability": relative_permeability, + "sat_pressure_relationship": sat_pressure_relationship, + # "starttime": starttime, + "number_of_timesteps": number_of_timesteps, + "number_of_timesteps_to_analyse": number_of_timesteps_to_analyse, + "plot_timestep_every": plot_timestep_every, + "timestep_size": timestep_size, + "source_expression": source_expression, + "initial_condition": initial_condition, + "dirichletBC": dirichletBC, + "exact_solution": exact_solution, + "densities": densities, + "include_gravity": include_gravity, + "gravity_acceleration": gravity_acceleration, + "write_to_file": write_to_file, + "analyse_condition": analyse_condition + } + for starttime in starttimes: + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + starttime, + mesh_resolution ) - else: - errors.to_csv( - filename, - sep='\t', - encoding='utf-8', - index=False ) + LDDsim.start() + LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # )