Skip to content

Commit

Permalink
Add cubed sphere test cases (#199)
Browse files Browse the repository at this point in the history
* add getter for number of dofs

* change trixi_load_cell_average

add parameter index and only get averaged of the variable at position index

* fix: ensure Int32

* adapt reference value

* add getter for all dofs values

* adapt next value

* add trixi_store_in_database

proof of concept for creating data vectors in libelixirs which can be
filled later by external applications

* add trixi_ndofs_element

* add tests

* add missing parts in tests

* fix tests

* deallocate first

* add trixi_load_prim to Fortran API

* reference value

* get doxygen right

* add trixi_get_time, trixi_load_node_coordinates

required in this PR, subject to change

* implement source terms controller

demonstrates source terms via database

* libelixir with source terms via database

* change Array to Vector

* spelling

* libelixir for baroclinic instability

* remove method only required by steady state correction

* Update LibTrixi.jl/examples/libelixir_p4est3d_euler_baroclinic_instability.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* make everything more consistent!

* use const double *

* missed merge conflict

* transpose calloc args

* add functions to get quadrature information

* update CI badge URL

* remove load_node_coordinates

* add baroclinic instability elixir and controller

* fix

* clean up

* add get_t8code_forest to Fortran interface

* format

* add trixi_store_in_database to Fortran interface

* add fortran controller for baroclinic test case

* remove parameter t

was not used and cause either omitted parameter error or unused parameter warning using older GCCs

* spelling

* adapt min version according to main CMakeLists.txt

* missed while merging

* change default argument for DataBase parameter

* remove deprecated export

* test for get_time

* add tests for trixi_store_in_database

* need Int32

* cannot compare Refs

* check based on address

* relax error tolerance

* remove deprecated example

* new cubed sphere examples

* remove deprecated libelixir

* unify naming of libelixirs

* fixes

* adapt index.md to README.md

* change to non-adaptive time integrator

* missing controlers in ci checks

* use zeros in initial data

* missing renamings

* time not needed anymore

* tuned simulation parameters

* itree and ielement need to be 0-based

* one more to free

* ouch...

* less output

* typo in docstring

* adopt new t8code ForestWrapper

* switch to newest t8code release

* cubed sphere stuff only works with latest Trixi release

* missed renaming

* enable t8code for package compiler as well

* review

* filter out baroclinic example

---------

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
Co-authored-by: Benedict <[email protected]>
Co-authored-by: Benedict Geihe <[email protected]>
  • Loading branch information
4 people authored Jan 16, 2025
1 parent 5131626 commit 0e0ea3d
Show file tree
Hide file tree
Showing 23 changed files with 620 additions and 258 deletions.
43 changes: 27 additions & 16 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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' }}
Expand Down Expand Up @@ -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

Expand 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

Expand 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

Expand 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
Expand Down
2 changes: 1 addition & 1 deletion LibTrixi.jl/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
4 changes: 2 additions & 2 deletions LibTrixi.jl/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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)

Expand All @@ -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
124 changes: 124 additions & 0 deletions LibTrixi.jl/examples/libelixir_t8code3d_euler_tracer.jl
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 0e0ea3d

Please sign in to comment.