RFC: Source terms via callback

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.5.29, 0.6"
+Trixi = "0.5.29, 0.6, 0.7"
 julia = "1.8"
 
 [preferences.OrdinaryDiffEq]

LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl

@@ -0,0 +1,84 @@
+using LibTrixi
+using Trixi
+using OrdinaryDiffEq
+using Libdl
+
+# TODO: hard-coded path for SO
+so_handle = dlopen("./lib/libsource_terms.so")
+println("Opened library ", dlpath(so_handle))
+source_term_fptr = dlsym(so_handle, "source_term_wave")
+println("Obtained function pointer ", source_term_fptr)
+
+# TODO: global buffer to avoid allocation of temporary storage in source_term_callback
+sources_tmp::Vector{Cdouble} = Vector{Cdouble}(undef, 4)
+
+function source_term_callback(u, x, t, equations::CompressibleEulerEquations2D)
+    @ccall $source_term_fptr(u::Ptr{Cdouble}, x::Ptr{Cdouble}, t::Cdouble,
+                                      equations.gamma::Cdouble,
+                                      sources_tmp::Ptr{Cdouble})::Cvoid
+    return SVector(sources_tmp[1], sources_tmp[2], sources_tmp[3], sources_tmp[4])
+end
+
+# The function to create the simulation state needs to be named `init_simstate`
+function init_simstate()
+
+    ###############################################################################
+    # semidiscretization of the compressible Euler equations
+
+    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_term_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

LibTrixi.jl/test/Project.toml

@@ -5,7 +5,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 
 [compat]
 OrdinaryDiffEq = "6.53.2"
-Trixi = "0.5.29, 0.6"
+Trixi = "0.5.29, 0.6, 0.7"
 
 [preferences.OrdinaryDiffEq]
 PrecompileAutoSpecialize = false

examples/CMakeLists.txt

@@ -11,6 +11,14 @@ if ( NOT T8CODE_FOUND )
     list( FILTER EXAMPLES EXCLUDE REGEX ".*(t|T)8(c|C)(o|O)(d|D)(e|E).*" )
 endif()
 
+# example standalone library with callback
+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);

src/api.f90

@@ -297,11 +297,14 @@ integer(c_int) function trixi_ndofs_global(handle) bind(c)
       integer(c_int), value, intent(in) :: handle
     end function
 
-    !! @anchor trixi_ndofs_element_api_c
+    !>
+    !! @fn LibTrixi::trixi_ndofs_element::trixi_ndofs_element(handle)
     !!
     !! @brief Return number of degrees of freedom per element (cell).
     !!
     !! @param[in]  handle  simulation handle
+    !!
+    !! @see @ref trixi_ndofs_element_api_c "trixi_ndofs_element (C API)"
     integer(c_int) function trixi_ndofs_element(handle) bind(c)
       use, intrinsic :: iso_c_binding, only: c_int
       integer(c_int), value, intent(in) :: handle