diff --git a/Usecases/InstallationTest/installTest.py b/Usecases/InstallationTest/installTest.py new file mode 100755 index 0000000000000000000000000000000000000000..450f6da0d7c07f0a14784a869677103e0bfd9af8 --- /dev/null +++ b/Usecases/InstallationTest/installTest.py @@ -0,0 +1,360 @@ +#!/usr/bin/python3 +"""TPR 2 patch soil simulation, Copyright 2020, David Seus + +This program runs an LDD simulation on a two-domain substructuring using +a TPR coupling. + +# LICENCE ##################################################################### +Copyright 2020, David Seus +david.seus[at]ians.uni-stuttgart.de +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see <http://www.gnu.org/licenses/>. +############################################################################### +""" +import dolfin as df +import sympy as sym +import functions as fts +import LDDsimulation as ldd +import helpers as hlp +import datetime +import os +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 = "InstalationTest_TPR-2-patch" +# The name of this very file. Needed for creating log output. +thisfile = "installTest.py" + +# GENERAL SOLVER CONFIG ###################################################### +# maximal iteration per timestep +max_iter_num = 5 +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.0} +timestep_size = 0.01 +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 + } + + +viscosity = {# +# subdom_num : viscosity + 1: {'wetting' :1, + 'nonwetting': 1/50}, # + 2: {'wetting' :1, + 'nonwetting': 1/50} +} + +porosity = {# +# subdom_num : porosity + 1: 0.22,# + 2: 0.22 +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting': 997, + 'nonwetting': 1.225}, + 2: {'wetting': 997, + 'nonwetting': 1.225} +} + +gravity_acceleration = 9.81 + +L = {# +# subdom_num : subdomain L for L-scheme + 1 : {'wetting' :Lw1, + 'nonwetting': Lnw1},# + 2 : {'wetting' :Lw2, + 'nonwetting': Lnw2} +} + + +lambda_param = {# +# subdom_num : lambda parameter for the L-scheme + 0 : {'wetting' :lambda_w, + 'nonwetting': lambda_nw},# +} + +intrinsic_permeability = { + 1: 0.01, + 2: 0.01, +} + +# RELATIVE PEMRMEABILITIES +rel_perm_definition = { + 1: {"wetting": "Spow2", + "nonwetting": "oneMinusSpow2"}, + 2: {"wetting": "Spow3", + "nonwetting": "oneMinusSpow3"}, +} + +rel_perm_dict = fts.generate_relative_permeability_dicts(rel_perm_definition) +relative_permeability = rel_perm_dict["ka"] +ka_prime = rel_perm_dict["ka_prime"] + +# 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 # +############################################################################### +x, y = sym.symbols('x[0], x[1]') # needed by UFL +t = sym.symbols('t', positive=True) + +p_e_sym = { + 1: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x + y*y))}, #*(1-x)**2*(1+x)**2*(1-y)**2}, + 2: {'wetting': (-7.0 - (1.0 + t*t)*(1.0 + x*x)), #*(1-x)**2*(1+x)**2*(1+y)**2, + '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 = hlp.generate_exact_symbolic_pc( + isRichards=isRichards, + symbolic_pressure=p_e_sym + ) + +symbols = {"x": x, + "y": y, + "t": t} +# turn above symbolic code into exact solution for dolphin and +# construct the rhs that matches the above exact solution. +exact_solution_example = hlp.generate_exact_solution_expressions( + symbols=symbols, + isRichards=isRichards, + symbolic_pressure=p_e_sym, + symbolic_capillary_pressure=pc_e_sym, + saturation_pressure_relationship=S_pc_sym, + 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, + gravity_acceleration=gravity_acceleration, + include_gravity=include_gravity, + ) +source_expression = exact_solution_example['source'] +exact_solution = exact_solution_example['exact_solution'] +initial_condition = exact_solution_example['initial_condition'] + +# 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() + +# 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) diff --git a/Usecases/InstallationTest/run-simulation b/Usecases/InstallationTest/run-simulation new file mode 100755 index 0000000000000000000000000000000000000000..0eb497502a082a0fec07a5449b1fe946d59c8cc7 --- /dev/null +++ b/Usecases/InstallationTest/run-simulation @@ -0,0 +1,16 @@ +#!/bin/bash + +[ $# -eq 0 ] && { echo "Usage: $0 simulation_file [logfile_name]"; exit 1; } + +SIMULATION_FILE=$1 +SIMULATION=${SIMULATION_FILE%.py} +LOGFILE_DEFAULT="$SIMULATION.log" + +DATE=$(date -I) +LOGFILE=${2:-$DATE-$LOGFILE_DEFAULT} + +GREETING="Simulation $SIMULATION is run on $DATE by $USER" + +echo $GREETING +echo "running $SIMULATION_FILE | tee $LOGFILE" +./$SIMULATION_FILE | tee $LOGFILE