RFC: source terms via global vector storage (#173)

* 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

* 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 get_t8code_forest to Fortran interface

* format

* add trixi_store_in_database to Fortran interface

* 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

* rename get_time to get_simulation_time

* rename database to registry

* spurious Ref access

---------

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
Co-authored-by: Benedict <135045760+bgeihe@users.noreply.github.com>

Author: benegee

LibTrixi.jl/examples/libelixir_p4est2d_dgsem_euler_sedov.jl

@@ -106,7 +106,9 @@ function init_simstate()
     ###############################################################################
     # Create simulation state
 
-    simstate = SimulationState(semi, integrator)
+    # registry only used for tests
+    registry = LibTrixiDataRegistry(undef, 1)
+    simstate = SimulationState(semi, integrator, registry)
 
     return simstate
 end

LibTrixi.jl/src/LibTrixi.jl

@@ -58,6 +58,9 @@ export trixi_load_primitive_vars,
 export trixi_load_element_averaged_primitive_vars,
        trixi_load_element_averaged_primitive_vars_cfptr,
        trixi_load_element_averaged_primitive_vars_jl
+export trixi_register_data,
+       trixi_register_data_cfptr,
+       trixi_register_data_jl
 export trixi_version_library,
        trixi_version_library_cfptr,
        trixi_version_library_jl
@@ -82,8 +85,13 @@ export trixi_get_t8code_forest,
 export trixi_eval_julia,
        trixi_eval_julia_cfptr,
        trixi_eval_julia_jl
+export trixi_get_simulation_time,
+       trixi_get_simulation_time_cfptr,
+       trixi_get_simulation_time_jl
+
 
 export SimulationState, store_simstate, load_simstate, delete_simstate!
+export LibTrixiDataRegistry
 
 
 # global storage of name and version information of loaded packages

LibTrixi.jl/src/api_c.jl

@@ -461,6 +461,53 @@ trixi_load_primitive_vars_cfptr() =
     @cfunction(trixi_load_primitive_vars, Cvoid, (Cint, Cint, Ptr{Cdouble}))
 
 
+"""
+    trixi_register_data(data::Ptr{Cdouble}, size::Cint, index::Cint,
+                        simstate_handle::Cint)::Cvoid
+
+Store data vector in current simulation's registry.
+
+A reference to the passed data array `data` will be stored in the registry of the simulation
+given by `simstate_handle` at given `index`. The registry object has to be created in
+`init_simstate()` of the running libelixir and can be used throughout the simulation.
+
+The registry object has to exist, has to be of type `LibTrixiDataRegistry`, and has to hold
+enough data references such that access at `index` is valid.
+Memory storage remains on the user side. It must not be deallocated as long as it might be
+accessed via the registry. The size of `data` has to match `size`.
+"""
+function trixi_register_data end
+
+Base.@ccallable function trixi_register_data(simstate_handle::Cint, index::Cint,
+                                             size::Cint, data::Ptr{Cdouble})::Cvoid
+    simstate = load_simstate(simstate_handle)
+
+    # convert C to Julia array
+    data_jl = unsafe_wrap(Array, data, size)
+
+    trixi_register_data_jl(simstate, index, data_jl)
+    return nothing
+end
+
+trixi_register_data_cfptr() =
+    @cfunction(trixi_register_data, Cvoid, (Cint, Cint, Cint, Ptr{Cdouble},))
+
+
+"""
+    trixi_get_simulation_time(simstate_handle::Cint)::Cdouble
+
+Return current physical time.
+"""
+function trixi_get_simulation_time end
+
+Base.@ccallable function trixi_get_simulation_time(simstate_handle::Cint)::Cdouble
+    simstate = load_simstate(simstate_handle)
+    return trixi_get_simulation_time_jl(simstate)
+end
+
+trixi_get_simulation_time_cfptr() = @cfunction(trixi_get_simulation_time, Cdouble, (Cint,))
+
+
 """
     trixi_load_element_averaged_primitive_vars(simstate_handle::Cint, variable_id::Cint,
                                             data::Ptr{Cdouble})::Cvoid

LibTrixi.jl/src/api_jl.jl

@@ -199,6 +199,20 @@ function trixi_load_element_averaged_primitive_vars_jl(simstate, variable_id, da
 end
 
 
+function trixi_register_data_jl(simstate, index, data)
+    simstate.registry[index] = data
+    if show_debug_output()
+        println("New data vector registered at index ", index)
+    end
+    return nothing
+end
+
+
+function trixi_get_simulation_time_jl(simstate)
+    return simstate.integrator.t
+end
+
+
 function trixi_get_t8code_forest_jl(simstate)
     mesh, _, _, _ = mesh_equations_solver_cache(simstate.semi)
     return mesh.forest

LibTrixi.jl/src/simulationstate.jl

@@ -1,12 +1,21 @@
+const LibTrixiDataRegistry = Vector{Vector{Float64}}
+
 """
     SimulationState
 
-Data structure to store a simulation state that consists of a semidiscretization
-plus the time integrator.
+Data structure to store a simulation state consisting of
+- a semidiscretization
+- the time integrator
+- an optional array of data vectors
 """
 mutable struct SimulationState{SemiType, IntegratorType}
     semi::SemiType
     integrator::IntegratorType
+    registry::LibTrixiDataRegistry
+
+    function SimulationState(semi, integrator, registry = LibTrixiDataRegistry())
+        return new{typeof(semi), typeof(integrator)}(semi, integrator, registry)
+    end
 end
 
 # Global variables to store different simulation states

LibTrixi.jl/test/test_interface.jl

@@ -59,9 +59,25 @@ end
     dt_jl = trixi_calculate_dt_jl(simstate_jl)
     @test dt_c == dt_jl
 
+    # compare time
+    time_c = trixi_get_simulation_time(handle)
+    time_jl = trixi_get_simulation_time_jl(simstate_jl)
+    @test time_c == time_jl
+
     # compare finished status
     @test trixi_is_finished(handle) == 0
     @test !trixi_is_finished_jl(simstate_jl)
+
+    # manually increase registries (for testing only!)
+    push!(simstate_jl.registry, Vector{Float64}())
+    push!(LibTrixi.simstates[handle].registry, Vector{Float64}())
+    # store a vector
+    test_data = [1.0, 2.0, 3.0]
+    trixi_register_data(handle, Int32(1), Int32(3), pointer(test_data))
+    trixi_register_data_jl(simstate_jl, 1, test_data)
+    # check that the same memory is referenced
+    @test pointer(simstate_jl.registry[1]) ==
+        pointer(LibTrixi.simstates[handle].registry[1])
 end
 
 

examples/external/CMakeLists.txt

@@ -1,5 +1,5 @@
 # Specify the minimum version
-cmake_minimum_required ( VERSION 3.9 )
+cmake_minimum_required ( VERSION 3.12 )
 
 # Specify a project name
 project(ExternalLibTrixi)

src/api.c

@@ -24,6 +24,7 @@ enum {
     TRIXI_FPTR_LOAD_NODE_WEIGHTS,
     TRIXI_FTPR_LOAD_PRIMITIVE_VARS,
     TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VARS,
+    TRIXI_FTPR_REGISTER_DATA,
     TRIXI_FTPR_VERSION_LIBRARY,
     TRIXI_FTPR_VERSION_LIBRARY_MAJOR,
     TRIXI_FTPR_VERSION_LIBRARY_MINOR,
@@ -32,6 +33,7 @@ enum {
     TRIXI_FTPR_VERSION_JULIA_EXTENDED,
     TRIXI_FTPR_EVAL_JULIA,
     TRIXI_FTPR_GET_T8CODE_FOREST,
+    TRIXI_FPTR_GET_SIMULATION_TIME,
 
     // The last one is for the array size
     TRIXI_NUM_FPTRS
@@ -60,14 +62,16 @@ static const char* trixi_function_pointer_names[] = {
     [TRIXI_FPTR_LOAD_NODE_WEIGHTS]                    = "trixi_load_node_weights_cfptr",
     [TRIXI_FTPR_LOAD_PRIMITIVE_VARS]                  = "trixi_load_primitive_vars_cfptr",
     [TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VARS] = "trixi_load_element_averaged_primitive_vars_cfptr",
+    [TRIXI_FTPR_REGISTER_DATA]                        = "trixi_register_data_cfptr",
     [TRIXI_FTPR_VERSION_LIBRARY]                      = "trixi_version_library_cfptr",
     [TRIXI_FTPR_VERSION_LIBRARY_MAJOR]                = "trixi_version_library_major_cfptr",
     [TRIXI_FTPR_VERSION_LIBRARY_MINOR]                = "trixi_version_library_minor_cfptr",
     [TRIXI_FTPR_VERSION_LIBRARY_PATCH]                = "trixi_version_library_patch_cfptr",
     [TRIXI_FTPR_VERSION_JULIA]                        = "trixi_version_julia_cfptr",
     [TRIXI_FTPR_VERSION_JULIA_EXTENDED]               = "trixi_version_julia_extended_cfptr",
     [TRIXI_FTPR_EVAL_JULIA]                           = "trixi_eval_julia_cfptr",
-    [TRIXI_FTPR_GET_T8CODE_FOREST]                    = "trixi_get_t8code_forest_cfptr"
+    [TRIXI_FTPR_GET_T8CODE_FOREST]                    = "trixi_get_t8code_forest_cfptr",
+    [TRIXI_FPTR_GET_SIMULATION_TIME]                  = "trixi_get_simulation_time_cfptr"
 };
 
 // Track initialization/finalization status to prevent unhelpful errors
@@ -691,6 +695,58 @@ void trixi_load_element_averaged_primitive_vars(int handle, int variable_id, dou
 }
 
 
+/**
+ * @anchor trixi_register_data_api_c
+ *
+ * @brief Store data vector in current simulation's registry
+ *
+ * A reference to the passed data array `data` will be stored in the registry of the
+ * simulation given by `simstate_handle` at given `index`. The registry object has to be
+ * created in `init_simstate()` of the running libelixir and can be used throughout the
+ * simulation.
+ *
+ * The registry object has to exist, has to be of type `LibTrixiDataRegistry`, and has to
+ * hold enough data references such that access at `index` is valid.
+ * Memory storage remains on the user side. It must not be deallocated as long as it might
+ * be accessed via the registry. The size of `data` has to match `size`.
+ *
+ * @param[in]  handle  simulation handle
+ * @param[in]  index   index in registry where data vector will be stored
+ * @param[in]  size    size of given data vector
+ * @param[in]  data    data vector to store
+ */
+void trixi_register_data(int handle, int index, int size, const double * data) {
+
+    // Get function pointer
+    void (*register_data)(int, int, int, const double *) =
+        trixi_function_pointers[TRIXI_FTPR_REGISTER_DATA];
+
+    // Call function
+    register_data(handle, index, size, data);
+}
+
+
+/**
+ * @anchor trixi_get_simulation_time_api_c
+ *
+ * @brief Return current physical time.
+ *
+ * @param[in]  handle  simulation handle
+ *
+ * @return physical time
+ */
+double trixi_get_simulation_time(int handle) {
+
+    // Get function pointer
+    double (*get_simulation_time)(int) =
+        trixi_function_pointers[TRIXI_FPTR_GET_SIMULATION_TIME];
+
+    // Call function
+    return get_simulation_time(handle);
+}
+
+
+
 /******************************************************************************************/
 /* T8code                                                                                 */
 /******************************************************************************************/

src/api.f90

@@ -388,9 +388,24 @@ subroutine trixi_load_primitive_vars(handle, variable_id, data) bind(c)
       use, intrinsic :: iso_c_binding, only: c_int, c_double
       integer(c_int), value, intent(in) :: handle
       integer(c_int), value, intent(in) :: variable_id
-      real(c_double), dimension(*), intent(in) :: data
+      real(c_double), dimension(*), intent(out) :: data
     end subroutine
 
+    !>
+    !! @fn LibTrixi::trixi_get_simulation_time::trixi_get_simulation_time(handle)
+    !!
+    !! @brief Return current physical time.
+    !!
+    !! @param[in]  handle  simulation handle
+    !!
+    !! @return  physical time
+    !!
+    !! @see @ref trixi_get_simulation_time_api_c "trixi_get_simulation_time (C API)"
+    real(c_double) function trixi_get_simulation_time(handle) bind(c)
+      use, intrinsic :: iso_c_binding, only: c_int, c_double
+      integer(c_int), value, intent(in) :: handle
+    end function
+
     !>
     !! @fn LibTrixi::trixi_load_element_averaged_primitive_vars::trixi_load_element_averaged_primitive_vars(handle, variable_id, data)
     !!
@@ -405,6 +420,25 @@ subroutine trixi_load_element_averaged_primitive_vars(handle, variable_id, data)
       use, intrinsic :: iso_c_binding, only: c_int, c_double
       integer(c_int), value, intent(in) :: handle
       integer(c_int), value, intent(in) :: variable_id
+      real(c_double), dimension(*), intent(out) :: data
+    end subroutine
+
+    !>
+    !! @fn LibTrixi::trixi_register_data::trixi_register_data(handle, variable_id, data)
+    !!
+    !! @brief Store data vector in current simulation's registry
+    !!
+    !! @param[in]  handle  simulation handle
+    !! @param[in]  index   index in registry where data vector will be stored
+    !! @param[in]  size    size of given data vector
+    !! @param[in]  data    data vector to store
+    !!
+    !! @see @ref trixi_register_data_api_c "trixi_register_data (C API)"
+    subroutine trixi_register_data(handle, index, size, data) bind(c)
+      use, intrinsic :: iso_c_binding, only: c_int, c_double
+      integer(c_int), value, intent(in) :: handle
+      integer(c_int), value, intent(in) :: index
+      integer(c_int), value, intent(in) :: size
       real(c_double), dimension(*), intent(in) :: data
     end subroutine
 

src/trixi.h

@@ -34,10 +34,12 @@ int trixi_ndofselement(int handle);
 int trixi_nvariables(int handle);
 int trixi_nnodes(int handle);
 double trixi_calculate_dt(int handle);
+double trixi_get_simulation_time(int handle);
 void trixi_load_node_reference_coordinates(int handle, double* node_coords);
 void trixi_load_node_weights(int handle, double* node_weights);
 void trixi_load_primitive_vars(int handle, int variable_id, double * data);
 void trixi_load_element_averaged_primitive_vars(int handle, int variable_id, double * data);
+void trixi_register_data(int handle, int index, int size, const double * data);
 
 // T8code
 #if !defined(T8_H) && !defined(T8_FOREST_GENERAL_H)