From eb9d6c94be6ec90e75f2bb430e012d8ddf435b30 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Wed, 24 Jun 2020 20:41:12 +0200
Subject: [PATCH] try to achieve parallelistation

---
 LDDsimulation/LDDsimulation.py                |  24 +-
 LDDsimulation/solutionFile.py                 |   4 +-
 ...layered_soil_with_inner_patch-realistic.py | 274 ++++++++++++------
 3 files changed, 199 insertions(+), 103 deletions(-)

diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 146ee46..43ff853 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -88,8 +88,8 @@ class LDDsimulation(object):
         # # no sponge zones
         # self.sponge = False
         # set up mpi for internal parallelisation.
-        self._mpi_rank = df.MPI.rank(df.MPI.comm_world)
-        self._mpi_communicator = df.MPI.comm_world
+        # self._mpi_rank = df.MPI.rank(df.MPI.comm_world)
+        # self._mpi_communicator = df.MPI.comm_world
 
 
         ## Public variables.
@@ -606,8 +606,8 @@ class LDDsimulation(object):
                 filename = self.output_dir+self.output_filename_parameter_part[ind]\
                 + "solution_iterations_at_timestep"\
                 + "{number}".format(number=self.timestep_num) +".xdmf"
-                solution_over_iteration_within_timestep.update( #
-                    {ind: SolutionFile(self._mpi_communicator, filename)}
+                solution_over_iteration_within_timestep.update( #self._mpi_communicator,
+                    {ind: SolutionFile(filename)}
                     )
 
             # reset all interface[has_interface].current_iteration[ind] = 0, for
@@ -881,8 +881,8 @@ class LDDsimulation(object):
             self.output_filename_parameter_part.update(
                 {subdom_ind: filename_param_part}
             )
-            self.solution_file.update( #
-                {subdom_ind: SolutionFile(self._mpi_communicator, filename)}
+            self.solution_file.update( #self._mpi_communicator, 
+                {subdom_ind: SolutionFile(filename)}
                 )
             # self.postprocessed_solution_file.update(
             #     {subdom_ind: SolutionFile(self._mpi_communicator, post_filename)}
@@ -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/LDDsimulation/solutionFile.py b/LDDsimulation/solutionFile.py
index e71450f..22bfdfd 100644
--- a/LDDsimulation/solutionFile.py
+++ b/LDDsimulation/solutionFile.py
@@ -6,8 +6,8 @@ class SolutionFile(df.XDMFFile):
     This class extends `df.XDMFFile` with some minor changes for convenience.
     """
     def __init__(self, mpicomm, filepath):
-        df.XDMFFile.__init__(self, mpicomm, filepath)
-
+        df.XDMFFile.__init__(self, filepath)
+        # mpicomm, 
         self.parameters["functions_share_mesh"] = True
         self.parameters["flush_output"] = True
         self.path = filepath # Mimic the file path attribute from a `file` returned by `open`
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 42c7614..b991ae3 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 = 1
 
 # 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,188 @@ 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(
+        **kwargs
+        # 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.
+    """
+    for key, value in kwargs.items():
+        print ("%s == %s" %(key, value))
+
+    # initialise LDD simulation class
+    simulation = ldd.LDDsimulation(
+        tol=kwargs["tol"],
+        LDDsolver_tol=kwargs["solver_tol"],
+        debug=kwargs["debugflag"],
+        max_iter_num=kwargs["max_iter_num"],
+        FEM_Lagrange_degree=kwargs["FEM_Lagrange_degree"],
+        mesh_study=kwargs["mesh_study"]
+        )
+
+    simulation.set_parameters(
+        use_case=kwargs["use_case"],
+        output_dir=kwargs["output_string"],
+        subdomain_def_points=kwargs["subdomain_def_points"],
+        isRichards=kwargs["isRichards"],
+        interface_def_points=kwargs["interface_def_points"],
+        outer_boundary_def_points=kwargs["outer_boundary_def_points"],
+        adjacent_subdomains=kwargs["adjacent_subdomains"],
+        mesh_resolution=kwargs["mesh_resolution"],
+        viscosity=kwargs["viscosity"],
+        porosity=kwargs["porosity"],
+        L=kwargs["L"],
+        lambda_param=kwargs["lambda_param"],
+        relative_permeability=kwargs["relative_permeability"],
+        saturation=kwargs["sat_pressure_relationship"],
+        starttime=kwargs["starttime"],
+        number_of_timesteps=kwargs["number_of_timesteps"],
+        number_of_timesteps_to_analyse=kwargs["number_of_timesteps_to_analyse"],
+        plot_timestep_every=kwargs["plot_timestep_every"],
+        timestep_size=kwargs["timestep_size"],
+        sources=kwargs["source_expression"],
+        initial_conditions=kwargs["initial_condition"],
+        dirichletBC_expression_strings=kwargs["dirichletBC"],
+        exact_solution=kwargs["exact_solution"],
+        densities=kwargs["densities"],
+        include_gravity=kwargs["include_gravity"],
+        gravity_acceleration=kwargs["gravity_acceleration"],
+        write2file=kwargs["write_to_file"],
+        )
+
+    simulation.initialise()
+    output_dir = simulation.output_dir
+    # simulation.write_exact_solution_to_xdmf()
+    output = simulation.run(analyse_condition=kwargs["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__':
+    # mp.set_start_method('spawn')
+    # 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,
+                 "mesh_resolution": mesh_resolution,
+                 "starttime": starttime}
+            )
+            # info(simulation_parameter["use_case"])
+            # LDDsim = mp.Process(
+            #             target=run_simulation,
+            #             args=(
+            #                 starttime,
+            #                 mesh_resolution,
+            #                 simulation_parameter
+            #                 )
+            #             )
+            # LDDsim.start()
+            # LDDsim.join()
+            run_simulation(
+                # simulation_parameter,
+                # starttime=starttime,
+                kwargs=mesh_resolution
+                )
-- 
GitLab