RFC: raw data access

.github/workflows/ci.yml

@@ -246,6 +246,7 @@ jobs:
           ../build/examples/trixi_controller_simple_f . libelixir_p4est2d_dgsem_euler_sedov.jl
           ../build/examples/trixi_controller_data_c . libelixir_t8code_2d_dgsem_advection_amr.jl
           ../build/examples/trixi_controller_data_f . libelixir_t8code_2d_dgsem_advection_amr.jl
+          ../build/examples/trixi_controller_raw_data_c . libelixir_t8code_2d_dgsem_advection_amr.jl
           ../build/examples/trixi_controller_t8code_c . libelixir_t8code_2d_dgsem_advection_amr.jl
           ../build/examples/trixi_controller_t8code_f . libelixir_t8code_2d_dgsem_advection_amr.jl
         env:
@@ -274,6 +275,8 @@ jobs:
                          "../build/examples/trixi_controller_data_c ."   \
                          "../build/examples/trixi_controller_data_f"     \
                          "../build/examples/trixi_controller_data_f ."   \
+                         "../build/examples/trixi_controller_raw_data_c"     \
+                         "../build/examples/trixi_controller_raw_data_c ."   \
                          "../build/examples/trixi_controller_t8code_c"   \
                          "../build/examples/trixi_controller_t8code_c ." \
                          "../build/examples/trixi_controller_t8code_f"   \

LibTrixi.jl/examples/libelixir_t8code_2d_dgsem_advection_amr.jl

@@ -51,18 +51,23 @@ function init_simstate()
     # The AMRCallback triggers adaptive mesh refinement
     amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
                                           base_level=2,
-                                          med_level=3, med_threshold=0.1,
-                                          max_level=4, max_threshold=0.6)
+                                          med_level=3, med_threshold=0.8,
+                                          max_level=4, max_threshold=1.2)
     amr_callback = AMRCallback(semi, amr_controller,
                                interval=10,
                                adapt_initial_condition=true,
                                adapt_initial_condition_only_refine=true)
 
+    save_solution = SaveSolutionCallback(interval=10,
+                                         save_initial_solution=true,
+                                         save_final_solution=true)
+
     # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
     callbacks = CallbackSet(summary_callback,
                             analysis_callback,
                             alive_callback,
                             amr_callback,
+                            save_solution,
                             stepsize_callback)
 
 

LibTrixi.jl/src/LibTrixi.jl

@@ -61,6 +61,9 @@ export trixi_load_element_averaged_primitive_vars,
 export trixi_register_data,
        trixi_register_data_cfptr,
        trixi_register_data_jl
+export trixi_get_data_pointer,
+       trixi_get_data_pointer_cfptr,
+       trixi_get_data_pointer_jl
 export trixi_version_library,
        trixi_version_library_cfptr,
        trixi_version_library_jl

LibTrixi.jl/src/api_c.jl

@@ -538,6 +538,25 @@ trixi_load_element_averaged_primitive_vars_cfptr() =
     @cfunction(trixi_load_element_averaged_primitive_vars, Cvoid, (Cint, Cint, Ptr{Cdouble}))
 
 
+"""
+    trixi_get_data_pointer(simstate_handle::Cint)::Ptr{Cdouble}
+
+Return pointer to internal data vector.
+"""
+function trixi_get_data_pointer end
+
+Base.@ccallable function trixi_get_data_pointer(simstate_handle::Cint)::Ptr{Cdouble}
+    simstate = load_simstate(simstate_handle)
+    return trixi_get_data_pointer_jl(simstate)
+end
+
+trixi_get_data_pointer_cfptr() = @cfunction(trixi_get_data_pointer, Ptr{Cdouble}, (Cint,))
+
+
+
+############################################################################################
+# t8code
+############################################################################################
 """
     trixi_get_t8code_forest(simstate_handle::Cint)::::Ptr{Trixi.t8_forest}
 
@@ -557,6 +576,8 @@ end
 trixi_get_t8code_forest_cfptr() =
     @cfunction(trixi_get_t8code_forest, Ptr{Trixi.t8_forest}, (Cint,))
 
+
+
 ############################################################################################
 # Auxiliary
 ############################################################################################

LibTrixi.jl/src/api_jl.jl

@@ -208,6 +208,11 @@ function trixi_register_data_jl(simstate, index, data)
 end
 
 
+function trixi_get_data_pointer_jl(simstate)
+    return pointer(simstate.integrator.u)
+end
+
+
 function trixi_get_simulation_time_jl(simstate)
     return simstate.integrator.t
 end

examples/CMakeLists.txt

@@ -5,6 +5,7 @@ set ( EXAMPLES
       trixi_controller_mpi.f90
       trixi_controller_data.c
       trixi_controller_data.f90
+      trixi_controller_raw_data.c
       trixi_controller_t8code.c
       trixi_controller_t8code.f90 )
 

examples/trixi_controller_raw_data.c

@@ -0,0 +1,63 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <trixi.h>
+
+int main ( int argc, char *argv[] ) {
+
+    if ( argc < 2 ) {
+        fprintf(stderr, "ERROR: missing arguments: PROJECT_DIR LIBELIXIR_PATH\n\n");
+        fprintf(stderr, "usage: %s PROJECT_DIR LIBELIXIR_PATH\n", argv[0]);
+        return 2;
+    } else if ( argc < 3 ) {
+        fprintf(stderr, "ERROR: missing argument: LIBELIXIR_PATH\n\n");
+        fprintf(stderr, "usage: %s PROJECT_DIR LIBELIXIR_PATH\n", argv[0]);
+        return 2;
+    }
+
+    // Initialize Trixi
+    printf("\n*** Trixi controller ***   Initialize Trixi\n");
+    trixi_initialize( argv[1], NULL );
+
+    // Set up the Trixi simulation
+    // We get a handle to use subsequently
+    printf("\n*** Trixi controller ***   Set up Trixi simulation\n");
+    int handle = trixi_initialize_simulation( argv[2] );
+
+    // Main loop
+    int steps = 0;
+
+    printf("\n*** Trixi controller ***   Entering main loop\n");
+    while ( !trixi_is_finished( handle ) ) {
+
+        trixi_step( handle );
+        steps++;
+
+        if (steps % 10 == 0) {
+
+            // Get number of degrees of freedom
+            int ndofs = trixi_ndofsglobal( handle );
+
+            // Get a pointer to Trixi's internal simulation data
+            double * raw_data = trixi_get_data_pointer(handle);
+
+            for (int i = 0; i < ndofs; ++i) {
+                // Compute local deviation from density 1
+                const double dev = raw_data[i] - 1.0;
+
+                // Apply 20% damping
+                raw_data[i] = 1.0 + 0.8 * dev;
+            }
+        }
+    }
+
+    // Finalize Trixi simulation
+    printf("\n*** Trixi controller ***   Finalize Trixi simulation\n");
+    trixi_finalize_simulation( handle );
+
+    // Finalize Trixi
+    printf("\n*** Trixi controller ***   Finalize Trixi\n");
+    trixi_finalize();
+
+    return 0;
+}

src/api.c

@@ -25,6 +25,7 @@ enum {
     TRIXI_FTPR_LOAD_PRIMITIVE_VARS,
     TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VARS,
     TRIXI_FTPR_REGISTER_DATA,
+    TRIXI_FPTR_GET_DATA_POINTER,
     TRIXI_FTPR_VERSION_LIBRARY,
     TRIXI_FTPR_VERSION_LIBRARY_MAJOR,
     TRIXI_FTPR_VERSION_LIBRARY_MINOR,
@@ -63,6 +64,7 @@ static const char* trixi_function_pointer_names[] = {
     [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_FPTR_GET_DATA_POINTER]                     = "trixi_get_data_pointer_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",
@@ -726,6 +728,29 @@ void trixi_register_data(int handle, int index, int size, const double * data) {
 }
 
 
+/**
+ * @anchor trixi_get_data_pointer_api_c
+ *
+ * @brief Return pointer to internal data vector.
+ *
+ * The returned pointer points to the beginning of the internal data array used in Trixi.jl.
+ * This array contains the conservative variables, i.e. density, momentum density in the
+ * three Cartesian coordinates, and energy density, in this sequence. The pointer can be
+ * used to read, but also to write these variables. The latter should be done with care.
+ * Writing while a time step in being performed will lead to undefined behavior.
+ *
+ * @param[in]  handle  simulation handle
+ */
+double * trixi_get_data_pointer(int handle) {
+
+    // Get function pointer
+    double * (*get_data_pointer)(int) = trixi_function_pointers[TRIXI_FPTR_GET_DATA_POINTER];
+
+    // Call function
+    return get_data_pointer(handle);
+}
+
+
 /**
  * @anchor trixi_get_simulation_time_api_c
  *

src/api.f90

@@ -442,6 +442,24 @@ subroutine trixi_register_data(handle, index, size, data) bind(c)
       real(c_double), dimension(*), intent(in) :: data
     end subroutine
 
+    !>
+    !! @anchor trixi_get_data_pointer_api_c
+    !!
+    !! @brief Return pointer to internal data vector.
+    !!
+    !! The returned pointer points to the beginning of the internal data array used in
+    !! Trixi.jl. This array contains the conservative variables, i.e. density, momentum
+    !! density in the three Cartesian coordinates, and energy density, in this sequence.
+    !! The pointer can be used to read, but also to write these variables. The latter
+    !! should be done with care. Writing while a time step in being performed will lead to
+    !! undefined behavior.
+    !!
+    !! @param[in]  handle  simulation handle
+    type (c_ptr) function trixi_get_data_pointer(handle) bind(c)
+      use, intrinsic :: iso_c_binding, only: c_int, c_ptr
+      integer(c_int), value, intent(in) :: handle
+    end function
+
 
 
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

src/trixi.h

@@ -40,6 +40,7 @@ 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);
+double * trixi_get_data_pointer(int handle);
 
 // T8code
 #if !defined(T8_H) && !defined(T8_FOREST_GENERAL_H)