First steps towards accessing Trixi data

.github/workflows/ci.yml

@@ -159,6 +159,10 @@ jobs:
               ../build/examples/simple_trixi_controller_c . libelixir_p4est2d_dgsem_euler_sedov.jl
           JULIA_DEPOT_PATH=~/.julia \
               ../build/examples/simple_trixi_controller_f . libelixir_p4est2d_dgsem_euler_sedov.jl
+          JULIA_DEPOT_PATH=~/.julia \
+              ../build/examples/trixi_controller_data_c . libelixir_p4est2d_dgsem_euler_sedov.jl
+          JULIA_DEPOT_PATH=~/.julia \
+              ../build/examples/trixi_controller_data_f . libelixir_p4est2d_dgsem_euler_sedov.jl
         env:
           LIBTRIXI_DEBUG: all
       - name: Run memory checks with Valgrind
@@ -174,6 +178,10 @@ jobs:
               ../build/examples/simple_trixi_controller_c . libelixir_p4est2d_dgsem_euler_sedov.jl
           JULIA_DEPOT_PATH=~/.julia valgrind --error-exitcode=1 -s \
               ../build/examples/simple_trixi_controller_f . libelixir_p4est2d_dgsem_euler_sedov.jl
+          JULIA_DEPOT_PATH=~/.julia \
+              ../build/examples/trixi_controller_data_c . libelixir_p4est2d_dgsem_euler_sedov.jl
+          JULIA_DEPOT_PATH=~/.julia \
+              ../build/examples/trixi_controller_data_f . libelixir_p4est2d_dgsem_euler_sedov.jl
       - name: Process coverage data
         if: ${{ matrix.test_type == 'coverage' }}
         run: |

LibTrixi.jl/src/LibTrixi.jl

@@ -1,7 +1,8 @@
 module LibTrixi
 
 using OrdinaryDiffEq: OrdinaryDiffEq, step!, check_error, DiscreteCallback
-using Trixi: Trixi, summary_callback
+using Trixi: Trixi, summary_callback, mesh_equations_solver_cache, nelements, nvariables,
+             nnodes, wrap_array, eachelement, cons2prim, get_node_vars, eachnode
 using MPI: MPI, run_init_hooks, set_default_error_handler_return
 
 export trixi_initialize_simulation,
@@ -19,6 +20,18 @@ export trixi_is_finished,
 export trixi_step,
        trixi_step_cfptr,
        trixi_step_jl
+export trixi_ndims,
+       trixi_ndims_cfptr,
+       trixi_ndims_jl
+export trixi_nelements,
+       trixi_nelements_cfptr,
+       trixi_nelements_jl
+export trixi_nvariables,
+       trixi_nvariables_cfptr,
+       trixi_nvariables_jl
+export trixi_load_cell_averages,
+       trixi_load_cell_averages_cfptr,
+       trixi_load_cell_averages_jl
 
 export SimulationState, store_simstate, load_simstate, delete_simstate!
 

LibTrixi.jl/src/api_c.jl

@@ -129,3 +129,75 @@ Base.@ccallable function trixi_step(simstate_handle::Cint)::Cvoid
 end
 
 trixi_step_cfptr() = @cfunction(trixi_step, Cvoid, (Cint,))
+
+
+"""
+    trixi_ndims(simstate_handle::Cint)::Cint
+
+Return number of spatial dimensions.
+"""
+function trixi_ndims end
+
+Base.@ccallable function trixi_ndims(simstate_handle::Cint)::Cint
+    simstate = load_simstate(simstate_handle)
+    return trixi_ndims_jl(simstate)
+end
+
+trixi_ndims_cfptr() = @cfunction(trixi_ndims, Cint, (Cint,))
+
+
+"""
+    trixi_nelements(simstate_handle::Cint)::Cint
+
+Return number of elements (cells).
+"""
+function trixi_nelements end
+
+Base.@ccallable function trixi_nelements(simstate_handle::Cint)::Cint
+    simstate = load_simstate(simstate_handle)
+    return trixi_nelements_jl(simstate)
+end
+
+trixi_nelements_cfptr() = @cfunction(trixi_nelements, Cint, (Cint,))
+
+
+"""
+    trixi_nvariables(simstate_handle::Cint)::Cint
+
+Return number of variables.
+"""
+function trixi_nvariables end
+
+Base.@ccallable function trixi_nvariables(simstate_handle::Cint)::Cint
+    simstate = load_simstate(simstate_handle)
+    return trixi_nvariables_jl(simstate)
+end
+
+trixi_nvariables_cfptr() = @cfunction(trixi_nvariables, Cint, (Cint,))
+
+
+"""
+    trixi_load_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid
+
+Return cell averaged solution state.
+
+Cell averaged values for each cell and each primitive variable are stored in a contiguous
+array, where cell values for the first variable appear first and values for the other
+variables subsequently (structure-of-arrays layout).
+
+The given array has to be of correct size and memory has to be allocated beforehand.
+"""
+function trixi_load_cell_averages end
+
+Base.@ccallable function trixi_load_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid
+    simstate = load_simstate(simstate_handle)
+
+    # convert C to julia array
+    size = trixi_nvariables_jl(simstate) * trixi_nelements_jl(simstate)
+    data_jl = unsafe_wrap(Array, data, size)
+
+    trixi_load_cell_averages_jl(data_jl, simstate)
+    return nothing
+end
+
+trixi_load_cell_averages_cfptr() = @cfunction(trixi_load_cell_averages, Cvoid, (Ptr{Cdouble}, Cint,))

LibTrixi.jl/src/api_jl.jl

@@ -49,3 +49,60 @@ function trixi_step_jl(simstate)
 
     return nothing
 end
+
+function trixi_ndims_jl(simstate)
+    mesh, _, _, _ = mesh_equations_solver_cache(simstate.semi)
+    return ndims(mesh)
+end
+
+function trixi_nelements_jl(simstate)
+    _, _, solver, cache = mesh_equations_solver_cache(simstate.semi)
+    return nelements(solver, cache)
+end
+
+function trixi_nvariables_jl(simstate)
+    _, equations, _, _ = mesh_equations_solver_cache(simstate.semi)
+    return nvariables(equations)
+end
+
+function trixi_load_cell_averages_jl(data, simstate)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi)
+    n_elements = nelements(solver, cache)
+    n_variables = nvariables(equations)
+    n_nodes = nnodes(solver)
+    n_dims = ndims(mesh)
+
+    u_ode = simstate.integrator.u
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+    # all permutations of nodes indices for arbitrary dimension
+    node_cis = CartesianIndices(ntuple(i -> n_nodes, n_dims))
+
+    # temporary storage for mean value on current element for all variables
+    u_mean = get_node_vars(u, equations, solver, node_cis[1], 1)
+
+    for element in eachelement(solver, cache)
+
+        # compute mean value using nodal dg values and quadrature
+        u_mean = zero(u_mean)
+        for node_ci in node_cis
+            u_node_prim = cons2prim(get_node_vars(u, equations, solver, node_ci, element), equations)
+            weight = 1.
+            for node_index in Tuple(node_ci)
+                weight *= solver.basis.weights[node_index]
+            end
+            u_mean += u_node_prim * weight
+        end
+
+        # normalize to unit element
+        u_mean = u_mean / 2^n_dims
+
+        # copy to provided array
+        # all element values for first variable, then for second, ...
+        for ivar = 0:n_variables-1
+            data[element + ivar * n_elements] = u_mean[ivar+1]
+        end
+    end
+
+    return nothing
+end

examples/CMakeLists.txt

@@ -1,6 +1,8 @@
 set ( EXAMPLES
       simple_trixi_controller.c
-      simple_trixi_controller.f90 )
+      simple_trixi_controller.f90
+      trixi_controller_data.c
+      trixi_controller_data.f90 )
 
 
 foreach ( EXAMPLE ${EXAMPLES} )

examples/trixi_controller_data.c

@@ -0,0 +1,64 @@
+#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
+    printf("\n*** Trixi controller ***   Entering main loop\n");
+    while ( !trixi_is_finished( handle ) ) {
+
+        trixi_step( handle );
+    }
+
+    // get number of elements
+    int nelements = trixi_nelements( handle );
+    printf("\n*** Trixi controller ***   nelements %d\n", nelements);
+
+    // get number of variables
+    int nvariables = trixi_nvariables( handle );
+    printf("\n*** Trixi controller ***   nvariables %d\n", nvariables);
+
+    // allocate memory
+    double* data = malloc( sizeof(double) * nelements * nvariables );
+
+    // get averaged cell values for each variable
+    trixi_load_cell_averages(data, handle);
+
+    // compute temperature
+    const double gas_constant = 0.287;
+
+    for (int i = 0; i < nelements; ++i) {
+        printf("T[cell %3d] = %f\n", i, data[i+3*nelements] / (gas_constant * data[i]) );
+    }
+
+    // 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;
+}

examples/trixi_controller_data.f90

@@ -0,0 +1,78 @@
+program simple_trixi_controller_f
+  use LibTrixi
+  use, intrinsic :: iso_fortran_env, only: error_unit
+  use, intrinsic :: iso_c_binding, only: c_int, c_double
+
+  implicit none
+
+  integer(c_int) :: handle, nelements, nvariables, i
+  character(len=256) :: argument
+  integer, parameter :: dp = selected_real_kind(12)
+  real(dp), dimension(:), pointer :: data
+  real(c_double) :: gas_constant
+
+
+  if (command_argument_count() < 1) then
+    call get_command_argument(0, argument)
+    write(error_unit, '(a)') "ERROR: missing arguments: PROJECT_DIR LIBELIXIR_PATH"
+    write(error_unit, '(a)') ""
+    write(error_unit, '(3a)') "usage: ", trim(argument), " PROJECT_DIR LIBELIXIR_PATH"
+    stop 1
+  else if (command_argument_count() < 2) then
+    call get_command_argument(0, argument)
+    write(error_unit, '(a)') "ERROR: missing argument: LIBELIXIR_PATH"
+    write(error_unit, '(a)') ""
+    write(error_unit, '(3a)') "usage: ", trim(argument), " PROJECT_DIR LIBELIXIR_PATH"
+    stop 1
+  end if
+
+
+  ! Initialize Trixi
+  call get_command_argument(1, argument)
+  call trixi_initialize(argument)
+
+  ! Set up the Trixi simulation
+  ! We get a handle to use subsequently
+  call get_command_argument(2, argument)
+  handle = trixi_initialize_simulation(argument)
+
+  ! Main loop
+  do
+    ! Exit loop once simulation is completed
+    if ( trixi_is_finished(handle) ) exit
+
+    call trixi_step(handle)
+  end do
+
+  ! get number of elements
+  nelements = trixi_nelements(handle);
+  write(*, '(a,i6)') "*** Trixi controller ***   nelements ", nelements
+  write(*, '(a)') ""
+
+  ! get number of variables
+  nvariables = trixi_nvariables( handle );
+  write(*, '(a,i6)') "*** Trixi controller ***   nvariables ", nvariables
+  write(*, '(a)') ""
+
+  ! allocate memory
+  allocate ( data(nelements*nvariables) )
+
+  ! get averaged cell values for each variable
+  call trixi_load_cell_averages(data, handle);
+
+  ! compute temperature
+  gas_constant = 0.287;
+
+  do i = 1,nelements
+    print "('T[cell  ', i4, '] = ', e14.8)", i, data(i+3*nelements)/(gas_constant*data(i))
+  end do
+
+  write(*, '(a,i6)') "*** Trixi controller ***   Finalize Trixi simulation "
+  write(*, '(a)') ""
+
+  ! Finalize Trixi simulation
+  call trixi_finalize_simulation(handle)
+
+  ! Finalize Trixi
+  call trixi_finalize()
+end program