Add cubed sphere test cases

.github/workflows/ci.yml

@@ -38,13 +38,13 @@ jobs:
           - '1.10'
           - '1.11'
         t8code_version:
-          - '3.0.0'
+          - '3.0.1'
         include:
           - os: ubuntu-latest
             test_type: coverage
             arch: x64
             julia_version: '1.11'
-            t8code_version: '3.0.0'
+            t8code_version: '3.0.1'
     env:
       # Necessary for HDF5 to play nice with Julia
       LD_PRELOAD: /lib/x86_64-linux-gnu/libcurl.so.4
@@ -170,9 +170,7 @@ jobs:
               --t8code-library ../t8code-local/prefix/lib/libt8.so \
               --julia-depot ~/.julia \
               --force
-          cp ../install/share/libtrixi/LibTrixi.jl/examples/libelixir_tree1d_dgsem_advection_basic.jl .
-          cp ../install/share/libtrixi/LibTrixi.jl/examples/libelixir_p4est2d_dgsem_euler_sedov.jl .
-          cp ../install/share/libtrixi/LibTrixi.jl/examples/libelixir_t8code_2d_dgsem_advection_amr.jl .
+          cp ../install/share/libtrixi/LibTrixi.jl/examples/libelixir_* .
 
       - name: Initialize project directory (test_type == 'package-compiler')
         if: ${{ matrix.test_type == 'package-compiler' }}
@@ -215,7 +213,7 @@ jobs:
           ./build.sh
           mpirun -n 2 ./build/trixi_controller_simple_c \
               ../../libtrixi-julia \
-              ../../LibTrixi.jl/examples/libelixir_p4est2d_dgsem_euler_sedov.jl
+              ../../LibTrixi.jl/examples/libelixir_p4est2d_euler_sedov.jl
         env:
           LIBTRIXI_DEBUG: all
 
@@ -240,14 +238,19 @@ jobs:
         if: ${{ matrix.test_type == 'regular' || matrix.test_type == 'coverage' }}
         run: |
           cd libtrixi-julia
-          ../build/examples/trixi_controller_simple_c . libelixir_tree1d_dgsem_advection_basic.jl
-          ../build/examples/trixi_controller_simple_f . libelixir_tree1d_dgsem_advection_basic.jl
-          ../build/examples/trixi_controller_simple_c . libelixir_p4est2d_dgsem_euler_sedov.jl
-          ../build/examples/trixi_controller_simple_f . libelixir_p4est2d_dgsem_euler_sedov.jl
-          ../build/examples/trixi_controller_data_c . libelixir_t8code_2d_dgsem_advection_amr.jl
-          ../build/examples/trixi_controller_data_f . libelixir_t8code_2d_dgsem_advection_amr.jl
-          ../build/examples/trixi_controller_t8code_c . libelixir_t8code_2d_dgsem_advection_amr.jl
-          ../build/examples/trixi_controller_t8code_f . libelixir_t8code_2d_dgsem_advection_amr.jl
+          # all controllers
+          ../build/examples/trixi_controller_simple_c          . libelixir_tree1d_advection_basic.jl
+          ../build/examples/trixi_controller_simple_f          . libelixir_tree1d_advection_basic.jl
+          ../build/examples/trixi_controller_data_c            . libelixir_t8code2d_advection_amr.jl
+          ../build/examples/trixi_controller_data_f            . libelixir_t8code2d_advection_amr.jl
+          ../build/examples/trixi_controller_t8code_c          . libelixir_t8code2d_advection_amr.jl
+          ../build/examples/trixi_controller_t8code_f          . libelixir_t8code2d_advection_amr.jl
+          ../build/examples/trixi_controller_baroclinic_c      . libelixir_t8code3d_euler_baroclinic_instability.jl
+          ../build/examples/trixi_controller_baroclinic_f      . libelixir_t8code3d_euler_baroclinic_instability.jl
+          mpirun -n 2 ../build/examples/trixi_controller_mpi_c . libelixir_p4est2d_euler_sedov.jl
+          mpirun -n 2 ../build/examples/trixi_controller_mpi_f . libelixir_p4est2d_euler_sedov.jl
+          # remaining libelixirs
+          ../build/examples/trixi_controller_simple_c          . libelixir_t8code3d_euler_tracer.jl
         env:
           LIBTRIXI_DEBUG: all
 
@@ -257,7 +260,7 @@ jobs:
           cd build/examples
           mpirun -n 2 trixi_controller_simple_c \
               ../../libtrixi-julia \
-              ../../LibTrixi.jl/examples/libelixir_p4est2d_dgsem_euler_sedov.jl
+              ../../LibTrixi.jl/examples/libelixir_p4est2d_euler_sedov.jl
         env:
           LIBTRIXI_DEBUG: all
 
@@ -277,7 +280,15 @@ jobs:
                          "../build/examples/trixi_controller_t8code_c"   \
                          "../build/examples/trixi_controller_t8code_c ." \
                          "../build/examples/trixi_controller_t8code_f"   \
-                         "../build/examples/trixi_controller_t8code_f ."
+                         "../build/examples/trixi_controller_t8code_f ." \
+                         "../build/examples/trixi_controller_baroclinic_c"   \
+                         "../build/examples/trixi_controller_baroclinic_c ." \
+                         "../build/examples/trixi_controller_baroclinic_f"   \
+                         "../build/examples/trixi_controller_baroclinic_f ." \
+                         "../build/examples/trixi_controller_mpi_c"   \
+                         "../build/examples/trixi_controller_mpi_c ." \
+                         "../build/examples/trixi_controller_mpi_f"   \
+                         "../build/examples/trixi_controller_mpi_f ."
             do
               $command
               if [ $? -ne 2 ]; then

LibTrixi.jl/Project.toml

@@ -13,7 +13,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 MPI = "0.20.13"
 OrdinaryDiffEq = "6.53.2"
 Pkg = "1.8"
-Trixi = "0.7.16, 0.8, 0.9"
+Trixi = "0.9.12"
 julia = "1.8"
 
 [preferences.OrdinaryDiffEq]

LibTrixi.jl/README.md

@@ -18,8 +18,8 @@ run, and finalize a simulation by running the following code:
 ```julia
 julia> using LibTrixi
 
-julia> libelixir = pkgdir(LibTrixi, "examples", "libelixir_tree1d_dgsem_advection_basic.jl")
-"/path/to/libtrixi/LibTrixi.jl/examples/libelixir_tree1d_dgsem_advection_basic.jl"
+julia> libelixir = pkgdir(LibTrixi, "examples", "libelixir_tree1d_advection_basic.jl")
+"/path/to/libtrixi/LibTrixi.jl/examples/libelixir_tree1d_advection_basic.jl"
 
 julia> handle = trixi_initialize_simulation(libelixir); # initialize a new simulation setup
 

LibTrixi.jl/examples/libelixir_p4est3d_euler_baroclinic_instability.jl → LibTrixi.jl/examples/libelixir_t8code3d_euler_baroclinic_instability.jl

@@ -15,6 +15,43 @@ using LinearAlgebra
 using LibTrixi
 
 
+# Callable struct holding vectors with source terms
+struct SourceTerm
+    nnodesdim::Int
+    registry::LibTrixiDataRegistry
+end
+
+# We overwrite Trixi.jl's internal method here such that it calls source_terms with indices
+function Trixi.calc_sources!(du, u, t, source_terms::SourceTerm,
+                             equations::CompressibleEulerEquations3D, dg::DG, cache)
+    @unpack node_coordinates = cache.elements
+    Trixi.@threaded for element in eachelement(dg, cache)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_local = Trixi.get_node_vars(u, equations, dg, i, j, k, element)
+            du_local = source_terms(u_local, i, j, k, element, t, equations)
+            #x_local = Trixi.get_node_coords(node_coordinates, equations, dg,
+            #                                i, j, k, element)
+            #du_local_ref = source_terms_baroclinic_instability(u_local, x_local, t,
+            #                                                   equations)
+            Trixi.add_to_node_vars!(du, du_local, equations, dg, i, j, k, element)
+        end
+    end
+    return nothing
+end
+
+@inline function (source::SourceTerm)(u, i, j, k, element, t,
+                                      equations::CompressibleEulerEquations3D)
+    @unpack nnodesdim = source
+    index_global = (element-1) * nnodesdim^3 + (k-1) * nnodesdim^2 + (j-1) * nnodesdim + i
+    du2::Vector{Float64} = source.registry[1]
+    du3::Vector{Float64} = source.registry[2]
+    du4::Vector{Float64} = source.registry[3]
+    du5::Vector{Float64} = source.registry[4]
+    return SVector(zero(eltype(u)), du2[index_global], du3[index_global],
+                   du4[index_global], du5[index_global])
+end
+
+
 # Initial condition for an idealized baroclinic instability test
 # https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A
 function initial_condition_baroclinic_instability(x, t,
@@ -193,38 +230,51 @@ end
 # The function to create the simulation state needs to be named `init_simstate`
 function init_simstate()
 
-    ###############################################################################
-    # Setup for the baroclinic instability test
+    # compressible euler equations
     gamma = 1.4
     equations = CompressibleEulerEquations3D(gamma)
 
-    ###############################################################################
-    # semidiscretization of the problem
-
+    # setup of the problem
     initial_condition = initial_condition_baroclinic_instability
 
     boundary_conditions = Dict(:inside => boundary_condition_slip_wall,
                                :outside => boundary_condition_slip_wall)
 
-    # This is a good estimate for the speed of sound in this example.
-    # Other values between 300 and 400 should work as well.
+    # estimate for the speed of sound
     surface_flux = FluxLMARS(340)
     volume_flux = flux_kennedy_gruber
     solver = DGSEM(polydeg = 5, surface_flux = surface_flux,
                    volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
-    trees_per_cube_face = (16, 8)
-    mesh = Trixi.P4estMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000.0,
-                                      polydeg = 5, initial_refinement_level = 0)
+    # for nice results, use 4 and 8 here
+    lat_lon_levels = 2
+    layers = 4
+    mesh = Trixi.T8codeMeshCubedSphere(lat_lon_levels, layers, 6.371229e6, 30000.0,
+                                       polydeg = 5, initial_refinement_level = 0)
+
+    # create the data registry and four vectors for the source terms
+    registry = LibTrixiDataRegistry(undef, 4)
+
+    nnodesdim = Trixi.nnodes(solver)
+    nnodes = nnodesdim^3
+    nelements = Trixi.ncells(mesh)
+
+    # provide some data because calc_sources! will already be called during initialization
+    # Note: the data pointers in the registry will be overwritten before the first real use
+    registry[1] = zeros(Float64, nelements*nnodes)
+    registry[2] = zeros(Float64, nelements*nnodes)
+    registry[3] = zeros(Float64, nelements*nnodes)
+    registry[4] = zeros(Float64, nelements*nnodes)
+
+    source_term_data_registry = SourceTerm(nnodesdim, registry)
 
     semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
-                                        source_terms = source_terms_baroclinic_instability,
+                                        source_terms = source_term_data_registry,
                                         boundary_conditions = boundary_conditions)
 
-    ###############################################################################
-    # ODE solvers, callbacks etc.
-
-    tspan = (0.0, 12 * 24 * 60 * 60.0) # time in seconds for 12 days#
+    # for nice results, use 10 days
+    days = 0.02
+    tspan = (0.0, days * 24 * 60 * 60.0)
 
     ode = semidiscretize(semi, tspan)
 
@@ -235,33 +285,24 @@ function init_simstate()
 
     alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-    save_solution = SaveSolutionCallback(interval = 5000,
+    save_solution = SaveSolutionCallback(interval = 500,
                                          save_initial_solution = true,
                                          save_final_solution = true,
                                          solution_variables = cons2prim,
-                                         output_directory = "out_baroclinic",)
+                                         output_directory = "out_baroclinic")
 
     callbacks = CallbackSet(summary_callback,
                             analysis_callback,
                             alive_callback,
                             save_solution)
 
-
-    ###############################################################################
-    # create the time integrator
-
-    # OrdinaryDiffEq's `integrator`
-    # Use a Runge-Kutta method with automatic (error based) time step size control
-    integrator = init(ode, RDPK3SpFSAL49();
+    # use a Runge-Kutta method with automatic (error based) time step size control
+    integrator = init(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.False());
                       abstol = 1.0e-6, reltol = 1.0e-6,
-                      ode_default_options()...,
-                      callback = callbacks,
-                      maxiters=5e5);
-
-    ###############################################################################
-    # Create simulation state
+                      ode_default_options()..., callback = callbacks, maxiters=1e7);
 
-    simstate = SimulationState(semi, integrator)
+    # create simulation state
+    simstate = SimulationState(semi, integrator, registry)
 
     return simstate
 end

LibTrixi.jl/examples/libelixir_t8code3d_euler_tracer.jl

@@ -0,0 +1,124 @@
+# A manufactured solution of a circular wind with constant angular velocity
+# on a planetary-sized cubed sphere mesh with a blob detected by AMR
+#
+# Note that this libelixir is based on an elixir by Erik Faulhaber for Trixi.jl
+# Source: https://github.com/trixi-framework/Trixi.jl/blob/main/examples/p4est_3d_dgsem/elixir_euler_circular_wind_nonconforming.jl
+
+using OrdinaryDiffEq
+using Trixi
+using LinearAlgebra
+using LibTrixi
+
+
+function initial_condition_circular_wind(x, t, equations::CompressibleEulerEquations3D)
+    radius_earth = 6.371229e6
+    lambda, phi, r = cart_to_sphere(x)
+
+    p = 1e5
+    v1 = -10 * x[2] / radius_earth
+    v2 = 10 * x[1] / radius_earth
+    v3 = 0.0
+    rho = 1.0 + 0.1  * exp(-50 * ((lambda-1.0)^2/2.0 + (phi-0.4)^2)) +
+                0.08 * exp(-100 * ((lambda-0.8)^2/4.0 + (phi-0.5)^2))
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+@inline function source_terms_circular_wind(u, x, t,
+                                            equations::CompressibleEulerEquations3D)
+    radius_earth = 6.371229e6
+    rho = u[1]
+
+    du1 = 0.0
+    du2 = -rho * (10 / radius_earth) * (10 * x[1] / radius_earth)
+    du3 = -rho * (10 / radius_earth) * (10 * x[2] / radius_earth)
+    du4 = 0.0
+    du5 = 0.0
+
+    return SVector(du1, du2, du3, du4, du5)
+end
+
+function cart_to_sphere(x)
+    r = norm(x)
+    lambda = atan(x[2], x[1])
+    if lambda < 0
+        lambda += 2 * pi
+    end
+    phi = asin(x[3] / r)
+
+    return lambda, phi, r
+end
+
+
+# The function to create the simulation state needs to be named `init_simstate`
+function init_simstate()
+
+    # compressible Euler equations
+    gamma = 1.4
+    equations = CompressibleEulerEquations3D(gamma)
+
+    # setup of the problem
+    initial_condition = initial_condition_circular_wind
+
+    boundary_conditions = Dict(:inside => boundary_condition_slip_wall,
+                               :outside => boundary_condition_slip_wall)
+
+    # estimate for speed of sound
+    surface_flux = FluxLMARS(374)
+    solver = DGSEM(polydeg = 3, surface_flux = surface_flux)
+
+    # increase trees_per_cube_face to 4 to get nicer results
+    lat_lon_levels = 2
+    layers = 1
+    mesh = Trixi.T8codeMeshCubedSphere(lat_lon_levels, layers, 6.371229e6, 30000.0,
+                                       polydeg = 3, initial_refinement_level = 0)
+
+    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                        source_terms = source_terms_circular_wind,
+                                        boundary_conditions = boundary_conditions)
+
+    # increase number of days
+    days = 0.1
+    tspan = (0.0, days * 24 * 60 * 60.0)
+
+    ode = semidiscretize(semi, tspan)
+
+    summary_callback = SummaryCallback()
+
+    analysis_interval = 5000
+    analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+    alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+    save_solution = SaveSolutionCallback(interval = 2000,
+                                         save_initial_solution = true,
+                                         save_final_solution = true,
+                                         solution_variables = cons2prim,
+                                         output_directory = "out_tracer")
+
+    amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                          base_level = 0,
+                                          med_level = 1, med_threshold = 1.004,
+                                          max_level = 3, max_threshold = 1.11)
+
+    amr_callback = AMRCallback(semi, amr_controller,
+                               interval = 2000,
+                               adapt_initial_condition = true,
+                               adapt_initial_condition_only_refine = true)
+
+    callbacks = CallbackSet(summary_callback,
+                            analysis_callback,
+                            alive_callback,
+                            amr_callback,
+                            save_solution)
+
+    # use a Runge-Kutta method with automatic (error based) time step size control
+    integrator = init(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.False());
+                      abstol = 1.0e-6, reltol = 1.0e-6,
+                      ode_default_options()..., callback = callbacks, maxiters=1e7);
+
+    # create simulation state
+    simstate = SimulationState(semi, integrator)
+
+    return simstate
+end