diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 146ee4634e5730dfb7d4bd06d5c2c530e88437df..bdbb20a198bef5e5060fd22fa8dae6872969f725 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -30,7 +30,7 @@ class LDDsimulation(object):
def __init__(self, dimension: int = 2,#
debug: bool = False,#
- tol: float = None,#
+ tol: float = 1E-14,#
LDDsolver_tol: float = 1E-8,
# maximal number of L-iterations that the LDD solver uses.
max_iter_num: int = 1000,
@@ -1117,9 +1117,9 @@ class LDDsimulation(object):
mesh = mshr.generate_mesh(domain, mesh_resolution)
if self.write2file['meshes_and_markers']:
- filepath = self.output_dir+"global_mesh.pvd"
+ filepath = self.output_dir+f"global_mesh_meshrez{mesh_resolution}.pvd"
df.File(filepath) << mesh
- filepath = self.output_dir+"global_mesh.xml.gz"
+ filepath = self.output_dir+f"global_mesh_meshrez{mesh_resolution}.xml.gz"
df.File(filepath) << mesh
self.mesh_subdomain.append(mesh)
@@ -1139,7 +1139,7 @@ class LDDsimulation(object):
if self.write2file['meshes_and_markers']:
- filepath = self.output_dir+"domain_marker.pvd"
+ filepath = self.output_dir+f"domain_marker_meshrez{mesh_resolution}.pvd"
df.File(filepath) << self.domain_marker
@@ -1214,7 +1214,7 @@ class LDDsimulation(object):
# print(self.interface[glob_interface_num])
if self.write2file['meshes_and_markers']:
- filepath = self.output_dir+"interface_marker_global.pvd"
+ filepath = self.output_dir+f"interface_marker_global_meshrez{self.mesh_resolution}.pvd"
df.File(filepath) << self.interface_marker
@@ -1281,9 +1281,9 @@ class LDDsimulation(object):
# print(f"parameter: {parameter} = {value}")
if self.write2file['meshes_and_markers']:
- filepath = self.output_dir+f"interface_marker_subdomain{subdom_num}.pvd"
+ filepath = self.output_dir+f"interface_marker_subdomain{subdom_num}_meshrez{self.mesh_resolution}.pvd"
df.File(filepath) << self.subdomain[subdom_num].interface_marker
- filepath = self.output_dir+f"outer_boundary_marker_subdomain{subdom_num}.pvd"
+ filepath = self.output_dir+f"outer_boundary_marker_subdomain{subdom_num}_meshrez{self.mesh_resolution}.pvd"
df.File(filepath) << self.subdomain[subdom_num].outer_boundary_marker
# print(self.subdomain[subdom_num])
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 42c7614031ea35a5d2ccd59e5ba2ed9eff6ed94e..405bc4251d86ee640a956a0b052909fcf8742387 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
@@ -12,6 +12,7 @@ import helpers as hlp
import datetime
import os
import pandas as pd
+import multiprocessing as mp
# init sympy session
@@ -37,18 +38,18 @@ thisfile = "TP-R-layered_soil_with_inner_patch-realistic.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
-max_iter_num = 1000
+max_iter_num = 1
FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
-mesh_study = False
+mesh_study = True
resolutions = {
- # 1: 2e-6, # h=2
+ 1: 2e-6, # h=2
# 2: 2e-6, # h=1.1180
# 4: 2e-6, # h=0.5590
- # 8: 2e-6, # h=0.2814
+ 8: 2e-6, # h=0.2814
# 16: 8e-6, # h=0.1412
- 32: 5e-6,
+ # 32: 5e-6,
# 64: 2e-6,
# 128: 2e-6
}
@@ -58,7 +59,7 @@ resolutions = {
# for each element t_0 in starttimes.
starttimes = [0.0]
timestep_size = 0.001
-number_of_timesteps = 1000
+number_of_timesteps = 5
# LDD scheme parameters ######################################################
@@ -108,7 +109,7 @@ lambda56_w = 4
lambda56_nw = 4
include_gravity = True
-debugflag = False
+debugflag = True
analyse_condition = False
# I/O CONFIG #################################################################
@@ -118,7 +119,7 @@ analyse_condition = False
plot_timestep_every = 1
# Decide how many timesteps you want analysed. Analysed means, that
# subsequent errors of the L-iteration within the timestep are written out.
-number_of_timesteps_to_analyse = 5
+number_of_timesteps_to_analyse = 1
# fine grained control over data to be written to disk in the mesh study case
# as well as for a regular simuation for a fixed grid.
@@ -945,93 +946,183 @@ for subdomain in isRichards.keys():
# 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()
+# f = open(thisfile, 'r')
+# 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}
+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
)
- 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:
+ # 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(
- filename,
+ 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())
+
+
+# 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})
+ simulation_parameter.update({"mesh_resolution": mesh_resolution})
+ simulation_parameter.update({"starttime": starttime})
+ info(simulation_parameter["use_case"])
+ LDDsim = mp.Process(
+ target=run_simulation,
+ args=(
+ simulation_parameter,
+ starttime,
+ mesh_resolution
+ )
+ )
+ LDDsim.start()
+ LDDsim.join()
+ # run_simulation(
+ # mesh_resolution=mesh_resolution,
+ # starttime=starttime,
+ # parameter=simulation_parameter
+ # )