Merge pull request #167 from trixi-framework/msl/mpi_parallel

Add support for parallel 2D simulations on static, uniform mesh

.travis.yml

@@ -32,11 +32,12 @@ env:
   global:
     - COVERALLS_PARALLEL=true
   jobs:
-    - TRIXI_TEST=1D
     - TRIXI_TEST=2D
     - TRIXI_TEST=3D
-    - TRIXI_TEST=misc
     - TRIXI_TEST=paper-self-gravitating-gas-dynamics
+    - TRIXI_TEST=parallel_2d
+    - TRIXI_TEST=1D
+    - TRIXI_TEST=misc
 notifications:
   webhooks: https://coveralls.io/webhook
   email: false

Project.toml

@@ -9,6 +9,8 @@ HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

docs/make.jl

@@ -40,6 +40,7 @@ makedocs(
         "Home" => "index.md",
         "Development" => "development.md",
         "Visualization" => "visualization.md",
+        "Parallelization" => "parallelization.md",
         "Style guide" => "styleguide.md",
         "GitHub & Git" => "github-git.md",
         "Reference" => [

docs/src/parallelization.md

@@ -0,0 +1,98 @@
+# Parallelization
+
+## Shared-memory parallelization with threads
+Many compute-intensive loops in Trixi.jl are parallelized using the
+[multi-threading](https://docs.julialang.org/en/v1/manual/multi-threading/)
+support provided by Julia. You can recognize those loops by the
+`Threads.@threads` macro prefixed to them, e.g.,
+```julia
+Threads.@threads for element_id in 1:dg.n_elements
+  ...
+end
+```
+This will statically assign an equal iteration count to each available thread.
+
+To use multi-threading, you need to tell Julia at startup how many threads you
+want to use by either setting the environment variable `JULIA_NUM_THREADS` or by
+providing the `-t/--threads` command line argument. For example, to start Julia
+with four threads, start Julia with
+```bash
+julia -t 4
+```
+If both the environment variable and the command line argument are specified at
+the same time, the latter takes precedence.
+
+
+## Distributed computing with MPI
+In addition to the shared memory parallelization with multi-threading, Trixi.jl
+supports distributed parallelism via
+[MPI.jl](https://github.com/JuliaParallel/MPI.jl), which leverages the Message
+Passing Interface (MPI). MPI.jl comes with its own MPI library binaries such
+that there is no need to install MPI yourself. However, it is also possible to
+instead use an existing MPI installation, which is recommended if you are
+running MPI programs on a cluster or supercomputer
+([see the MPI.jl docs](https://juliaparallel.github.io/MPI.jl/stable/configuration/)
+to find out how to select the employed MPI library).
+
+To start Trixi in parallel with MPI, there are three options:
+
+1. **Run from the REPL with `mpiexec()`:** You can start a parallel execution directly from the
+   REPL by executing
+   ```julia
+   julia> using MPI
+
+   julia> mpiexec() do cmd
+            run(`$cmd -n 3 $(Base.julia_cmd()) --project=. -e 'using Trixi; Trixi.run("examples/2d/parameters.toml")'`)
+          end
+   ```
+   The parameter `-n 3` specifies that Trixi should run with three processes (or
+   *ranks* in MPI parlance) and should be adapted to your available
+   computing resources and problem size. The `$(Base.julia_cmd())` argument
+   ensures that Julia is executed in parallel with the same optimization level
+   etc. as you used for the REPL; if this is unnecessary or undesired, you can
+   also just use `julia`.  Further, if you are not running Trixi from a local
+   clone but have installed it as a package, you need to omit the `--project=.`.
+2. **Run from the command line with `mpiexecjl`:** Alternatively, you can
+   use the `mpiexecjl` script provided by MPI.jl, which allows you to start
+   Trixi in parallel directly from the command line. As a preparation, you need to
+   install the script *once* by running
+   ```julia
+   julia> using MPI
+
+   julia> MPI.install_mpiexecjl(destdir="/somewhere/in/your/PATH")
+   ```
+   Then, to execute Trixi in parallel, execute the following command from your
+   command line:
+   ```bash
+   mpiexecjl -n 3 julia --project=. -e 'using Trixi; Trixi.run("examples/2d/parameters.toml")'
+   ```
+3. **Run interactively with `tmpi` (Linux/MacOS only):** If you are on a
+   Linux/macOS system, you have a third option which lets you run Julia in
+   parallel interactively from the REPL. This comes in handy especially during
+   development, as in contrast to the first two options, it allows to reuse the
+   compilation cache and thus facilitates much faster startup times after the
+   first execution. It requires [tmux](https://github.com/tmux/tmux) and the
+   [OpenMPI](https://www.open-mpi.org) library to be installed before, both of
+   which are usually available through a package manager. Once you have
+   installed both tools, you need to configure MPI.jl to use the OpenMPI for
+   your system, which is explained
+   [here](https://juliaparallel.github.io/MPI.jl/stable/configuration/#Using-a-system-provided-MPI).
+   Then, you can download and install the
+   [tmpi](https://github.com/Azrael3000/tmpi)
+   script by executing
+   ```bash
+   curl https://raw.githubusercontent.com/Azrael3000/tmpi/master/tmpi -o /somewhere/in/your/PATH/tmpi
+   ```
+   Finally, you can start and control multiple Julia REPLs simultaneously by
+   running
+   ```bash
+   tmpi 3 julia --project=.
+   ```
+   This will start Julia inside `tmux` three times and multiplexes all commands
+   you enter in one REPL to all other REPLs (try for yourself to understand what
+   it means). If you have no prior experience with `tmux`, handling the REPL
+   this way feels slightly weird in the beginning. However, there is a lot of
+   documentation for `tmux`
+   [available](https://github.com/tmux/tmux/wiki/Getting-Started) and once you
+   get the hang of it, developing Trixi in parallel becomes much smoother this
+   way.

src/Trixi.jl

@@ -15,13 +15,15 @@ module Trixi
 # Include other packages that are used in Trixi
 # (standard library packages first, other packages next, all of them sorted alphabetically)
 using LinearAlgebra: dot
-using Pkg.TOML: parsefile
+using Pkg.TOML: parsefile, parse
 using Printf: @printf, @sprintf, println
 using Profile: clear_malloc_data
 using Random: seed!
 
 using EllipsisNotation
 using HDF5: h5open, attrs
+import MPI
+using OffsetArrays: OffsetArray, OffsetVector
 using StaticArrays: @MVector, @SVector, MVector, MMatrix, MArray, SVector, SMatrix, SArray
 using TimerOutputs: @notimeit, @timeit, TimerOutput, print_timer, reset_timer!
 using UnPack: @unpack
@@ -37,6 +39,7 @@ export globals
 
 # Include all top-level source files
 include("auxiliary/auxiliary.jl")
+include("parallel/parallel.jl")
 include("equations/equations.jl")
 include("mesh/mesh.jl")
 include("solvers/solvers.jl")
@@ -58,4 +61,9 @@ export flux_central, flux_lax_friedrichs, flux_hll,
 export examples_dir, get_examples, default_example
 
 
+function __init__()
+  init_mpi()
+end
+
+
 end

src/auxiliary/auxiliary.jl

@@ -13,10 +13,24 @@ const parameters = Dict{Symbol,Any}()
 
 
 # Parse parameters file into global dict
-function parse_parameters_file(filename)
+parse_parameters_file(filename) = parse_parameters_file(filename, mpi_parallel())
+function parse_parameters_file(filename, mpi_parallel::Val{false})
   parameters[:default] = parsefile(filename)
   parameters[:default]["parameters_file"] = filename
 end
+function parse_parameters_file(filename, mpi_parallel::Val{true})
+  if mpi_isroot()
+    buffer = read(filename)
+    MPI.Bcast!(Ref(length(buffer)), mpi_root(), mpi_comm())
+    MPI.Bcast!(buffer, mpi_root(), mpi_comm())
+  else
+    count = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())
+    buffer = Vector{UInt8}(undef, count[])
+    MPI.Bcast!(buffer, mpi_root(), mpi_comm())
+  end
+  parameters[:default] = parse(String(buffer))
+  parameters[:default]["parameters_file"] = filename
+end
 
 
 # Return parameter by name, optionally taking a default value and a range of valid values.
@@ -118,7 +132,7 @@ function print_startup_message()
        ██║   ██║  ██║██║██╔╝ ██╗██║
        ╚═╝   ╚═╝  ╚═╝╚═╝╚═╝  ╚═╝╚═╝
     """
-  println(s)
+  mpi_println(s)
 end
 
 

src/auxiliary/containers.jl

@@ -307,3 +307,15 @@ function clear!(c::AbstractContainer)
 
   return c
 end
+
+
+# Helpful overloads for `raw_copy`
+function raw_copy!(c::AbstractContainer, first::Int, last::Int, destination::Int)
+  raw_copy!(c, c, first, last, destination)
+end
+function raw_copy!(target::AbstractContainer, source::AbstractContainer, from::Int, destination::Int)
+  raw_copy!(target, source, from, from, destination)
+end
+function raw_copy!(c::AbstractContainer, from::Int, destination::Int)
+  raw_copy!(c, c, from, from, destination)
+end

src/io/io.jl

@@ -1,19 +1,19 @@
+include("parallel.jl")
 
 # Load restart file and store solution in solver
-function load_restart_file!(dg::AbstractDg, restart_filename)
+load_restart_file!(dg, restart_filename) = load_restart_file!(dg, restart_filename, mpi_parallel())
+function load_restart_file!(dg::AbstractDg, restart_filename, mpi_parallel::Val{false})
   # Create variables to be returned later
   time = NaN
   step = -1
 
   # Open file
   h5open(restart_filename, "r") do file
-    equation = equations(dg)
-
     # Read attributes to perform some sanity checks
     if read(attrs(file)["ndims"]) != ndims(dg)
       error("restart mismatch: ndims in solver differs from value in restart file")
     end
-    if read(attrs(file)["equations"]) != get_name(equation)
+    if read(attrs(file)["equations"]) != get_name(equations(dg))
       error("restart mismatch: equations in solver differs from value in restart file")
     end
     if read(attrs(file)["polydeg"]) != polydeg(dg)
@@ -28,7 +28,7 @@ function load_restart_file!(dg::AbstractDg, restart_filename)
     step = read(attrs(file)["timestep"])
 
     # Read data
-    varnames = varnames_cons(equation)
+    varnames = varnames_cons(equations(dg))
     for v in 1:nvariables(dg)
       # Check if variable name matches
       var = file["variables_$v"]
@@ -37,7 +37,6 @@ function load_restart_file!(dg::AbstractDg, restart_filename)
       end
 
       # Read variable
-      println("Reading variables_$v ($name)...")
       dg.elements.u[v, .., :] = read(file["variables_$v"])
     end
   end
@@ -48,7 +47,10 @@ end
 
 # Save current DG solution with some context information as a HDF5 file for
 # restarting.
-function save_restart_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep)
+save_restart_file(dg, mesh, time, dt, timestep) = save_restart_file(dg, mesh, time, dt, timestep,
+                                                                    mpi_parallel())
+function save_restart_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep,
+                           mpi_parallel::Val{false})
   # Create output directory (if it does not exist)
   output_directory = parameter("output_directory", "out")
   mkpath(output_directory)
@@ -62,22 +64,20 @@ function save_restart_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep)
 
   # Open file (clobber existing content)
   h5open(filename * ".h5", "w") do file
-    equation = equations(dg)
-
     # Add context information as attributes
     attrs(file)["ndims"] = ndims(dg)
-    attrs(file)["equations"] = get_name(equation)
+    attrs(file)["equations"] = get_name(equations(dg))
     attrs(file)["polydeg"] = polydeg(dg)
     attrs(file)["n_vars"] = nvariables(dg)
-    attrs(file)["n_elements"] = dg.n_elements
+    attrs(file)["n_elements"] = dg.n_elements_global
     attrs(file)["mesh_file"] = splitdir(mesh.current_filename)[2]
     attrs(file)["time"] = time
     attrs(file)["dt"] = dt
     attrs(file)["timestep"] = timestep
 
     # Restart files always store conservative variables
     data = dg.elements.u
-    varnames = varnames_cons(equation)
+    varnames = varnames_cons(equations(dg))
 
     # Store each variable of the solution
     for v in 1:nvariables(dg)
@@ -95,6 +95,11 @@ end
 # Save current DG solution with some context information as a HDF5 file for
 # postprocessing.
 function save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep, system="")
+  return save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep, system,
+                            mpi_parallel())
+end
+function save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep, system,
+                            mpi_parallel::Val{false})
   # Create output directory (if it does not exist)
   output_directory = parameter("output_directory", "out")
   mkpath(output_directory)
@@ -112,11 +117,9 @@ function save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep,
 
   # Open file (clobber existing content)
   h5open(filename * ".h5", "w") do file
-    equation = equations(dg)
-
     # Add context information as attributes
     attrs(file)["ndims"] = ndims(dg)
-    attrs(file)["equations"] = get_name(equation)
+    attrs(file)["equations"] = get_name(equations(dg))
     attrs(file)["polydeg"] = polydeg(dg)
     attrs(file)["n_vars"] = nvariables(dg)
     attrs(file)["n_elements"] = dg.n_elements
@@ -130,15 +133,15 @@ function save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep,
                                    valid=["conservative", "primitive"])
     if solution_variables == "conservative"
       data = dg.elements.u
-      varnames = varnames_cons(equation)
+      varnames = varnames_cons(equations(dg))
     else
       # Reinterpret the solution array as an array of conservative variables,
       # compute the primitive variables via broadcasting, and reinterpret the
       # result as a plain array of floating point numbers
       data = Array(reinterpret(eltype(dg.elements.u),
              cons2prim.(reinterpret(SVector{nvariables(dg),eltype(dg.elements.u)}, dg.elements.u),
                         Ref(equations(dg)))))
-      varnames = varnames_prim(equation)
+      varnames = varnames_prim(equations(dg))
     end
 
     # Store each variable of the solution
@@ -165,7 +168,8 @@ end
 
 
 # Save current mesh with some context information as an HDF5 file.
-function save_mesh_file(mesh::TreeMesh, timestep=-1)
+save_mesh_file(mesh, timestep=-1) = save_mesh_file(mesh, timestep, mpi_parallel())
+function save_mesh_file(mesh::TreeMesh, timestep, mpi_parallel::Val{false})
   # Create output directory (if it does not exist)
   output_directory = parameter("output_directory", "out")
   mkpath(output_directory)