diff --git a/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-different-intrinsic-perm.py b/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-different-intrinsic-perm.py index 6e211c020c0b968bcec38331c3ac1a32f8224c77..13c575d944a84911b3b65bda8bb3c3ae546b3c47 100755 --- a/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-different-intrinsic-perm.py +++ b/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-different-intrinsic-perm.py @@ -30,7 +30,7 @@ date = datetime.datetime.now() datestr = date.strftime("%Y-%m-%d") # Name of the usecase that will be printed during simulation. -use_case = "TP-TP-2P-realistic-different-intrinsic-perm-smaller-L" +use_case = "TP-TP-2P-realistic-different-intrinsic-perm-params1a-higher-tol" # The name of this very file. Needed for creating log output. thisfile = "TP-TP-2-patch-different-intrinsic-perm.py" @@ -61,14 +61,14 @@ timestep_size = 0.001 number_of_timesteps = 1000 # LDD scheme parameters ###################################################### -Lw1 = 0.01 #/timestep_size -Lnw1= 0.01 +Lw1 = 0.1 #/timestep_size +Lnw1= 0.1 -Lw2 = 0.001 #/timestep_size -Lnw2= 0.001 +Lw2 = 0.1 #/timestep_size +Lnw2= 0.1 -lambda_w = 2 -lambda_nw = 2 +lambda_w = 40 +lambda_nw = 40 include_gravity = False debugflag = False analyse_condition = False diff --git a/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-same-intrinsic-perm.py b/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-same-intrinsic-perm.py index e1069be0dbb6991d03211d7423491def69011792..75e9b5242ac94883198bcaaadb7a164b23a484d1 100755 --- a/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-same-intrinsic-perm.py +++ b/Two-phase-Two-phase/two-patch/TP-TP-2-patch-test-case/TP-TP-2-patch-same-intrinsic-perm.py @@ -29,7 +29,7 @@ date = datetime.datetime.now() datestr = date.strftime("%Y-%m-%d") # Name of the usecase that will be printed during simulation. -use_case = "TP-TP-2P-realistic-same-intrinsic0.1-control-for-case160620" +use_case = "TP-TP-2P-realistic-same-intrinsic" # The name of this very file. Needed for creating log output. thisfile = "TP-TP-2-patch-same-intrinsic-perm.py" @@ -46,7 +46,7 @@ resolutions = { # 4: 1e-6, # 8: 1e-6, # 16: 5e-6, - 32: 3e-6, + 32: 1e-5, # 64: 2e-6, # 128: 1e-6, # 256: 1e-6, @@ -56,20 +56,21 @@ resolutions = { # 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.0} +# starttimes = {0: 0.0, 1:0.3, 2:0.6, 3:0.9} timestep_size = 0.001 number_of_timesteps = 1000 # LDD scheme parameters ###################################################### -Lw1 = 0.25 #/timestep_size -Lnw1= 0.25 +Lw1 = 0.5 #/timestep_size +Lnw1= 0.2 -Lw2 = 0.25 #/timestep_size -Lnw2= 0.25 +Lw2 = 0.5 #/timestep_size +Lnw2= 0.2 -lambda_w = 4 -lambda_nw = 4 +lambda_w = 2 +lambda_nw = 1 -include_gravity = False +include_gravity = True debugflag = False analyse_condition = False @@ -80,7 +81,7 @@ analyse_condition = False 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 +number_of_timesteps_to_analyse = 10 # 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. @@ -176,8 +177,8 @@ lambda_param = {# } intrinsic_permeability = { - 1: 0.1, - 2: 0.1, + 1: 0.01, + 2: 0.01, } # RELATIVE PEMRMEABILITIES @@ -339,7 +340,7 @@ if __name__ == '__main__': # parameter=simulation_parameter # ) - LDDsim.join() + # LDDsim.join() if mesh_study: simulation_output_dir = processQueue.get() hlp.merge_spacetime_errornorms(isRichards=isRichards, diff --git a/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py b/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py index bd6ee937074f6c998921da8e032fea90f675c629..19595959e172816c5835b019246a6a2c604867fe 100755 --- a/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py +++ b/Two-phase-Two-phase/two-patch/mesh_studies/TP-TP-2-patch-mesh-study.py @@ -1,78 +1,109 @@ #!/usr/bin/python3 +"""TP-TP 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 functions as fts import LDDsimulation as ldd -import functools as ft 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 = "TP-TP-2-patch-realistic" -# solver_tol = 6E-7 -max_iter_num = 300 +# The name of this very file. Needed for creating log output. +thisfile = "TP-TP-2-patch-mesh-study.py" +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep +max_iter_num = 400 FEM_Lagrange_degree = 1 -mesh_study = False + +# GRID AND MESH STUDY SPECIFICATIONS ######################################### +mesh_study = True resolutions = { - # 1: 1e-6, - # 2: 1e-6, - # 4: 1e-6, - # 8: 1e-6, - # 16: 1e-6, - 32: 1e-6, - # 64: 1e-6, - # 128: 1e-6, + 1: 1e-5, + 2: 1e-5, + 4: 1e-5, + 8: 1e-5, + 16: 5e-6, + 32: 3e-6, + 64: 1e-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.0} timestep_size = 0.001 -number_of_timesteps = 20 -plot_timestep_every = 2 -# decide how many timesteps you want analysed. Analysed means, that we write out -# subsequent errors of the L-iteration within the timestep. -number_of_timesteps_to_analyse = 6 -starttimes = [0.0] +number_of_timesteps = 1000 -Lw = 0.025 #/timestep_size -Lnw=Lw +# LDD scheme parameters ###################################################### +Lw1 = 0.01 #/timestep_size +Lnw1= 0.01 -lambda_w = 4 -lambda_nw = 4 +Lw2 = 0.001 #/timestep_size +Lnw2= 0.001 -include_gravity = True -debugflag = True -analyse_condition = True +lambda_w = 2 +lambda_nw = 2 -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) +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 = 2 +# 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 = 8 -# 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, - 'L_iterations_per_timestep': 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: @@ -86,70 +117,19 @@ 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 ####################################################### +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 -##### 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]] +# MODEL CONFIGURATION ######################################################### isRichards = { 1: False, # 2: False @@ -158,16 +138,16 @@ isRichards = { viscosity = {# # subdom_num : viscosity - 1 : {'wetting' :1, + 1: {'wetting' :1, 'nonwetting': 1/50}, # - 2 : {'wetting' :1, + 2: {'wetting' :1, 'nonwetting': 1/50} } porosity = {# # subdom_num : porosity - 1 : 0.22,# - 2 : 0.0022 + 1: 0.22,# + 2: 0.022 } # Dict of the form: { subdom_num : density } @@ -175,206 +155,53 @@ densities = { 1: {'wetting': 997, 'nonwetting': 1.225}, 2: {'wetting': 997, - 'nonwetting': 1.225}, -} - -intrinsic_permeability = { - 1: {"wetting": 1, - "nonwetting": 1}, - 2: {"wetting": 1, - "nonwetting": 1}, + 'nonwetting': 1.225} } gravity_acceleration = 9.81 L = {# # subdom_num : subdomain L for L-scheme - 1 : {'wetting' :Lw, - 'nonwetting': Lnw},# - 2 : {'wetting' :Lw, - 'nonwetting': Lnw} + 1 : {'wetting' :Lw1, + 'nonwetting': Lnw1},# + 2 : {'wetting' :Lw2, + 'nonwetting': Lnw2} } lambda_param = {# # subdom_num : lambda parameter for the L-scheme - 1 : {'wetting' :lambda_w, + 0 : {'wetting' :lambda_w, 'nonwetting': lambda_nw},# - 2 : {'wetting' :lambda_w, - 'nonwetting': lambda_nw} -} - -## relative permeabilty functions on subdomain 1 -def rel_perm1w(s): - # relative permeabilty wetting on subdomain1 - return intrinsic_permeability[1]["wetting"]*s**2 - -def rel_perm1nw(s): - # relative permeabilty nonwetting on subdomain1 - return intrinsic_permeability[1]["nonwetting"]*(1-s)**2 - -_rel_perm1w = ft.partial(rel_perm1w) -_rel_perm1nw = ft.partial(rel_perm1nw) -subdomain1_rel_perm = { - 'wetting': _rel_perm1w,# - 'nonwetting': _rel_perm1nw -} -## relative permeabilty functions on subdomain 2 -def rel_perm2w(s): - # relative permeabilty wetting on subdomain2 - return intrinsic_permeability[2]["wetting"]*s**3 - -def rel_perm2nw(s): - # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2 - return intrinsic_permeability[2]["nonwetting"]*(1-s)**3 - -_rel_perm2w = ft.partial(rel_perm2w) -_rel_perm2nw = ft.partial(rel_perm2nw) - -subdomain2_rel_perm = { - 'wetting': _rel_perm2w,# - 'nonwetting': _rel_perm2nw -} - -## dictionary of relative permeabilties on all domains. -relative_permeability = {# - 1: subdomain1_rel_perm, - 2: subdomain2_rel_perm -} - - -# definition of the derivatives of the relative permeabilities -# relative permeabilty functions on subdomain 1 -def rel_perm1w_prime(s): - # relative permeabilty on subdomain1 - return intrinsic_permeability[1]["wetting"]*2*s - -def rel_perm1nw_prime(s): - # relative permeabilty on subdomain1 - return -intrinsic_permeability[1]["nonwetting"]*2*(1-s) - -# definition of the derivatives of the relative permeabilities -# relative permeabilty functions on subdomain 1 -def rel_perm2w_prime(s): - # relative permeabilty on subdomain1 - return intrinsic_permeability[2]["wetting"]*3*s**2 - -def rel_perm2nw_prime(s): - # relative permeabilty on subdomain1 - return -intrinsic_permeability[2]["nonwetting"]*3*(1-s)**2 - -_rel_perm1w_prime = ft.partial(rel_perm1w_prime) -_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime) -_rel_perm2w_prime = ft.partial(rel_perm2w_prime) -_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime) - -subdomain1_rel_perm_prime = { - 'wetting': _rel_perm1w_prime, - 'nonwetting': _rel_perm1nw_prime -} - - -subdomain2_rel_perm_prime = { - 'wetting': _rel_perm2w_prime, - 'nonwetting': _rel_perm2nw_prime } -# dictionary of relative permeabilties on all domains. -ka_prime = { - 1: subdomain1_rel_perm_prime, - 2: subdomain2_rel_perm_prime, +intrinsic_permeability = { + 1: 0.01, + 2: 0.001, } - -# def saturation1(pc, subdomain_index): -# # inverse capillary pressure-saturation-relationship -# return df.conditional(pc > 0, 1/((1 + pc)**(1/(subdomain_index + 1))), 1) -# -# def saturation2(pc, n_index, alpha): -# # inverse capillary pressure-saturation-relationship -# return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1) -# -# # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where -# # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw -# def saturation1_sym(pc, subdomain_index): -# # inverse capillary pressure-saturation-relationship -# return 1/((1 + pc)**(1/(subdomain_index + 1))) -# -# -# def saturation2_sym(pc, n_index, alpha): -# # inverse capillary pressure-saturation-relationship -# #df.conditional(pc > 0, -# return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)) -# -# -# # derivative of S-pc relationship with respect to pc. This is needed for the -# # construction of a analytic solution. -# def saturation1_sym_prime(pc, subdomain_index): -# # inverse capillary pressure-saturation-relationship -# return -(1/(subdomain_index + 1))*(1 + pc)**((-subdomain_index - 2)/(subdomain_index + 1)) -# -# -# def saturation2_sym_prime(pc, n_index, alpha): -# # inverse capillary pressure-saturation-relationship -# return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) ) -# -# # note that the conditional definition of S-pc in the nonsymbolic part will be -# # incorporated in the construction of the exact solution below. -# S_pc_sym = { -# 1: ft.partial(saturation1_sym, subdomain_index = 1), -# 2: ft.partial(saturation2_sym, n_index=3, alpha=0.001), -# } -# -# S_pc_sym_prime = { -# 1: ft.partial(saturation1_sym_prime, subdomain_index = 1), -# 2: ft.partial(saturation2_sym_prime, n_index=3, alpha=0.001), -# } -# -# sat_pressure_relationship = { -# 1: ft.partial(saturation1, subdomain_index = 1),#, -# 2: ft.partial(saturation2, n_index=3, alpha=0.001), -# } - -def saturation(pc, index): - # inverse capillary pressure-saturation-relationship - return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1) - - -def saturation_sym(pc, index): - # inverse capillary pressure-saturation-relationship - return 1/((1 + pc)**(1/(index + 1))) - - -# derivative of S-pc relationship with respect to pc. This is needed for the -# construction of a analytic solution. -def saturation_sym_prime(pc, index): - # inverse capillary pressure-saturation-relationship - return -1/((index+1)*(1 + pc)**((index+2)/(index+1))) - - -# note that the conditional definition of S-pc in the nonsymbolic part will be -# incorporated in the construction of the exact solution below. -S_pc_sym = { - 1: ft.partial(saturation_sym, index=1), - 2: ft.partial(saturation_sym, index=2), - # 3: ft.partial(saturation_sym, index=2), - # 4: ft.partial(saturation_sym, index=1) +# RELATIVE PEMRMEABILITIES +rel_perm_definition = { + 1: {"wetting": "Spow2", + "nonwetting": "oneMinusSpow2"}, + 2: {"wetting": "Spow3", + "nonwetting": "oneMinusSpow3"}, } -S_pc_sym_prime = { - 1: ft.partial(saturation_sym_prime, index=1), - 2: ft.partial(saturation_sym_prime, index=2), - # 3: ft.partial(saturation_sym_prime, index=2), - # 4: ft.partial(saturation_sym_prime, index=1) -} +rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition) +relative_permeability = rel_perm_dict["ka"] +ka_prime = rel_perm_dict["ka_prime"] -sat_pressure_relationship = { - 1: ft.partial(saturation, index=1), - 2: ft.partial(saturation, index=2), - # 3: ft.partial(saturation, index=2), - # 4: ft.partial(saturation, index=1) +# S-pc relation +Spc_on_subdomains = { + 1: {"testSpc": {"index": 1}}, + 2: {"testSpc": {"index": 2}}, } +Spc = fts.generate_Spc_dicts(Spc_on_subdomains) +S_pc_sym = Spc["symbolic"] +S_pc_sym_prime = Spc["prime_symbolic"] +sat_pressure_relationship = Spc["dolfin"] ############################################# # Manufacture source expressions with sympy # @@ -384,20 +211,15 @@ t = sym.symbols('t', positive=True) p_e_sym = { 1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)), - 'nonwetting': -1 -t*(1.1 + y*y)}, #*(1-x)**2*(1+x)**2*(1-y)**2}, - 2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2, - 'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**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()}) + 'nonwetting': (-1 -t*(1.1+ y*y))}, + 2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)), + 'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)}, +} +pc_e_sym = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) symbols = {"x": x, "y": y, @@ -413,6 +235,7 @@ exact_solution_example = hlp.generate_exact_solution_expressions( saturation_pressure_relationship_prime=S_pc_sym_prime, viscosity=viscosity, porosity=porosity, + intrinsic_permeability=intrinsic_permeability, relative_permeability=relative_permeability, relative_permeability_prime=ka_prime, densities=densities, @@ -423,111 +246,96 @@ 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]} - ) - - -# def saturation(pressure, subdomain_index): -# # inverse capillary pressure-saturation-relationship -# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) -# -# sa - -f = open('TP-TP-2-patch-mesh-study.py', 'r') +# 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 +# log file with ./TP-R-layered_soil.py | tee simulation.log +f = open(thisfile, 'r') print(f.read()) f.close() -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, - 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, different_errornorms in subdomain_output['errornorm'].items(): - filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase) - # for errortype, errornorm in different_errornorms.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 error_type, errornorms in different_errornorms.items(): - data_dict.update( - {error_type: errornorms} - ) - errors = pd.DataFrame(data_dict, index=[mesh_resolution]) - # check if file exists - if os.path.isfile(filename) == True: - with open(filename, 'a') as f: - errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False) - else: - errors.to_csv(filename, sep='\t', encoding='utf-8', index=False) +# 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, + "intrinsic_permeability": intrinsic_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 number_shift, starttime in starttimes.items(): + simulation_parameter.update( + {"starttime_timestep_number_shift": number_shift} + ) + for mesh_resolution, solver_tol in resolutions.items(): + simulation_parameter.update({"solver_tol": solver_tol}) + hlp.info(simulation_parameter["use_case"]) + processQueue = mp.Queue() + LDDsim = mp.Process( + target=hlp.run_simulation, + args=( + simulation_parameter, + processQueue, + starttime, + mesh_resolution + ) + ) + LDDsim.start() + # LDDsim.join() + # hlp.run_simulation( + # mesh_resolution=mesh_resolution, + # starttime=starttime, + # parameter=simulation_parameter + # ) + + LDDsim.join() + if mesh_study: + simulation_output_dir = processQueue.get() + hlp.merge_spacetime_errornorms(isRichards=isRichards, + resolutions=resolutions, + output_dir=simulation_output_dir)