diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index bdbb20a198bef5e5060fd22fa8dae6872969f725..e250eeb30fc9fc34f09d81eb1a77e241375de1af 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -15,6 +15,7 @@ import helpers as hlp
import pandas as pd
import sys
import copy
+# import multiprocessing as mp
# import errno
import h5py
from termcolor import colored
@@ -204,37 +205,37 @@ class LDDsimulation(object):
# minimal error has been achieved, i.e. the calculation can be considered
# done or not.
- def set_parameters(self,#
- use_case: str,
- output_dir: str,#
- subdomain_def_points: tp.List[tp.List[df.Point]],#
- # this is actually an array of bools
- isRichards: tp.Dict[int, bool],#
- interface_def_points: tp.List[tp.List[df.Point]],#
- outer_boundary_def_points: tp.Dict[int, tp.Dict[int, tp.List[df.Point]]],#
- adjacent_subdomains: tp.List[np.ndarray],#
- mesh_resolution: float,#
- viscosity: tp.Dict[int, tp.List[float]],#
- porosity: tp.Dict[int, tp.Dict[str, float]],#
- L: tp.Dict[int, tp.Dict[str, float]],#
- lambda_param: tp.Dict[int, tp.Dict[str, float]],#
- relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
- saturation: tp.Dict[int, tp.Callable[...,None]],#
- starttime: float,#
- timestep_size: tp.Dict[int, float],#
- number_of_timesteps: int,#
- number_of_timesteps_to_analyse: int,#
- sources: tp.Dict[int, tp.Dict[str, str]],#
- initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
- dirichletBC_expression_strings: tp.Dict[int, tp.Dict[str, str]],#
- write2file: tp.Dict[str, bool],#
- exact_solution: tp.Dict[int, tp.Dict[str, str]] = None, #
- densities: tp.Dict[int, tp.Dict[str, float]] = None,#
- include_gravity: bool = False,#
- gravity_acceleration: float = 9.81,
- plot_timestep_every: int = 1,
-
- )-> None:
+ def set_parameters(
+ self,#
+ use_case: str,
+ output_dir: str,#
+ subdomain_def_points: tp.List[tp.List[df.Point]],#
+ # this is actually an array of bools
+ isRichards: tp.Dict[int, bool],#
+ interface_def_points: tp.List[tp.List[df.Point]],#
+ outer_boundary_def_points: tp.Dict[int, tp.Dict[int, tp.List[df.Point]]],#
+ adjacent_subdomains: tp.List[np.ndarray],#
+ mesh_resolution: float,#
+ viscosity: tp.Dict[int, tp.List[float]],#
+ porosity: tp.Dict[int, tp.Dict[str, float]],#
+ L: tp.Dict[int, tp.Dict[str, float]],#
+ lambda_param: tp.Dict[int, tp.Dict[str, float]],#
+ relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
+ saturation: tp.Dict[int, tp.Callable[...,None]],#
+ starttime: float,#
+ timestep_size: tp.Dict[int, float],#
+ number_of_timesteps: int,#
+ number_of_timesteps_to_analyse: int,#
+ sources: tp.Dict[int, tp.Dict[str, str]],#
+ initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
+ dirichletBC_expression_strings: tp.Dict[int, tp.Dict[str, str]],#
+ write2file: tp.Dict[str, bool],#
+ exact_solution: tp.Dict[int, tp.Dict[str, str]] = None, #
+ densities: tp.Dict[int, tp.Dict[str, float]] = None,#
+ include_gravity: bool = False,#
+ gravity_acceleration: float = 9.81,
+ plot_timestep_every: int = 1,
+ )-> None:
""" set parameters of an instance of class LDDsimulation"""
self.use_case = use_case
@@ -486,7 +487,12 @@ class LDDsimulation(object):
{phase: error_calculated}
)
if debug:
- print(f"time = {time} (at timestep {self.timestep_num}): subsequent error on subdomain {sd_index} for {phase} phase in iteration {iteration} = ", subsequent_iter_error[phase])
+ debugstring = "time = {} (at timestep".format(time) +\
+ "{}): subsequent ".formal(self.timestep_num) +\
+ "error on subdomain {} ".format(sd_index) +\
+ "for {} phase in ".format(phase) +\
+ "iteration {} = ".format(iteration)
+ print(debugstring, subsequent_iter_error[phase])
# calculate the maximum over phase of subsequent errors for the
# stopping criterion.
diff --git a/LDDsimulation/helpers.py b/LDDsimulation/helpers.py
index efcaf68185d27d7afc70960804cb3fe38adfdcf2..a2a937406b04da916c5062c8404bb87f65072742 100644
--- a/LDDsimulation/helpers.py
+++ b/LDDsimulation/helpers.py
@@ -10,10 +10,125 @@ import matplotlib.pyplot as plt
import numpy as np
import typing as tp
import sympy as sym
+import os
+import pandas as pd
+import LDDsimulation as ldd
+
+
+# RUN #########################################################################
+def run_simulation(
+ parameter,
+ starttime: float = 0.0,
+ mesh_resolution: float = 32,
+ ) -> None:
+ """Setup and run LDD simulation
+
+ Setup and run LDD simulation for starttime starttime and mesh resolution
+ mesh_resolution.
+
+ PARAMETERS
+ starttime #type float number >= 0 specifying the starting point
+ of the simulation.
+ mesh_resolution #type float mesh resolution parameter for Fenics.
+ parameter #type tp.Dict Dictionary of all parameters of the LDD
+ simulation.
+ """
+ # parameter = kwargs
+ # initialise LDD simulation class
+ simulation = ldd.LDDsimulation(
+ tol=parameter["tol"],
+ LDDsolver_tol=parameter["solver_tol"],
+ debug=parameter["debugflag"],
+ max_iter_num=parameter["max_iter_num"],
+ FEM_Lagrange_degree=parameter["FEM_Lagrange_degree"],
+ mesh_study=parameter["mesh_study"]
+ )
+
+ simulation.set_parameters(
+ use_case=parameter["use_case"],
+ output_dir=parameter["output_string"],
+ subdomain_def_points=parameter["subdomain_def_points"],
+ isRichards=parameter["isRichards"],
+ interface_def_points=parameter["interface_def_points"],
+ outer_boundary_def_points=parameter["outer_boundary_def_points"],
+ adjacent_subdomains=parameter["adjacent_subdomains"],
+ mesh_resolution=mesh_resolution, # parameter["mesh_resolution"],
+ viscosity=parameter["viscosity"],
+ porosity=parameter["porosity"],
+ L=parameter["L"],
+ lambda_param=parameter["lambda_param"],
+ relative_permeability=parameter["relative_permeability"],
+ saturation=parameter["sat_pressure_relationship"],
+ starttime=starttime, # parameter["starttime"],
+ number_of_timesteps=parameter["number_of_timesteps"],
+ number_of_timesteps_to_analyse=parameter["number_of_timesteps_to_analyse"],
+ plot_timestep_every=parameter["plot_timestep_every"],
+ timestep_size=parameter["timestep_size"],
+ sources=parameter["source_expression"],
+ initial_conditions=parameter["initial_condition"],
+ dirichletBC_expression_strings=parameter["dirichletBC"],
+ exact_solution=parameter["exact_solution"],
+ densities=parameter["densities"],
+ include_gravity=parameter["include_gravity"],
+ gravity_acceleration=parameter["gravity_acceleration"],
+ write2file=parameter["write_to_file"],
+ )
+
+ simulation.initialise()
+ output_dir = simulation.output_dir
+ # simulation.write_exact_solution_to_xdmf()
+ output = simulation.run(analyse_condition=parameter["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_meshrez{}.csv".format(
+ phase,
+ mesh_resolution
+ )
+ # 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
+ )
+ else:
+ errors.to_csv(
+ filename,
+ sep='\t',
+ encoding='utf-8',
+ index=False
+ )
+
+
+def info(title):
+ """show process ids of multiprocessing simulation"""
+ print(title)
+ print('module name:', __name__)
+ print('parent process:', os.getppid())
+ print('process id:', os.getpid())
-# from termcolor import colored
-# from colormap import Colormap
-#import fenicstools as ft
def print_once(*args,**kwargs):
"""When running in mpi only print arguments only once."""
@@ -23,22 +138,22 @@ def print_once(*args,**kwargs):
def generate_exact_solution_expressions(
- symbols: tp.Dict[int, int],
- isRichards: tp.Dict[int, bool],
- symbolic_pressure: tp.Dict[int, tp.Dict[str, str]],
- symbolic_capillary_pressure: tp.Dict[int, str],
- viscosity: tp.Dict[int, tp.List[float]],#
- porosity: tp.Dict[int, tp.Dict[str, float]],#
- relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
- relative_permeability_prime: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],
- saturation_pressure_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
- saturation_pressure_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
- symbolic_S_pc_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
- symbolic_S_pc_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
- densities: tp.Dict[int, tp.Dict[str, float]] = None,#
- gravity_acceleration: float = 9.81,
- include_gravity: bool = False,
- ) -> tp.Dict[str, tp.Dict[int, tp.Dict[str, str]] ]:
+ symbols: tp.Dict[int, int],
+ isRichards: tp.Dict[int, bool],
+ symbolic_pressure: tp.Dict[int, tp.Dict[str, str]],
+ symbolic_capillary_pressure: tp.Dict[int, str],
+ viscosity: tp.Dict[int, tp.List[float]],#
+ porosity: tp.Dict[int, tp.Dict[str, float]],#
+ relative_permeability: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],#
+ relative_permeability_prime: tp.Dict[int, tp.Dict[str, tp.Callable[...,None]] ],
+ saturation_pressure_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
+ saturation_pressure_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
+ symbolic_S_pc_relationship: tp.Dict[int, tp.Callable[...,None]] = None,#
+ symbolic_S_pc_relationship_prime: tp.Dict[int, tp.Callable[...,None]] = None,#
+ densities: tp.Dict[int, tp.Dict[str, float]] = None,#
+ gravity_acceleration: float = 9.81,
+ include_gravity: bool = False,
+ ) -> tp.Dict[str, tp.Dict[int, tp.Dict[str, str]] ]:
""" Generate exact solution, source and initial condition code from symbolic expressions"""
output = {
@@ -133,3 +248,60 @@ def generate_exact_solution_expressions(
)
return output
+
+
+def generate_exact_DirichletBC(
+ isRichards: tp.Dict[int, bool],
+ outer_boundary_def_points: tp.Dict[int, tp.Dict[int, tp.List[df.Point]]],
+ exact_solution: tp.Dict[int, tp.Dict[str, str]]
+ ) -> tp.Dict[int, tp.Dict[str, str]]:
+ """generate Dirichlet BC from exact solution"""
+ 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]}
+ )
+
+ return dirichletBC
+
+
+def generate_exact_symbolic_pc(
+ isRichards: tp.Dict[int, bool],
+ symbolic_pressure: tp.Dict[int, tp.Dict[str, str]]
+ ) -> tp.Dict[str, tp.Dict[str, str]]:
+ """generate symbolic capillary pressure from exact solution"""
+
+ symbolic_pc = dict()
+ for subdomain, isR in isRichards.items():
+ if isR:
+ symbolic_pc.update(
+ {subdomain: -symbolic_pressure[subdomain]['wetting'].copy()}
+ )
+ else:
+ symbolic_pc.update(
+ {subdomain: symbolic_pressure[subdomain]['nonwetting'].copy()
+ - symbolic_pressure[subdomain]['wetting'].copy()}
+ )
+
+ return symbolic_pc
diff --git a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
index 405bc4251d86ee640a956a0b052909fcf8742387..fb77415e1106e64f7b4456a745984e499f1307dd 100755
--- a/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
+++ b/Two-phase-Richards/multi-patch/layered_soil_with_inner_patch/TP-R-layered_soil_with_inner_patch-realistic.py
@@ -13,6 +13,7 @@ import datetime
import os
import pandas as pd
import multiprocessing as mp
+import domainSubstructuring as dss
# init sympy session
@@ -162,216 +163,12 @@ output_string = "./output/{}-{}_timesteps{}_P{}".format(
datestr, use_case, number_of_timesteps, FEM_Lagrange_degree
)
-
# DOMAIN AND INTERFACES #######################################################
-# global domain
-subdomain0_vertices = [df.Point(-1.0,-1.0), #
- df.Point(1.0,-1.0),#
- df.Point(1.0,1.0),#
- df.Point(-1.0,1.0)]
-
-interface12_vertices = [df.Point(-1.0, 0.8),
- df.Point(0.3, 0.8),
- df.Point(0.5, 0.9),
- df.Point(0.8, 0.7),
- df.Point(1.0, 0.65)]
-
-
- # interface23
-interface23_vertices = [df.Point(-1.0, 0.0),
- df.Point(-0.35, 0.0),
- # df.Point(6.5, 4.5),
- df.Point(0.0, 0.0)]
-
-interface24_vertices = [interface23_vertices[2],
- df.Point(0.6, 0.0),
- ]
-
-interface25_vertices = [interface24_vertices[1],
- df.Point(1.0, 0.0)
- ]
-
-
-interface32_vertices = [interface23_vertices[2],
- interface23_vertices[1],
- interface23_vertices[0]]
-
-
-interface36_vertices = [df.Point(-1.0, -0.6),
- df.Point(-0.6, -0.45)]
-
-
-interface46_vertices = [interface36_vertices[1],
- df.Point(0.3, -0.25)]
-
-interface56_vertices = [interface46_vertices[1],
- df.Point(0.65, -0.6),
- df.Point(1.0, -0.7)]
-
-
-
-
-interface34_vertices = [interface36_vertices[1],
- interface23_vertices[2]]
-
-# Interface 45 needs to be split, because of the shape. There can be triangles
-# with two facets on the interface and this creates a rogue dof type error when
-# integrating over that particular interface. Accordingly, the lambda_param
-# dictionary has two entries for that interface.
-interface45_vertices_a = [interface56_vertices[0],
- df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
- ]
-interface45_vertices_b = [df.Point(0.7, -0.2),#df.Point(0.7, -0.2),
- interface25_vertices[0]
- ]
-
-# interface_vertices introduces a global numbering of interfaces.
-interface_def_points = [interface12_vertices,
- interface23_vertices,
- interface24_vertices,
- interface25_vertices,
- interface34_vertices,
- interface36_vertices,
- interface45_vertices_a,
- interface45_vertices_b,
- interface46_vertices,
- interface56_vertices,
- ]
-adjacent_subdomains = [[1,2],
- [2,3],
- [2,4],
- [2,5],
- [3,4],
- [3,6],
- [4,5],
- [4,5],
- [4,6],
- [5,6]
- ]
-
-# subdomain1.
-subdomain1_vertices = [interface12_vertices[0],
- interface12_vertices[1],
- interface12_vertices[2],
- interface12_vertices[3],
- interface12_vertices[4], # southern boundary, 12 interface
- subdomain0_vertices[2], # eastern boundary, outer boundary
- subdomain0_vertices[3]] # northern boundary, outer on_boundary
-
-# vertex coordinates of the outer boundaries. If it can not be specified as a
-# polygon, use an entry per boundary polygon. This information is used for defining
-# the Dirichlet boundary conditions. If a domain is completely internal, the
-# dictionary entry should be 0: None
-subdomain1_outer_boundary_verts = {
- 0: [subdomain1_vertices[4], #
- subdomain1_vertices[5], # eastern boundary, outer boundary
- subdomain1_vertices[6],
- subdomain1_vertices[0]]
-}
-
-#subdomain1
-subdomain2_vertices = [interface23_vertices[0],
- interface23_vertices[1],
- interface23_vertices[2],
- interface24_vertices[1],
- interface25_vertices[1], # southern boundary, 23 interface
- subdomain1_vertices[4], # eastern boundary, outer boundary
- subdomain1_vertices[3],
- subdomain1_vertices[2],
- subdomain1_vertices[1],
- subdomain1_vertices[0] ] # northern boundary, 12 interface
-
-subdomain2_outer_boundary_verts = {
- 0: [subdomain2_vertices[9],
- subdomain2_vertices[0]],
- 1: [subdomain2_vertices[4],
- subdomain2_vertices[5]]
-}
-
-
-subdomain3_vertices = [interface36_vertices[0],
- interface36_vertices[1],
- # interface34_vertices[0],
- interface34_vertices[1],
- # interface32_vertices[0],
- interface32_vertices[1],
- interface32_vertices[2]
- ]
-
-subdomain3_outer_boundary_verts = {
- 0: [subdomain3_vertices[4],
- subdomain3_vertices[0]]
-}
-
-
-# subdomain3
-subdomain4_vertices = [interface46_vertices[0],
- interface46_vertices[1],
- # interface45_vertices[1],
- interface45_vertices_a[1],
- interface24_vertices[1],
- interface24_vertices[0],
- interface34_vertices[1]
- ]
-
-subdomain4_outer_boundary_verts = None
-
-subdomain5_vertices = [interface56_vertices[0],
- interface56_vertices[1],
- interface56_vertices[2],
- interface25_vertices[1],
- interface25_vertices[0],
- interface45_vertices_b[1],
- interface45_vertices_b[0]
-]
-
-
-subdomain5_outer_boundary_verts = {
- 0: [subdomain5_vertices[2],
- subdomain5_vertices[3]]
-}
-
-
-
-subdomain6_vertices = [subdomain0_vertices[0],
- subdomain0_vertices[1], # southern boundary, outer boundary
- interface56_vertices[2],
- interface56_vertices[1],
- interface56_vertices[0],
- interface36_vertices[1],
- interface36_vertices[0]
- ]
-
-subdomain6_outer_boundary_verts = {
- 0: [subdomain6_vertices[6],
- subdomain6_vertices[0],
- subdomain6_vertices[1],
- subdomain6_vertices[2]]
-}
-
-
-subdomain_def_points = [subdomain0_vertices,#
- subdomain1_vertices,#
- subdomain2_vertices,#
- subdomain3_vertices,#
- subdomain4_vertices,
- subdomain5_vertices,
- subdomain6_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,
- 3: subdomain3_outer_boundary_verts,
- 4: subdomain4_outer_boundary_verts,
- 5: subdomain5_outer_boundary_verts,
- 6: subdomain6_outer_boundary_verts
-}
+substructuring = dss.layeredSoilInnerPatch()
+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 = {
@@ -879,15 +676,10 @@ p_e_sym = {
'nonwetting': (-1.0 -t*(1.0 + x**2) - sym.sqrt(2+t**2)**2)*y**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,
@@ -914,154 +706,24 @@ 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
# log file with ./TP-R-layered_soil.py | tee simulation.log
-# f = open(thisfile, 'r')
-# print(f.read())
-# f.close()
-
-
-# RUN #########################################################################
-def run_simulation(
- parameter,
- starttime: float = 0.0,
- mesh_resolution: float = 32,
- ) -> None:
- """Setup and run LDD simulation
-
- Setup and run LDD simulation for starttime starttime and mesh resolution
- mesh_resolution.
-
- PARAMETERS
- starttime #type float number >= 0 specifying the starting point
- of the simulation.
- mesh_resolution #type float mesh resolution parameter for Fenics.
- """
- # parameter = kwargs
- # initialise LDD simulation class
- simulation = ldd.LDDsimulation(
- tol=parameter["tol"],
- LDDsolver_tol=parameter["solver_tol"],
- debug=parameter["debugflag"],
- max_iter_num=parameter["max_iter_num"],
- FEM_Lagrange_degree=parameter["FEM_Lagrange_degree"],
- mesh_study=parameter["mesh_study"]
- )
-
- simulation.set_parameters(
- use_case=parameter["use_case"],
- output_dir=parameter["output_string"],
- subdomain_def_points=parameter["subdomain_def_points"],
- isRichards=parameter["isRichards"],
- interface_def_points=parameter["interface_def_points"],
- outer_boundary_def_points=parameter["outer_boundary_def_points"],
- adjacent_subdomains=parameter["adjacent_subdomains"],
- mesh_resolution=mesh_resolution, # parameter["mesh_resolution"],
- viscosity=parameter["viscosity"],
- porosity=parameter["porosity"],
- L=parameter["L"],
- lambda_param=parameter["lambda_param"],
- relative_permeability=parameter["relative_permeability"],
- saturation=parameter["sat_pressure_relationship"],
- starttime=starttime, # parameter["starttime"],
- number_of_timesteps=parameter["number_of_timesteps"],
- number_of_timesteps_to_analyse=parameter["number_of_timesteps_to_analyse"],
- plot_timestep_every=parameter["plot_timestep_every"],
- timestep_size=parameter["timestep_size"],
- sources=parameter["source_expression"],
- initial_conditions=parameter["initial_condition"],
- dirichletBC_expression_strings=parameter["dirichletBC"],
- exact_solution=parameter["exact_solution"],
- densities=parameter["densities"],
- include_gravity=parameter["include_gravity"],
- gravity_acceleration=parameter["gravity_acceleration"],
- write2file=parameter["write_to_file"],
- )
-
- simulation.initialise()
- output_dir = simulation.output_dir
- # simulation.write_exact_solution_to_xdmf()
- output = simulation.run(analyse_condition=parameter["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_meshrez{}.csv".format(
- phase,
- mesh_resolution
- )
- # 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
- )
- else:
- errors.to_csv(
- filename,
- sep='\t',
- encoding='utf-8',
- index=False
- )
-
-
-def info(title):
- """show process ids of multiprocessing simulation"""
- print(title)
- print('module name:', __name__)
- print('parent process:', os.getppid())
- print('process id:', os.getpid())
+f = open(thisfile, 'r')
+print(f.read())
+f.close()
# MAIN ########################################################################
@@ -1108,11 +770,9 @@ if __name__ == '__main__':
for starttime in starttimes:
for mesh_resolution, solver_tol in resolutions.items():
simulation_parameter.update({"solver_tol": solver_tol})
- simulation_parameter.update({"mesh_resolution": mesh_resolution})
- simulation_parameter.update({"starttime": starttime})
- info(simulation_parameter["use_case"])
+ hlp.info(simulation_parameter["use_case"])
LDDsim = mp.Process(
- target=run_simulation,
+ target=hlp.run_simulation,
args=(
simulation_parameter,
starttime,