Add an MPI capable test case (#94)

* run the warm bubble test using a P4estMesh

and compare to Trixi's TreeMesh version

* switch to TrixiAtmo elixir

* need to use Trixi initial condition

* tolerances

* add MPI test

* !fixup

* tolerances again

* switch to recommended mpiexec syntax

* add some safety measures from Trixi.jl

* tolerance

* fixed expected errors

* remove surface_integral

* remove dry bubble for serial tests

* Revert "remove dry bubble for serial tests"

This reverts commit fba9b73.

* re-add test, but without polydeg override

* check tighter tolerances

* revert tolerance

---------

Co-authored-by: Andrés Rueda-Ramírez <[email protected]>

Author: benegee

examples/elixir_moist_euler_dry_bubble.jl

@@ -1,4 +1,4 @@
-using OrdinaryDiffEq
+using OrdinaryDiffEqLowStorageRK
 using Trixi
 using TrixiAtmo
 using TrixiAtmo: source_terms_geopotential, cons2drypot
@@ -48,10 +48,8 @@ end
 
 initial_condition = initial_condition_warm_bubble
 
-boundary_condition = (x_neg = boundary_condition_periodic,
-                      x_pos = boundary_condition_periodic,
-                      y_neg = boundary_condition_slip_wall,
-                      y_pos = boundary_condition_slip_wall)
+boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall,
+                           :y_pos => boundary_condition_slip_wall)
 
 # Gravity source since Q_ph=0
 source_term = source_terms_geopotential
@@ -70,18 +68,18 @@ volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
 # and the EC flux (chandrashekar) as volume flux
 solver = DGSEM(basis, surface_flux, volume_integral)
 
-coordinates_min = (0.0, 0.0)
-coordinates_max = (20000.0, 10000.0)
+coordinates_min = (0.0, -5000.0)
+coordinates_max = (20000.0, 15000.0)
 
-cells_per_dimension = (64, 32)
-
-# Create curved mesh with 64 x 32 elements
-mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
-                      periodicity = (true, false))
+trees_per_dimension = (2, 2)
+mesh = P4estMesh(trees_per_dimension, polydeg = 1,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 initial_refinement_level = 5,
+                 periodicity = (true, false))
 
 # A semidiscretization collects data structures and functions for the spatial discretization
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
-                                    boundary_conditions = boundary_condition,
+                                    boundary_conditions = boundary_conditions,
                                     source_terms = source_term)
 
 ###############################################################################

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d"
@@ -8,6 +9,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 TrixiTest = "0a316866-cbd0-4425-8bcb-08103b2c1f26"
 
 [compat]
+MPI = "0.20.6"
 NLsolve = "4.5.1"
 OrdinaryDiffEq = "6.91"
 OrdinaryDiffEqLowStorageRK = "1.2"

test/runtests.jl

@@ -1,4 +1,5 @@
 using Test
+using MPI: mpiexec
 
 # Trixi.jl runs tests in parallel with CI jobs setting the `TRIXI_TEST` environment
 # variable to determine the subset of tests to execute.
@@ -9,6 +10,26 @@ const TRIXI_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3)
 const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3)
 
 @time @testset verbose=true showtiming=true "TrixiAtmo.jl tests" begin
+    # This is placed first since tests error out otherwise if `TRIXI_TEST == "all"`,
+    # at least on some systems.
+    @time if TRIXI_TEST == "all" || TRIXI_TEST == "mpi"
+        # Do a dummy `@test true`:
+        # If the process errors out the testset would error out as well,
+        # cf. https://github.com/JuliaParallel/MPI.jl/pull/391
+        @test true
+
+        run(`$(mpiexec()) -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --threads=1 --check-bounds=yes $(abspath("test_mpi.jl"))`)
+    end
+
+    @time if TRIXI_TEST == "all" || TRIXI_TEST == "threaded"
+        # Do a dummy `@test true`:
+        # If the process errors out the testset would error out as well,
+        # cf. https://github.com/JuliaParallel/MPI.jl/pull/391
+        @test true
+
+        run(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none $(abspath("test_threaded.jl"))`)
+    end
+
     @time if TRIXI_TEST == "all" || TRIXI_TEST == "trixi_consistency"
         include("test_trixi_consistency.jl")
     end
@@ -29,15 +50,6 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3)
         include("test_2d_shallow_water_covariant.jl")
     end
 
-    @time if TRIXI_TEST == "all" || TRIXI_TEST == "threaded"
-        # Do a dummy `@test true`:
-        # If the process errors out the testset would error out as well,
-        # cf. https://github.com/JuliaParallel/MPI.jl/pull/391
-        @test true
-
-        run(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none $(abspath("test_threaded.jl"))`)
-    end
-
     @time if TRIXI_TEST == "upstream"
         include("test_trixi_consistency.jl")
     end

test/test_2d_moist_euler.jl

@@ -11,22 +11,21 @@ EXAMPLES_DIR = TrixiAtmo.examples_dir()
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_moist_euler_dry_bubble.jl"),
                         l2=[
-                            1.300428671901329e-6,
-                            2.601090012108739e-5,
-                            0.0006660124630171347,
-                            0.008969786054960861,
+                            9.103834949215729e-7,
+                            1.8209333828866736e-5,
+                            0.0004709417153612775,
+                            0.006342004383628925,
                             0.0,
                             0.0
                         ],
                         linf=[
-                            1.0312042909910168e-5,
-                            0.00020625488871672815,
-                            0.006392107590872236,
-                            0.07612038028310053,
+                            1.0258282803210506e-5,
+                            0.0002051932980897675,
+                            0.006394867661494521,
+                            0.076401537633501,
                             0.0,
                             0.0
                         ],
-                        polydeg=3,
                         tspan=(0.0, 0.1))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)

test/test_mpi.jl

@@ -0,0 +1,58 @@
+module TestMPI
+
+using Test
+using TrixiAtmo
+
+include("test_trixiatmo.jl")
+
+# Start with a clean environment: remove Trixi.jl output directory if it exists
+outdir = "out"
+Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive = true)
+Trixi.MPI.Barrier(Trixi.mpi_comm())
+
+# CI with MPI and some tests fails often on Windows. Thus, we check whether this
+# is the case here. We use GitHub Actions, so we can check whether we run CI
+# in the cloud with Windows as follows, see also
+# https://docs.github.com/en/actions/learn-github-actions/environment-variables
+CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows()
+
+EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples")
+
+@testset "MPI tests" begin
+#! format: noindent
+
+@trixiatmo_testset "elixir_moist_euler_moist_bubble" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_moist_euler_dry_bubble.jl"),
+                        l2=[
+                            9.103834949215729e-7,
+                            1.8209333828866736e-5,
+                            0.0004709417153612775,
+                            0.006342004383628925,
+                            0.0,
+                            0.0
+                        ],
+                        linf=[
+                            1.0258282803210506e-5,
+                            0.0002051932980897675,
+                            0.006394867661494521,
+                            0.076401537633501,
+                            0.0,
+                            0.0
+                        ],
+                        tspan=(0.0, 0.1))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+end
+
+# Clean up afterwards: delete Trixi.jl output directory
+Trixi.mpi_isroot() && @test_nowarn rm(outdir, recursive = true)
+Trixi.MPI.Barrier(Trixi.mpi_comm())
+
+end # module

test/test_trixi_consistency.jl

@@ -11,13 +11,15 @@ isdir(outdir) && rm(outdir, recursive = true)
 @trixiatmo_testset "Dry air consistency check" begin
     # Dry air warm bubble test case in Trixi.jl
     maxiters = 100
-    trixi_elixir = joinpath(Trixi.examples_dir(), "structured_2d_dgsem",
+    trixi_elixir = joinpath(Trixi.examples_dir(), "tree_2d_dgsem",
                             "elixir_euler_warm_bubble.jl")
 
-    # Override maxiter and fluxes
+    # Override fluxes, polydeg, cfl, maxiters
     trixi_include(trixi_elixir;
                   volume_flux = Trixi.flux_chandrashekar,
                   surface_flux = Trixi.FluxLMARS(360.0),
+                  polydeg = 4,
+                  stepsize_callback = Trixi.StepsizeCallback(cfl = 0.2),
                   maxiters = maxiters)
 
     # Save errors
@@ -27,37 +29,30 @@ isdir(outdir) && rm(outdir, recursive = true)
     equations_trixi = Trixi.CompressibleEulerEquations2D(Main.warm_bubble_setup.gamma)
     add_zeros = SVector(zero(eltype(Main.sol)), zero(eltype(Main.sol)))
 
-    # Now use moist equations instead
-    equations_moist = CompressibleMoistEulerEquations2D()
-
-    # Redefine source terms for CompressibleMoistEulerEquations2D
-    @inline function (setup::Main.WarmBubbleSetup)(u, x, t,
-                                                   equations_moist::CompressibleMoistEulerEquations2D)
-        ret_trixi = setup(u, x, t, equations_trixi)
-        return vcat(ret_trixi, add_zeros)
-    end
-
-    # Redefine initial condition for CompressibleMoistEulerEquations2D
+    # Redefine initial condition in Trixi.jl for CompressibleMoistEulerEquations2D
+    # Different formulae were used!
     @inline function (setup::Main.WarmBubbleSetup)(x, t,
-                                                   equations_moist::CompressibleMoistEulerEquations2D)
+                                                   ::CompressibleMoistEulerEquations2D)
         ret_trixi = setup(x, t, equations_trixi)
         return vcat(ret_trixi, add_zeros)
     end
 
-    # Run again with overrides
-    trixi_include(trixi_elixir,
-                  equations = equations_moist,
-                  volume_flux = flux_chandrashekar,
-                  surface_flux = FluxLMARS(360.0),
+    # Now use the elixir in TrixiAtmo
+    elixir_atmo = joinpath(TrixiAtmo.examples_dir(), "elixir_moist_euler_dry_bubble.jl")
+
+    # Override initial condition, maxiters
+    trixi_include(elixir_atmo,
+                  initial_condition = Main.warm_bubble_setup,
                   maxiters = maxiters)
 
+    # Save errors
     errors_atmo = Main.analysis_callback(Main.sol)
 
     for (error_trixi, error_atmo) in zip(errors_trixi.l2, errors_atmo.l2)
         @test isapprox(error_trixi, error_atmo, rtol = 1e-12)
     end
     for (error_trixi, error_atmo) in zip(errors_trixi.linf, errors_atmo.linf)
-        @test isapprox(error_trixi, error_atmo, rtol = 3e-12)
+        @test isapprox(error_trixi, error_atmo, rtol = 1.1e-10)
     end
 end
 end