change trixi_load_cell_average

add parameter index and only get averaged of the variable at position index

LibTrixi.jl/src/api_c.jl

@@ -358,32 +358,31 @@ trixi_nvariables_cfptr() = @cfunction(trixi_nvariables, Cint, (Cint,))
 
 
 """
-    trixi_load_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid
+    trixi_load_cell_averages(data::Ptr{Cdouble}, index::Cint, 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).
+Cell averaged values for the primitive variable at position `index` for each cell are stored
+in the given array `data`.
 
 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},
+Base.@ccallable function trixi_load_cell_averages(data::Ptr{Cdouble}, index::Cint,
                                                   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)
+    trixi_load_cell_averages_jl(data_jl, index, simstate)
     return nothing
 end
 
 trixi_load_cell_averages_cfptr() =
-    @cfunction(trixi_load_cell_averages, Cvoid, (Ptr{Cdouble}, Cint,))
+    @cfunction(trixi_load_cell_averages, Cvoid, (Ptr{Cdouble}, Cint, Cint,))
 
 
 """

LibTrixi.jl/src/api_jl.jl

@@ -100,10 +100,8 @@ function trixi_nvariables_jl(simstate)
 end
 
 
-function trixi_load_cell_averages_jl(data, simstate)
+function trixi_load_cell_averages_jl(data, index, 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)
 
@@ -113,15 +111,13 @@ function trixi_load_cell_averages_jl(data, simstate)
     # 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)
+        u_mean = zero(eltype(u))
         for node_ci in node_cis
-            u_node_prim = cons2prim(get_node_vars(u, equations, solver, node_ci, element), equations)
+            u_node_prim = cons2prim(get_node_vars(u, equations, solver, node_ci, element),
+                                    equations)[index]
             weight = 1.
             for node_index in Tuple(node_ci)
                 weight *= solver.basis.weights[node_index]
@@ -132,11 +128,8 @@ function trixi_load_cell_averages_jl(data, simstate)
         # 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
+        # write to provided array
+        data[element] = u_mean
     end
 
     return nothing

LibTrixi.jl/test/test_interface.jl

@@ -87,10 +87,10 @@ end
     @test nvariables_c == nvariables_jl
 
     # compare cell averaged values
-    data_c = zeros(nvariables_c * nelements_c)
-    trixi_load_cell_averages(pointer(data_c), handle)
-    data_jl = zeros(nvariables_jl * nelements_jl)
-    trixi_load_cell_averages_jl(data_jl, simstate_jl)
+    data_c = zeros(nelements_c)
+    trixi_load_cell_averages(pointer(data_c), 1, handle)
+    data_jl = zeros(nelements_jl)
+    trixi_load_cell_averages_jl(data_jl, 1, simstate_jl)
     @test data_c == data_jl
 end
 

examples/trixi_controller_data.c

@@ -46,10 +46,10 @@ int main ( int argc, char *argv[] ) {
             printf("\n*** Trixi controller ***   nelements %d\n", nelements);
 
             // Allocate memory
-            data = realloc( data, sizeof(double) * nelements * nvariables );
+            data = realloc( data, sizeof(double) * nelements );
 
-            // Get averaged cell values for each variable
-            trixi_load_cell_averages(data, handle);
+            // Get averaged cell values for first variable
+            trixi_load_cell_averages(data, 1, handle);
         }
     }
 

examples/trixi_controller_data.f90

@@ -63,10 +63,10 @@ program trixi_controller_data_f
 
       ! allocate memory
       if ( associated(data) ) deallocate(data)
-      allocate( data(nelements*nvariables) )
+      allocate( data(nelements) )
 
-      ! get averaged cell values for each variable
-      call trixi_load_cell_averages(data, handle)
+      ! get averaged cell values for first variable
+      call trixi_load_cell_averages(data, 1, handle)
     end if
   end do
 

src/api.c

@@ -554,29 +554,30 @@ int trixi_nvariables(int handle) {
 /** 
  * @anchor trixi_load_cell_averages_api_c
  *
- * @brief Return cell averaged values
+ * @brief 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).
+ * Cell averaged values for the primitive variable at position index for each cell are
+ * stored in the given array `data`.
  *
  * The given array has to be of correct size and memory has to be allocated beforehand.
  *
  * @param[in]  handle  simulation handle
+ * @param[in]  index   index of variable
  * @param[out] data    cell averaged values for all cells and all primitive variables
  */
-void trixi_load_cell_averages(double * data, int handle) {
+void trixi_load_cell_averages(double * data, int index, int handle) {
 
     // Get function pointer
-    void (*load_cell_averages)(double *, int) =
+    void (*load_cell_averages)(double *, int, int) =
         trixi_function_pointers[TRIXI_FTPR_LOAD_CELL_AVERAGES];
 
     // Call function
-    load_cell_averages(data, handle);
+    load_cell_averages(data, index, handle);
 }
 
 
 
+
 /******************************************************************************************/
 /* T8code                                                                                 */
 /******************************************************************************************/

src/api.f90

@@ -311,17 +311,19 @@ integer(c_int) function trixi_nvariables(handle) bind(c)
     end function
 
     !>
-    !! @fn LibTrixi::trixi_load_cell_averages::trixi_load_cell_averages(data, handle)
+    !! @fn LibTrixi::trixi_load_cell_averages::trixi_load_cell_averages(data, index, handle)
     !!
     !! @brief Return cell averaged values
     !!
     !! @param[in]  handle  simulation handle
+    !! @param[in]  index   index of variable
     !! @param[out] data    cell averaged values for all cells and all variables
     !!
     !! @see @ref trixi_load_cell_averages_api_c "trixi_load_cell_averages (C API)"
-    subroutine trixi_load_cell_averages(data, handle) bind(c)
+    subroutine trixi_load_cell_averages(data, index, handle) bind(c)
       use, intrinsic :: iso_c_binding, only: c_int, c_double
       real(c_double), dimension(*), intent(in) :: data
+      integer(c_int), value, intent(in) :: index
       integer(c_int), value, intent(in) :: handle
     end subroutine
 

src/trixi.h

@@ -32,7 +32,7 @@ int trixi_ndofs(int handle);
 int trixi_ndofs_global(int handle);
 int trixi_nvariables(int handle);
 double trixi_calculate_dt(int handle);
-void trixi_load_cell_averages(double * data, int handle);
+void trixi_load_cell_averages(double * data, int index, int handle);
 
 // T8code
 #if !defined(T8_H) && !defined(T8_FOREST_GENERAL_H)

test/c/simulation.cpp

@@ -67,49 +67,54 @@ TEST(CInterfaceTest, SimulationRun) {
     EXPECT_EQ(nvariables, 4);
 
     // Check cell averaged values
-    int size = nelements * nvariables;
-    std::vector<double> cell_averages(size);
-    trixi_load_cell_averages(cell_averages.data(), handle);
+    std::vector<double> rho_averages(nelements);
+    std::vector<double> v1_averages(nelements);
+    std::vector<double> v2_averages(nelements);
+    std::vector<double> e_averages(nelements);
+    trixi_load_cell_averages(rho_averages.data(), 1, handle);
+    trixi_load_cell_averages(v1_averages.data(), 2, handle);
+    trixi_load_cell_averages(v2_averages.data(), 3, handle);
+    trixi_load_cell_averages(e_averages.data(), 4, handle);
     if (nranks == 1) {
         // check memory boarders (densities at the beginning, energies at the end)
-        EXPECT_DOUBLE_EQ(cell_averages[0],              1.0);
-        EXPECT_DOUBLE_EQ(cell_averages[size-1],         1.0e-5);
+        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)
         // densities
-        EXPECT_NEAR(cell_averages[93],                  0.88263491354796, 1e-14);
-        EXPECT_NEAR(cell_averages[94],                  0.88263491354796, 1e-14);
-        EXPECT_NEAR(cell_averages[161],                 0.88263491354796, 1e-14);
-        EXPECT_NEAR(cell_averages[162],                 0.88263491354796, 1e-14);
+        EXPECT_NEAR(rho_averages[ 93],            0.88263491354796, 1e-14);
+        EXPECT_NEAR(rho_averages[ 94],            0.88263491354796, 1e-14);
+        EXPECT_NEAR(rho_averages[161],            0.88263491354796, 1e-14);
+        EXPECT_NEAR(rho_averages[162],            0.88263491354796, 1e-14);
         // velocities
-        EXPECT_NEAR(cell_averages[ 93+  nelements],    -0.14037267400591, 1e-14);
-        EXPECT_NEAR(cell_averages[ 94+2*nelements],    -0.14037267400591, 1e-14);
-        EXPECT_NEAR(cell_averages[161+2*nelements],     0.14037267400591, 1e-14);
-        EXPECT_NEAR(cell_averages[162+  nelements],     0.14037267400591, 1e-14);
+        EXPECT_NEAR(v1_averages[ 93],            -0.14037267400591, 1e-14);
+        EXPECT_NEAR(v2_averages[ 94],            -0.14037267400591, 1e-14);
+        EXPECT_NEAR(v2_averages[161],             0.14037267400591, 1e-14);
+        EXPECT_NEAR(v1_averages[162],             0.14037267400591, 1e-14);
     }
     else if (nranks == 2) {
         if (rank == 0) {
             // check memory boarders (densities at the beginning, energies at the end)
-            EXPECT_DOUBLE_EQ(cell_averages[0],          1.0);
-            EXPECT_DOUBLE_EQ(cell_averages[size-1],     1.0e-5);
+            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)
             // densities
-            EXPECT_NEAR(cell_averages[93],              0.88263491354796, 1e-14);
-            EXPECT_NEAR(cell_averages[94],              0.88263491354796, 1e-14);
+            EXPECT_NEAR(rho_averages[93],             0.88263491354796, 1e-14);
+            EXPECT_NEAR(rho_averages[94],             0.88263491354796, 1e-14);
             // velocities
-            EXPECT_NEAR(cell_averages[93+  nelements], -0.14037267400591, 1e-14);
-            EXPECT_NEAR(cell_averages[94+2*nelements], -0.14037267400591, 1e-14);
+            EXPECT_NEAR(v1_averages[93],             -0.14037267400591, 1e-14);
+            EXPECT_NEAR(v2_averages[94],             -0.14037267400591, 1e-14);
         }
         else {
             // check memory boarders (densities at the beginning, energies at the end)
-            EXPECT_DOUBLE_EQ(cell_averages[0],          1.0);
-            EXPECT_DOUBLE_EQ(cell_averages[size-1],     1.0e-5);
+            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)
             // densities
-            EXPECT_NEAR(cell_averages[33],              0.88263491354796, 1e-14);
-            EXPECT_NEAR(cell_averages[34],              0.88263491354796, 1e-14);
+            EXPECT_NEAR(rho_averages[33],             0.88263491354796, 1e-14);
+            EXPECT_NEAR(rho_averages[34],             0.88263491354796, 1e-14);
             // velocities
-            EXPECT_NEAR(cell_averages[33+2*nelements],  0.14037267400591, 1e-14);
-            EXPECT_NEAR(cell_averages[34+  nelements],  0.14037267400591, 1e-14);
+            EXPECT_NEAR(v2_averages[33],              0.14037267400591, 1e-14);
+            EXPECT_NEAR(v1_averages[34],              0.14037267400591, 1e-14);
         }
     }
     else {

test/fortran/simulationRun_suite.f90

@@ -63,13 +63,13 @@ subroutine test_simulationRun(error)
     call check(error, nvariables, 4)
 
     ! Check cell averaged values
-    size = nelements*nvariables
+    size = nelements
     allocate(data(size))
-    call trixi_load_cell_averages(data, handle)
+    call trixi_load_cell_averages(data, 1, handle)
     call check(error, data(1), 1.0_dp)
-    call check(error, data(929), 2.6605289164377273_dp)
+    call check(error, data(94), 0.88263491354796_dp)
     call check(error, data(size), 1e-5_dp)
-    
+
     ! Finalize Trixi simulation
     call trixi_finalize_simulation(handle)