diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 9c34840afc9f7c0e16cc3b92e686f3bdfe3229cb..032e7a1697ae60ce405d64134d7a62206520a172 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -464,14 +464,16 @@ class LDDsimulation(object):
# step on each subdomain.
# Here it should be possible to parallelise in
# the future.
- for sd_index in self.subdomain.keys():
+ for sd_index, subdomain in self.subdomain.items():
# set subdomain iteration to current iteration
+ # mp.set_start_method('fork')
+ processQueue = mp.Queue()
Lsolver_step = mp.Process(
target=self.run_Lsolver_step,
args=(
iteration,
sd_index,
- max_subsequent_error,
+ processQueue,
solution_over_iteration_within_timestep,
debug,
analyse_timestep,
@@ -479,16 +481,23 @@ class LDDsimulation(object):
)
)
Lsolver_step.start()
+ max_subsequent_error[sd_index-1] = processQueue.get()
+ print(Lsolver_step.exitcode)
Lsolver_step.join()
- # self.encapsulated_Lsolver_step(
- # analyse_timestep=analyse_timestep,
- # analyse_condition=analyse_condition,
- # debug=debug,
+ print(Lsolver_step.exitcode)
+ # update interface iteration numbers
+ # for index in subdomain.has_interface:
+ # interface = self.interface[index]
+ # interface.current_iteration[sd_index] = iteration
+ # Lsolver_step.terminate()
+ # self.run_Lsolver_step(
# iteration=iteration,
# subdomain_index=sd_index,
- # subdomain=subdomain,
# max_subsequent_error=max_subsequent_error,
# sol_over_iter_file=solution_over_iteration_within_timestep,
+ # analyse_timestep=analyse_timestep,
+ # analyse_condition=analyse_condition,
+ # debug=debug
# )
# end loop over subdomains
# stopping criterion for the solver.
@@ -513,7 +522,7 @@ class LDDsimulation(object):
self,
iteration: int,
subdomain_index: int,
- max_subsequent_error: np.ndarray,
+ processQueue,
sol_over_iter_file: tp.Dict[int, tp.Type[SolutionFile]],
debug: bool = False,
analyse_timestep: bool = False,
@@ -557,10 +566,6 @@ class LDDsimulation(object):
"iteration {} = ".format(iteration)
print(debugstring, subsequent_iter_error[phase])
- # calculate the maximum over phase of subsequent errors for the
- # stopping criterion.
- max_subsequent_error[sd_index-1] = max(subsequent_iter_error.values())
-
if analyse_timestep:
# save the calculated L2 errors to csv file for plotting
# with pgfplots
@@ -608,6 +613,12 @@ class LDDsimulation(object):
subdomain.pressure[phase]
)
+ # calculate the maximum over phase of subsequent errors for the
+ # stopping criterion.
+ # max_subsequent_error[sd_index-1] = max(subsequent_iter_error.values())
+ processQueue.put(max(subsequent_iter_error.values()))
+
+
def prepare_LDDsolver(
self,
time: float = None,
@@ -651,23 +662,23 @@ class LDDsimulation(object):
# before the iteration gets started all iteration numbers must be
# nulled
for subdomain_index in self.subdomain.keys():
- prepare_subdomain = mp.Process(
- target=self.prepare_subdomain,
- args=(
- subdomain_index,
- solution_over_iteration_within_timestep,
- debug,
- analyse_timestep,
- )
- )
- prepare_subdomain.start()
- prepare_subdomain.join()
- # self.prepare_subdomain(
- # subdomain_index=ind,
- # sol_over_iter_file=solution_over_iteration_within_timestep,
- # debug=debug,
- # analyse_timestep=analyse_timestep,
- # )
+ # prepare_subdomain = mp.Process(
+ # target=self.prepare_subdomain,
+ # args=(
+ # subdomain_index,
+ # solution_over_iteration_within_timestep,
+ # debug,
+ # analyse_timestep,
+ # )
+ # )
+ # prepare_subdomain.start()
+ # prepare_subdomain.join()
+ self.prepare_subdomain(
+ subdomain_index=subdomain_index,
+ sol_over_iter_file=solution_over_iteration_within_timestep,
+ debug=debug,
+ analyse_timestep=analyse_timestep,
+ )
return solution_over_iteration_within_timestep
@@ -772,7 +783,6 @@ class LDDsimulation(object):
f"pressure_prev_timestep[{phase}].vector()[:]" + " =\n" + \
f"{subdomain.pressure_prev_timestep[phase].vector()[:]}"
print(debugstr5)
-
# we need to update the previous pressure dictionaries of all
# interfaces of the subdomain with the current solution, so that the
# subdomain.calc_gli_term method works properly.
@@ -802,7 +812,6 @@ class LDDsimulation(object):
to write out each iteration.
"""
subdomain = self.subdomain[subdomain_index]
- # iteration = subdomain.iteration_number
condition_number = None
if analyse_condition:
condition_number = dict()
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 5dd7b983bd292ded3b33fab1d014bea264770977..f9aaec372c6005a0fbe29bd29e617ebe096dcb83 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
@@ -60,7 +60,7 @@ resolutions = {
# for each element t_0 in starttimes.
starttimes = [0.0]
timestep_size = 0.001
-number_of_timesteps = 500
+number_of_timesteps = 10
# LDD scheme parameters ######################################################
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 b807a25c4c6aebf3e0b21051e74af6f545ae06ed..9360a1fd371fc7fc2f0f9dd0001bfcc17f2bfb41 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
@@ -12,6 +12,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()
@@ -36,20 +38,20 @@ thisfile = "TP-R-2-patch-realistic.py"
# GENERAL SOLVER CONFIG ######################################################
# maximal iteration per timestep
-max_iter_num = 500
+max_iter_num = 50
FEM_Lagrange_degree = 1
# GRID AND MESH STUDY SPECIFICATIONS #########################################
-mesh_study = True
+mesh_study = False
resolutions = {
- 1: 1e-5,
- 2: 1e-5,
- 4: 1e-5,
- 8: 1e-5,
+ # 1: 1e-5,
+ # 2: 1e-5,
+ # 4: 1e-5,
+ # 8: 1e-5,
16: 5e-6,
- 32: 5e-6,
- 64: 2e-6,
- 128: 2e-6,
+ # 32: 5e-6,
+ # 64: 2e-6,
+ # 128: 2e-6,
# 256: 1e-6,
}
@@ -58,19 +60,19 @@ resolutions = {
# for each element t_0 in starttimes.
starttimes = [0.0]
timestep_size = 0.001
-number_of_timesteps = 800
+number_of_timesteps = 5
# LDD scheme parameters ######################################################
-Lw1 = 0.25 #/timestep_size
-Lnw1= 0.25
+Lw1 = 0.025 #/timestep_size
+Lnw1= 0.025
-Lw2 = 0.5 #/timestep_size
-Lnw2= 0.25
+Lw2 = 0.025 #/timestep_size
+Lnw2= 0.025
-lambda_w = 40
-lambda_nw = 40
+lambda_w = 44
+lambda_nw = 4
-include_gravity = True
+include_gravity = False
debugflag = True
analyse_condition = False
@@ -431,14 +433,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,
@@ -464,35 +462,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
@@ -502,88 +483,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
- )
- 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,
+ "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
+ # )
+ # )
+ # LDDsim.start()
+ # LDDsim.join()
+ hlp.run_simulation(
+ mesh_resolution=mesh_resolution,
+ starttime=starttime,
+ parameter=simulation_parameter
+ )