diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b8cfac1a..948d194b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -165,6 +165,7 @@ jobs: mkdir libtrixi-julia cd libtrixi-julia cp ../install/share/libtrixi/LibTrixi.jl/examples/libelixir_* . + cp ../install/lib/libsource_terms.so . mkdir julia-depot ../install/bin/libtrixi-init-julia .. \ --hdf5-library /usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so \ @@ -243,6 +244,7 @@ jobs: # 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_simple_c . libelixir_structured2d_source_terms_callback.jl ../build/examples/trixi_controller_data_load_c . libelixir_t8code2d_euler_tracer_amr.jl ../build/examples/trixi_controller_data_load_f . libelixir_t8code2d_euler_tracer_amr.jl ../build/examples/trixi_controller_data_store_c . libelixir_t8code2d_euler_tracer_amr.jl diff --git a/LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl b/LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl new file mode 100644 index 00000000..05437e7b --- /dev/null +++ b/LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl @@ -0,0 +1,94 @@ +using LibTrixi +using Trixi +using OrdinaryDiffEq +using Libdl + +struct SourceTermsCallback{NVARS} + source_term_fptr::Ptr{Nothing} + buffer::Base.RefValue{NTuple{NVARS, Float64}} + + function SourceTermsCallback(; n_vars, so_path) + so_handle = dlopen(so_path) + @info "Opened library ", dlpath(so_handle) + source_term_fptr = dlsym(so_handle, "source_term_wave") + @info "Obtained function pointer ", source_term_fptr + + new{n_vars}(source_term_fptr, Ref(ntuple(i -> 0.0, n_vars))) + end +end + +function (callback::SourceTermsCallback)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack source_term_fptr, buffer = callback + @ccall $source_term_fptr(u::Ptr{Cdouble}, x::Ptr{Cdouble}, t::Cdouble, + equations.gamma::Cdouble, buffer::Ptr{Cdouble})::Cvoid + return SVector(buffer[]) +end + +# The function to create the simulation state needs to be named `init_simstate` +function init_simstate() + + ############################################################################### + # semidiscretization of the compressible Euler equations + + # use external function in locally compiled shared object + source_terms_callback = SourceTermsCallback(n_vars = 4, + so_path = joinpath(@__DIR__, + "libsource_terms.so")) + + equations = CompressibleEulerEquations2D(1.4) + + initial_condition = initial_condition_convergence_test + + solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + + coordinates_min = (0.0, 0.0) + coordinates_max = (2.0, 2.0) + + cells_per_dimension = (16, 16) + + mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_callback) + + ############################################################################### + # ODE solvers, callbacks etc. + + tspan = (0.0, 2.0) + ode = semidiscretize(semi, tspan) + + summary_callback = SummaryCallback() + + analysis_interval = 100 + analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + + alive_callback = AliveCallback(analysis_interval = analysis_interval) + + save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + + stepsize_callback = StepsizeCallback(cfl = 1.0) + + callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + + ############################################################################### + # create OrdinaryDiffEq's `integrator` + + integrator = init(ode, + CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # will be overwritten by the stepsize_callback + save_everystep=false, + callback=callbacks); + + ############################################################################### + # Create simulation state + + simstate = SimulationState(semi, integrator) + + return simstate +end diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 508fb4bf..6910f2ec 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -17,6 +17,14 @@ if ( NOT T8CODE_FOUND ) list( FILTER EXAMPLES EXCLUDE REGEX "trixi_controller_baroclinic.*" ) endif() +# exemplary external library providing function used to compute source terms +add_library ( source_terms SHARED + source_terms.h + source_terms.c +) +install( TARGETS source_terms ) + + foreach ( EXAMPLE ${EXAMPLES} ) get_filename_component ( EXAMPLE_EXT ${EXAMPLE} EXT ) diff --git a/examples/source_terms.c b/examples/source_terms.c new file mode 100644 index 00000000..3357b758 --- /dev/null +++ b/examples/source_terms.c @@ -0,0 +1,27 @@ +#include +#include + +/* + * Example for an external source term evaluation + * Diagonally running wave, taken from Trixi.jl's source_terms_convergence_test + */ +void source_term_wave(const double * u, const double * x, const double t, + const double gamma, double * sourceterm) { + + const double c = 2.0; + const double A = 0.1; + const double L = 2.0; + const double f = 1.0 / L; + const double omega = 2 * M_PI * f; + + const double si = sin(omega * (x[0] + x[1] - t)); + const double co = cos(omega * (x[0] + x[1] - t)); + const double rho = c + A * si; + const double rho_x = omega * A * co; + const double tmp = (2 * rho - 1) * (gamma - 1); + + sourceterm[0] = rho_x; + sourceterm[1] = rho_x * (1 + tmp); + sourceterm[2] = sourceterm[1]; + sourceterm[3] = 2 * rho_x * (rho + tmp); +} diff --git a/examples/source_terms.h b/examples/source_terms.h new file mode 100644 index 00000000..cdddc196 --- /dev/null +++ b/examples/source_terms.h @@ -0,0 +1,2 @@ +void source_term_wave(const double * u, const double * x, const double t, + const double gamma, double * sourceterm);