RFC: Source terms via callback

.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

LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl

@@ -0,0 +1,94 @@
+using LibTrixi
+using Trixi
+using OrdinaryDiffEq
+using Libdl
+
+struct SourceTermsCallback
+    source_term_fptr::Ptr{Nothing}
+    buffer::Vector{Cdouble}
+
+    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(source_term_fptr, Vector{Cdouble}(undef, 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[1], buffer[2], buffer[3], buffer[4])
+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

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 )

examples/source_terms.c

@@ -0,0 +1,27 @@
+#include <stdio.h>
+#include <math.h>
+
+/* 
+ * 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);
+}

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);