Merge MPI into Taal

.github/workflows/TagBot.yml

@@ -1,10 +1,12 @@
 name: TagBot
 on:
-  schedule:
-    - cron: 0 0 * * *
+  issue_comment:
+    types:
+      - created
   workflow_dispatch:
 jobs:
   TagBot:
+    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
     runs-on: ubuntu-latest
     steps:
       - uses: JuliaRegistries/TagBot@v1

.gitignore

@@ -10,6 +10,7 @@
 **/Manifest.toml
 out*/
 docs/build
+docs/src/authors.md
 public/
 coverage/
 coverage_report/

.travis.yml

@@ -37,6 +37,7 @@ env:
     - TRIXI_TEST=3D
     - TRIXI_TEST=3D_OLD
     - TRIXI_TEST=paper-self-gravitating-gas-dynamics
+    - TRIXI_TEST=parallel_2d
     - TRIXI_TEST=1D
     - TRIXI_TEST=misc
 notifications:

Project.toml

@@ -11,6 +11,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"
@@ -27,6 +29,8 @@ EllipsisNotation = "1.0"
 HDF5 = "0.13"
 LinearMaps = "2.7"
 LoopVectorization = "0.8.25"
+MPI = "0.15"
+OffsetArrays = "1.3"
 StaticArrays = "0.12"
 TimerOutputs = "0.5"
 Tullio = "0.2.5"

docs/make.jl

@@ -6,6 +6,11 @@ using Trixi2Img
 # Get Trixi root directory
 trixi_root_dir = dirname(@__DIR__)
 
+# Copy list of authors to not need to synchronize it manually
+authors_text = read(joinpath(trixi_root_dir, "AUTHORS.md"), String)
+authors_text = replace(authors_text, "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)")
+write(joinpath(@__DIR__, "src", "authors.md"), authors_text)
+
 # Install dependencies and import modules...
 # ...Trixi
 Pkg.activate(trixi_root_dir)
@@ -42,6 +47,7 @@ makedocs(
         "Time integration" => "time_integration.md",
         "Development" => "development.md",
         "Visualization" => "visualization.md",
+        "Parallelization" => "parallelization.md",
         "Style guide" => "styleguide.md",
         "GitHub & Git" => "github-git.md",
         "Reference" => [

docs/src/index.md

@@ -24,7 +24,7 @@ installation and postprocessing procedures. Its features include:
 * Multi-physics simulations
   * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics)
 * Shared-memory parallelization via multithreading
-* Visualization of results with Julia-only tools (2D) or ParaView/VisIt (2D/3D)
+* Visualization of results with Julia-only tools ([Trixi2Img](https://github.com/trixi-framework/Trixi2Img.jl)) or ParaView/VisIt ([Trixi2Vtk](https://github.com/trixi-framework/Trixi2Vtk.jl))
 
 
 ## Installation

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.

docs/src/visualization.md

@@ -85,6 +85,15 @@ that can be opened with any image viewer:
 
 !["solution_000040_scalar_resized"](https://github.com/trixi-framework/Trixi2Img.jl/raw/master/docs/src/assets/solution_000040_scalar_resized.png)
 
+To visualize 3D data generated by Trixi.jl, a 2D slice must be extracted.
+A slice can only lie in a plane orthogonal to one of the coordinate axes.
+The slice plane is defined by the axis to which it is orthogonal and an axis intercept.
+
+For example, to create a 2D slice in the xy-plane at the axis intercept `z = 0.5` from 3D data, execute
+```julia
+julia> trixi2img(joinpath("out", "solution_000000.h5"), output_directory="out", grid_lines=true, slice_axis=:z, slice_axis_intercept=0.5)
+```
+
 Similar to [Trixi2Vtk](@ref), you can also provide multiple files to
 `trixi2img` or use file globbing, e.g.,
 ```julia

src/Trixi.jl

@@ -15,7 +15,7 @@ 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!
@@ -25,6 +25,8 @@ using DiffEqCallbacks: CallbackSet, DiscreteCallback
 using EllipsisNotation # ..
 using HDF5: h5open, attrs
 using LinearMaps: LinearMap
+import MPI
+using OffsetArrays: OffsetArray, OffsetVector
 using StaticArrays: @MVector, @SVector, MVector, MMatrix, MArray, SVector, SMatrix, SArray
 using TimerOutputs: @notimeit, @timeit, @timeit_debug, TimerOutput, print_timer, reset_timer!
 using UnPack: @unpack
@@ -50,6 +52,7 @@ include("basic_types.jl")
 
 # 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")
@@ -116,4 +119,9 @@ export trixi_include, examples_dir, get_examples, default_example
 export convergence_test, jacobian_fd, linear_structure
 
 
+function __init__()
+  init_mpi()
+end
+
+
 end

src/auxiliary/auxiliary.jl

@@ -49,10 +49,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.
@@ -154,7 +168,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/callbacks/save_restart.jl

@@ -97,8 +97,8 @@ function load_mesh(restart_file::AbstractString; n_cells_max)
     read(attrs(file)["ndims"])
   end
 
-  # Note: In the future, we should also read the mesh type
-  mesh = TreeMesh{ndims_}(n_cells_max)
+  # TODO: MPI, we should also read the mesh type
+  mesh = TreeMesh(SerialTree{ndims_}, n_cells_max)
   load_mesh!(mesh, restart_file)
 end
 

src/equations/2d/linear_scalar_advection.jl

@@ -31,11 +31,18 @@ get_name(::LinearScalarAdvectionEquation2D) = "LinearScalarAdvectionEquation2D"
 varnames_cons(::LinearScalarAdvectionEquation2D) = SVector("scalar")
 varnames_prim(::LinearScalarAdvectionEquation2D) = SVector("scalar")
 
+# Calculates translated coordinates `x` for a periodic domain
+function x_trans_periodic_2d(x, domain_length = SVector(2, 2), center = SVector(0, 0))
+  x_normalized = x .- center
+  x_shifted = x_normalized .% domain_length
+  x_offset = ((x_shifted .< -0.5*domain_length) - (x_shifted .> 0.5*domain_length)) .* domain_length
+  return center + x_shifted + x_offset
+end
 
 # Set initial conditions at physical location `x` for time `t`
 function initial_condition_gauss(x, t, equation::LinearScalarAdvectionEquation2D)
   # Store translated coordinate for easy use of exact solution
-  x_trans = x - equation.advectionvelocity * t
+  x_trans = x_trans_periodic_2d(x - equation.advectionvelocity * t)
 
   return @SVector [exp(-(x_trans[1]^2 + x_trans[2]^2))]
 end
@@ -62,9 +69,6 @@ function initial_condition_sin_sin(x, t, equation::LinearScalarAdvectionEquation
 end
 
 function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation2D)
-  # Store translated coordinate for easy use of exact solution
-  x_trans = x - equation.advectionvelocity * t
-
   return @SVector [2.0]
 end