RFC: source terms via global vector storage

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)

test/c/interface_c.cpp

@@ -9,11 +9,6 @@ extern "C" {
 // Julia project path defined via cmake
 const char * julia_project_path = JULIA_PROJECT_PATH;
 
-// Example libexlixir
-const char * libelixir_path =
-  "../../../LibTrixi.jl/examples/libelixir_p4est2d_dgsem_euler_sedov.jl";
-
-
 TEST(CInterfaceTest, JuliaProject) {
 
     // be evil

test/c/simulation.cpp

@@ -40,6 +40,12 @@ TEST(CInterfaceTest, SimulationRun) {
     EXPECT_DEATH(trixi_is_finished(42),
                  "the provided handle was not found in the stored simulation states: 42");
 
+    // Store a vector in registry
+    std::vector<double> test_data(3);
+    trixi_register_data(handle, 1, 3, test_data.data());
+    EXPECT_DEATH(trixi_register_data(handle, 2, 3, test_data.data()),
+                 "BoundsError");
+
     // Do 10 simulation steps
     for (int i = 0; i < 10; ++i) {
         trixi_step(handle);
@@ -48,6 +54,10 @@ TEST(CInterfaceTest, SimulationRun) {
     // Check time step length
     double dt = trixi_calculate_dt(handle);
     EXPECT_NEAR(dt, 0.0028566952356658794, 1e-17);
+
+    // Check time
+    double time = trixi_get_simulation_time(handle);
+    EXPECT_NEAR(time, 0.0304927240859461, 1e-16);
 
     // Check finished status
     int finished_status = trixi_is_finished(handle);
@@ -95,7 +105,7 @@ TEST(CInterfaceTest, SimulationRun) {
     std::vector<double> energy(ndofs);
     trixi_load_primitive_vars(handle, 1, rho.data());
     trixi_load_primitive_vars(handle, 4, energy.data());
-    // check memory boarders
+    // check memory borders
     EXPECT_DOUBLE_EQ(rho[0],          1.0);
     EXPECT_DOUBLE_EQ(rho[ndofs-1],    1.0);
     EXPECT_DOUBLE_EQ(energy[0],       1.0e-5);
@@ -111,7 +121,7 @@ TEST(CInterfaceTest, SimulationRun) {
     trixi_load_element_averaged_primitive_vars(handle, 3, v2_averages.data());
     trixi_load_element_averaged_primitive_vars(handle, 4, e_averages.data());
     if (nranks == 1) {
-        // check memory boarders (densities at the beginning, energies at the end)
+        // check memory borders (densities at the beginning, energies at the end)
         EXPECT_DOUBLE_EQ(rho_averages[0],         1.0);
         EXPECT_DOUBLE_EQ(e_averages[nelements-1], 1.0e-5);
         // check values somewhere near the center (expect symmetries)
@@ -128,7 +138,7 @@ TEST(CInterfaceTest, SimulationRun) {
     }
     else if (nranks == 2) {
         if (rank == 0) {
-            // check memory boarders (densities at the beginning, energies at the end)
+            // check memory borders (densities at the beginning, energies at the end)
             EXPECT_DOUBLE_EQ(rho_averages[0],         1.0);
             EXPECT_DOUBLE_EQ(e_averages[nelements-1], 1.0e-5);
             // check values somewhere near the center (expect symmetries)
@@ -140,7 +150,7 @@ TEST(CInterfaceTest, SimulationRun) {
             EXPECT_NEAR(v2_averages[94],             -0.14037267400591, 1e-14);
         }
         else {
-            // check memory boarders (densities at the beginning, energies at the end)
+            // check memory borders (densities at the beginning, energies at the end)
             EXPECT_DOUBLE_EQ(rho_averages[0],         1.0);
             EXPECT_DOUBLE_EQ(e_averages[nelements-1], 1.0e-5);
             // check values somewhere near the center (expect symmetries)