From 70b19bd89f2d97dd483b44a31a206b91a99d82dc Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 5 Jul 2023 13:17:44 +0200 Subject: [PATCH 01/25] first steps towards accessing Trixi data --- .../examples/libelixir_tree2d_dgsem_euler.jl | 74 +++++++++++++++++++ examples/trixi_controller_data.c | 65 ++++++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl create mode 100644 examples/trixi_controller_data.c diff --git a/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl b/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl new file mode 100644 index 00000000..ed1875c7 --- /dev/null +++ b/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl @@ -0,0 +1,74 @@ + +using LibTrixi +using Trixi +using OrdinaryDiffEq + +# The function to create the simulation state needs to be named `init_simstate` +function init_simstate() + + ############################################################################### + # semidiscretization of the compressible Euler equations + + equations = CompressibleEulerEquations2D(1.4) + + initial_condition = initial_condition_convergence_test + + source_terms = source_terms_convergence_test + + boundary_condition = BoundaryConditionDirichlet(initial_condition) + + boundary_conditions = (x_neg=boundary_condition, + x_pos=boundary_condition, + y_neg=boundary_condition, + y_pos=boundary_condition,) + + solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + + coordinates_min = (0.0, 0.0) + coordinates_max = (2.0, 2.0) + + mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=10_000) + + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms=source_terms, + boundary_conditions=boundary_conditions) + + + ############################################################################### + # ODE solvers, callbacks etc. + + tspan = (0.0, 1.0) + ode = semidiscretize(semi, tspan) + + summary_callback = SummaryCallback() + + analysis_interval = 100 + #analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + + alive_callback = AliveCallback(analysis_interval=analysis_interval) + + stepsize_callback = StepsizeCallback(cfl=0.8) + + callbacks = CallbackSet(summary_callback, + alive_callback, + stepsize_callback) + + + ############################################################################### + # create the time integrator + + integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); + + + ############################################################################### + # Create simulation state + + simstate = SimulationState(semi, integrator) + + return simstate + +end diff --git a/examples/trixi_controller_data.c b/examples/trixi_controller_data.c new file mode 100644 index 00000000..1ca7fbff --- /dev/null +++ b/examples/trixi_controller_data.c @@ -0,0 +1,65 @@ +#include +#include + +#include + +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] ); + + // 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_get_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; +} From 8affd33cc5153894f5f29322d5d5dedaf0db26e1 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Wed, 5 Jul 2023 13:19:51 +0200 Subject: [PATCH 02/25] the real changes --- LibTrixi.jl/src/LibTrixi.jl | 4 ++ LibTrixi.jl/src/api_c.jl | 61 ++++++++++++++++++++++++++++++ LibTrixi.jl/src/api_jl.jl | 50 +++++++++++++++++++++++++ examples/CMakeLists.txt | 22 +++++++++++ src/trixi.c | 75 +++++++++++++++++++++++++++++++++++-- src/trixi.h | 21 +++++++++-- 6 files changed, 227 insertions(+), 6 deletions(-) diff --git a/LibTrixi.jl/src/LibTrixi.jl b/LibTrixi.jl/src/LibTrixi.jl index ade28539..d559d5b7 100644 --- a/LibTrixi.jl/src/LibTrixi.jl +++ b/LibTrixi.jl/src/LibTrixi.jl @@ -18,6 +18,10 @@ export trixi_is_finished, export trixi_step, trixi_step_cfptr, trixi_step_jl +export trixi_ndims_cfptr +export trixi_nelements_cfptr +export trixi_nvariables_cfptr +export trixi_get_cell_averages_cfptr export SimulationState, store_simstate, load_simstate, delete_simstate! diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index fcee4965..cc2793a5 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -129,3 +129,64 @@ 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 elements (cells) +""" +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_get_cell_averages(simstate_handle::Cint)::Cint + +Return number of elements (cells) +""" +function trixi_get_cell_averages end + +Base.@ccallable function trixi_get_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid + simstate = load_simstate(simstate_handle) + trixi_get_cell_averages_jl(data, simstate) + return nothing +end + +trixi_get_cell_averages_cfptr() = @cfunction(trixi_get_cell_averages, Cvoid, (Ptr{Cdouble}, Cint,)) \ No newline at end of file diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index 518d3ad4..1f2186f0 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -49,3 +49,53 @@ function trixi_step_jl(simstate) return nothing end + +function trixi_ndims_jl(simstate) + mesh, _, _, _ = Trixi.mesh_equations_solver_cache(simstate.semi) + return ndims(mesh) +end + +function trixi_nelements_jl(simstate) + _, _, solver, cache = Trixi.mesh_equations_solver_cache(simstate.semi) + return Trixi.nelements(solver, cache) +end + +function trixi_nvariables_jl(simstate) + _, equations, _, api_jl = Trixi.mesh_equations_solver_cache(simstate.semi) + return Trixi.nvariables(equations) +end + +function trixi_get_cell_averages_jl(data_, simstate) + + mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(simstate.semi) + nelements = Trixi.nelements(solver, cache) + nvariables = Trixi.nvariables(equations) + + # convert C to julia + data = unsafe_wrap(Array, data_, nelements*nvariables) + + u_ode = simstate.integrator.u + u = Trixi.wrap_array(u_ode, mesh, equations, solver, cache) + + for (index, element) in enumerate(Trixi.eachelement(solver, cache)) + + # temporary storage for mean value on current element for all variables + u_mean = zero(Trixi.get_node_vars(u, equations, solver, 1, 1, element)) + + # compute mean value using nodal dg values and quadrature + for j in Trixi.eachnode(solver), i in Trixi.eachnode(solver) + u_node_prim = Trixi.cons2prim(Trixi.get_node_vars(u, equations, solver, i, j, element), equations) + u_mean += u_node_prim * solver.basis.weights[i] * solver.basis.weights[j] + end + + # normalize to unit element + u_mean = u_mean / 2^ndims(mesh) + + # copy to provided array + for ivar = 0:nvariables-1 + data[index + ivar * nelements] = u_mean[ivar+1] + end + end + + return nothing +end diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index af6a2f88..e4dd0b8d 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -41,3 +41,25 @@ target_include_directories( set_target_properties( simple_trixi_controller_f PROPERTIES INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") install( TARGETS simple_trixi_controller_f ) + + + +# +# Trixi Controller with Data (C version) +# +add_executable(trixi_controller_data_c trixi_controller_data.c) + +target_link_libraries( + trixi_controller_data_c + PRIVATE ${MPI_C_LIBRARIES} ${PROJECT_NAME} ${PROJECT_NAME}_tls +) + +target_include_directories( + trixi_controller_data_c + PRIVATE ${CMAKE_SOURCE_DIR}/src ${MPI_C_INCLUDE_DIRS} +) + +# Set runtime path for installed binaries +set_target_properties( trixi_controller_data_c PROPERTIES INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") + +install( TARGETS trixi_controller_data_c ) diff --git a/src/trixi.c b/src/trixi.c index ba6b66f5..bfad54d5 100644 --- a/src/trixi.c +++ b/src/trixi.c @@ -23,6 +23,10 @@ enum { TRIXI_FTPR_IS_FINISHED, TRIXI_FTPR_STEP, TRIXI_FTPR_FINALIZE_SIMULATION, + TRIXI_FTPR_NDIMS, + TRIXI_FTPR_NELEMENTS, + TRIXI_FTPR_NVARIABLES, + TRIXI_FTPR_GET_CELL_AVERAGES, // The last one is for the array size TRIXI_NUM_FPTRS @@ -37,13 +41,17 @@ static const char* trixi_function_pointer_names[] = { [TRIXI_FTPR_IS_FINISHED] = "trixi_is_finished_cfptr", [TRIXI_FTPR_STEP] = "trixi_step_cfptr", [TRIXI_FTPR_FINALIZE_SIMULATION] = "trixi_finalize_simulation_cfptr", + [TRIXI_FTPR_NDIMS] = "trixi_ndims_cfptr", + [TRIXI_FTPR_NELEMENTS] = "trixi_nelements_cfptr", + [TRIXI_FTPR_NVARIABLES] = "trixi_nvariables_cfptr", + [TRIXI_FTPR_GET_CELL_AVERAGES] = "trixi_get_cell_averages_cfptr", }; /** Initialize Julia runtime environment * - * \todo Path is still hard-coded + * \param[in] project_directory path to julia project directory */ void trixi_initialize(const char * project_directory) { @@ -172,8 +180,6 @@ void trixi_finalize_simulation(int handle) { /** Finalize Julia runtime environment - * - * \param[in] handle simulation handle to release */ void trixi_finalize() { @@ -190,6 +196,69 @@ void trixi_finalize() { } +/** Return number of spatial dimensions + * + * \param[in] handle simulation handle to release + */ +int trixi_ndims(int handle) { + + // Get function pointer + int (*ndims)(int) = trixi_function_pointers[TRIXI_FTPR_NDIMS]; + + // Call function + return ndims(handle); +} + + +/** Return number of elements (cells) + * + * \param[in] handle simulation handle to release + */ +int trixi_nelements(int handle) { + + // Get function pointer + int (*nelements)(int) = trixi_function_pointers[TRIXI_FTPR_NELEMENTS]; + + // Call function + return nelements(handle); +} + + +/** Return number of (conservative) variables + * + * \param[in] handle simulation handle to release + */ +int trixi_nvariables(int handle) { + + // Get function pointer + int (*nvariables)(int) = trixi_function_pointers[TRIXI_FTPR_NVARIABLES]; + + // Call function + return nvariables(handle); +} + +// int trixi_polydeg(int handle); // Return polynomial degree of DGSEM approximation +// int trixi_ndofs(int handle); // Return total number of degrees of freedom +// int trixi_ndofs_element(int handle); // Return number of degrees of freedom for one element + + +/** Return cell averaged values + * + * Cell averaged values for each cell and each variable are stored in a contiguous array. + * The given array has to be of correct size and memory has to be allocated beforehand. + * + * \param[in] handle simulation handle to release + * \param[out] data cell averaged values for all cells and all variables + */ +void trixi_get_cell_averages(double * data, int handle) { + + // Get function pointer + void (*get_cell_averages)(double *, int) = trixi_function_pointers[TRIXI_FTPR_GET_CELL_AVERAGES]; + + // Call function + get_cell_averages(data, handle); +} + void julia_eval_string(const char * code) { diff --git a/src/trixi.h b/src/trixi.h index ab0f34c4..af01d75e 100644 --- a/src/trixi.h +++ b/src/trixi.h @@ -1,14 +1,29 @@ #ifndef TRIXI_H_ #define TRIXI_H_ +// Setup void trixi_initialize(const char * project_directory); int trixi_initialize_simulation(const char * libelixir); -double trixi_calculate_dt(int handle); -int trixi_is_finished(int handle); -void trixi_step(int handle); void trixi_finalize_simulation(int handle); void trixi_finalize(); +// Flow control +int trixi_is_finished(int handle); +void trixi_step(int handle); + +// Basic querying +int trixi_ndims(int handle); // Return number of spatial dimensions +int trixi_nelements(int handle); // Return number of elements (cells) +// int trixi_polydeg(int handle); // Return polynomial degree of DGSEM approximation +int trixi_nvariables(int handle); // Return number of (conservative) variables +// int trixi_ndofs(int handle); // Return total number of degrees of freedom +// int trixi_ndofs_element(int handle); // Return number of degrees of freedom for one element + +// Data +void trixi_get_cell_averages(double * data, int handle); + +// Misc +double trixi_calculate_dt(int handle); void julia_eval_string(const char * code); #endif // ifndef LIBTRIXI_H_ From b6f1858daa36e556a8b9a1660be7a28882d55c99 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 12:43:11 +0200 Subject: [PATCH 03/25] adapt to current main --- LibTrixi.jl/src/LibTrixi.jl | 16 ++++++++++++---- examples/CMakeLists.txt | 2 +- examples/trixi_controller_data.c | 2 +- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/LibTrixi.jl/src/LibTrixi.jl b/LibTrixi.jl/src/LibTrixi.jl index d559d5b7..852de171 100644 --- a/LibTrixi.jl/src/LibTrixi.jl +++ b/LibTrixi.jl/src/LibTrixi.jl @@ -18,10 +18,18 @@ export trixi_is_finished, export trixi_step, trixi_step_cfptr, trixi_step_jl -export trixi_ndims_cfptr -export trixi_nelements_cfptr -export trixi_nvariables_cfptr -export trixi_get_cell_averages_cfptr +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_get_cell_averages, + trixi_get_cell_averages_cfptr, + trixi_get_cell_averages_jl export SimulationState, store_simstate, load_simstate, delete_simstate! diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index a3d4ca88..7b09d516 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,6 +1,6 @@ set ( EXAMPLES simple_trixi_controller.c - simple_trixi_controller.f9 + simple_trixi_controller.f90 trixi_controller_data.c ) diff --git a/examples/trixi_controller_data.c b/examples/trixi_controller_data.c index 1ca7fbff..3f120faf 100644 --- a/examples/trixi_controller_data.c +++ b/examples/trixi_controller_data.c @@ -17,7 +17,7 @@ int main ( int argc, char *argv[] ) { // Initialize Trixi printf("\n*** Trixi controller *** Initialize Trixi\n"); - trixi_initialize( argv[1] ); + trixi_initialize( argv[1], NULL ); // Set up the Trixi simulation // We get a handle to use subsequently From 7a39eea49252e79ee9e5abe18981e5b16c692ef4 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 12:43:46 +0200 Subject: [PATCH 04/25] add queries and data access to fortran interface --- src/trixi.f90 | 54 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/trixi.f90 b/src/trixi.f90 index 05cdc684..0d85be41 100644 --- a/src/trixi.f90 +++ b/src/trixi.f90 @@ -114,6 +114,60 @@ subroutine trixi_finalize_simulation(handle) bind(c) subroutine trixi_finalize() bind(c) end subroutine + !> + !! @fn LibTrixi::trixi_ndims::trixi_ndims(handle) + !! + !! @brief Return number of spatial dimensions + !! + !! @param[in] handle simulation handle + !! + !! @see @ref trixi_ndims_api_c "trixi_ndims (C API)" + integer(c_int) function trixi_ndims(handle) bind(c) + use, intrinsic :: iso_c_binding, only: c_int + integer(c_int), value, intent(in) :: handle + end function + + !> + !! @fn LibTrixi::trixi_nelements::trixi_nelements(handle) + !! + !! @brief Return number of elements (cells) + !! + !! @param[in] handle simulation handle + !! + !! @see @ref trixi_nelements_api_c "trixi_nelements (C API)" + integer(c_int) function trixi_nelements(handle) bind(c) + use, intrinsic :: iso_c_binding, only: c_int + integer(c_int), value, intent(in) :: handle + end function + + !> + !! @fn LibTrixi::trixi_nvariables::trixi_nvariables(handle) + !! + !! @brief Return number of (conservative) variables + !! + !! @param[in] handle simulation handle + !! + !! @see @ref trixi_nvariables_api_c "trixi_nvariables (C API" + integer(c_int) function trixi_nvariables(handle) bind(c) + use, intrinsic :: iso_c_binding, only: c_int + integer(c_int), value, intent(in) :: handle + end function + + !> + !! @fn LibTrixi::trixi_get_cell_averages::trixi_get_cell_averages(data, handle) + !! + !! @brief Return cell averaged values + !! + !! @param[in] handle simulation handle + !! @param[out] data cell averaged values for all cells and all variables + !! + !! @see @ref trixi_get_cell_averages_api_c "trixi_get_cell_averages (C API)" + subroutine trixi_get_cell_averages(data, handle) bind(c) + use, intrinsic :: iso_c_binding, only: c_int + integer(c_int), intent(in) :: data(*) + integer(c_int), value, intent(in) :: handle + end subroutine + !> !! @fn LibTrixi::julia_eval_string_c::julia_eval_string_c(code) !! From 9ad09582db15e319331a14659a1c48c7e9679db8 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 13:46:52 +0200 Subject: [PATCH 05/25] start with Fortran example for data handling --- examples/CMakeLists.txt | 3 +- examples/trixi_controller_data.f90 | 74 ++++++++++++++++++++++++++++++ src/trixi.f90 | 4 +- 3 files changed, 78 insertions(+), 3 deletions(-) create mode 100644 examples/trixi_controller_data.f90 diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 7b09d516..c4499852 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,7 +1,8 @@ set ( EXAMPLES simple_trixi_controller.c simple_trixi_controller.f90 - trixi_controller_data.c ) + trixi_controller_data.c + trixi_controller_data.f90 ) foreach ( EXAMPLE ${EXAMPLES} ) diff --git a/examples/trixi_controller_data.f90 b/examples/trixi_controller_data.f90 new file mode 100644 index 00000000..977d2da8 --- /dev/null +++ b/examples/trixi_controller_data.f90 @@ -0,0 +1,74 @@ +program simple_trixi_controller_f + use LibTrixi + use, intrinsic :: iso_fortran_env, only: error_unit + use, intrinsic :: iso_c_binding, only: c_int, c_null_char, c_double + + implicit none + + integer(c_int) :: handle, nelements, nvariables, i + character(len=256) :: argument + real(c_double), allocatable, target :: u(:) + 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 ", nelements + write(*, '(a)') "" + + ! allocate memory + allocate ( u(0:nelements*nvariables) ) + + ! get averaged cell values for each variable + call trixi_get_cell_averages(u,handle); + + ! compute temperature + gas_constant = 0.287; + + do i = 1,size(u) + print "('T[cell ' i5, '] = ', e10.10)", i, u(i+3*nelements)/(gas_constant*u(i)) + end do + + ! Finalize Trixi simulation + call trixi_finalize_simulation(handle) + + ! Finalize Trixi + call trixi_finalize() +end program diff --git a/src/trixi.f90 b/src/trixi.f90 index 0d85be41..e04272c0 100644 --- a/src/trixi.f90 +++ b/src/trixi.f90 @@ -163,8 +163,8 @@ integer(c_int) function trixi_nvariables(handle) bind(c) !! !! @see @ref trixi_get_cell_averages_api_c "trixi_get_cell_averages (C API)" subroutine trixi_get_cell_averages(data, handle) bind(c) - use, intrinsic :: iso_c_binding, only: c_int - integer(c_int), intent(in) :: data(*) + use, intrinsic :: iso_c_binding, only: c_int, c_double + real(c_double), intent(in) :: data(*) integer(c_int), value, intent(in) :: handle end subroutine From 155e895362888a77a8b2cebe89e6cecc43bf88da Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 16:02:48 +0200 Subject: [PATCH 06/25] set refinement level to 2 for testing --- LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl b/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl index ed1875c7..92677e73 100644 --- a/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl +++ b/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl @@ -28,7 +28,7 @@ function init_simstate() coordinates_max = (2.0, 2.0) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=4, + initial_refinement_level=2, n_cells_max=10_000) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, From b5ddfcbd2833dc2a1f63ae4e73a6f359d5736d09 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 16:03:11 +0200 Subject: [PATCH 07/25] enumerate is not necessary here --- LibTrixi.jl/src/api_jl.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index 1f2186f0..6de8dcea 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -77,7 +77,7 @@ function trixi_get_cell_averages_jl(data_, simstate) u_ode = simstate.integrator.u u = Trixi.wrap_array(u_ode, mesh, equations, solver, cache) - for (index, element) in enumerate(Trixi.eachelement(solver, cache)) + for element in Trixi.eachelement(solver, cache) # temporary storage for mean value on current element for all variables u_mean = zero(Trixi.get_node_vars(u, equations, solver, 1, 1, element)) @@ -93,7 +93,7 @@ function trixi_get_cell_averages_jl(data_, simstate) # copy to provided array for ivar = 0:nvariables-1 - data[index + ivar * nelements] = u_mean[ivar+1] + data[element + ivar * nelements] = u_mean[ivar+1] end end From 5f99e56c1279a051d8918f3af85080c751cdb24b Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 16:03:48 +0200 Subject: [PATCH 08/25] get Fortran example working --- examples/trixi_controller_data.f90 | 17 ++++++++++------- src/trixi.f90 | 4 ++-- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/examples/trixi_controller_data.f90 b/examples/trixi_controller_data.f90 index 977d2da8..284d6fe0 100644 --- a/examples/trixi_controller_data.f90 +++ b/examples/trixi_controller_data.f90 @@ -1,13 +1,13 @@ program simple_trixi_controller_f use LibTrixi use, intrinsic :: iso_fortran_env, only: error_unit - use, intrinsic :: iso_c_binding, only: c_int, c_null_char, c_double + use, intrinsic :: iso_c_binding, only: c_int, c_null_char, c_double, c_loc implicit none integer(c_int) :: handle, nelements, nvariables, i character(len=256) :: argument - real(c_double), allocatable, target :: u(:) + real(c_double), allocatable, target :: data(:) real(c_double) :: gas_constant @@ -50,22 +50,25 @@ program simple_trixi_controller_f ! get number of variables nvariables = trixi_nvariables( handle ); - write(*, '(a,i6)') "*** Trixi controller *** nvariables ", nelements + write(*, '(a,i6)') "*** Trixi controller *** nvariables ", nvariables write(*, '(a)') "" ! allocate memory - allocate ( u(0:nelements*nvariables) ) + allocate ( data(0:nelements*nvariables) ) ! get averaged cell values for each variable - call trixi_get_cell_averages(u,handle); + call trixi_get_cell_averages(c_loc(data),handle); ! compute temperature gas_constant = 0.287; - do i = 1,size(u) - print "('T[cell ' i5, '] = ', e10.10)", i, u(i+3*nelements)/(gas_constant*u(i)) + do i = 0,nelements-1 + 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) diff --git a/src/trixi.f90 b/src/trixi.f90 index e04272c0..a1ca416e 100644 --- a/src/trixi.f90 +++ b/src/trixi.f90 @@ -163,8 +163,8 @@ integer(c_int) function trixi_nvariables(handle) bind(c) !! !! @see @ref trixi_get_cell_averages_api_c "trixi_get_cell_averages (C API)" subroutine trixi_get_cell_averages(data, handle) bind(c) - use, intrinsic :: iso_c_binding, only: c_int, c_double - real(c_double), intent(in) :: data(*) + use, intrinsic :: iso_c_binding, only: c_int, c_ptr + type(c_ptr), value :: data integer(c_int), value, intent(in) :: handle end subroutine From 5db2864cb8439a8d1f40aefece76c381647e8bed Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 16:16:01 +0200 Subject: [PATCH 09/25] use load instead of get --- LibTrixi.jl/src/LibTrixi.jl | 6 +++--- LibTrixi.jl/src/api_c.jl | 17 +++++++++++------ LibTrixi.jl/src/api_jl.jl | 2 +- examples/trixi_controller_data.c | 2 +- examples/trixi_controller_data.f90 | 2 +- src/trixi.c | 18 ++++++++++-------- src/trixi.f90 | 6 +++--- src/trixi.h | 2 +- 8 files changed, 31 insertions(+), 24 deletions(-) diff --git a/LibTrixi.jl/src/LibTrixi.jl b/LibTrixi.jl/src/LibTrixi.jl index 852de171..5c01200b 100644 --- a/LibTrixi.jl/src/LibTrixi.jl +++ b/LibTrixi.jl/src/LibTrixi.jl @@ -27,9 +27,9 @@ export trixi_nelements, export trixi_nvariables, trixi_nvariables_cfptr, trixi_nvariables_jl -export trixi_get_cell_averages, - trixi_get_cell_averages_cfptr, - trixi_get_cell_averages_jl +export trixi_load_cell_averages, + trixi_load_cell_averages_cfptr, + trixi_load_cell_averages_jl export SimulationState, store_simstate, load_simstate, delete_simstate! diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index afe88d91..5941b6d3 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -177,16 +177,21 @@ trixi_nvariables_cfptr() = @cfunction(trixi_nvariables, Cint, (Cint,)) """ - trixi_get_cell_averages(simstate_handle::Cint)::Cint + trixi_load_cell_averages(simstate_handle::Cint)::Cint -Return number of elements (cells) +Return cell averaged values + +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. +The given array has to be of correct size and memory has to be allocated beforehand. """ -function trixi_get_cell_averages end +function trixi_load_cell_averages end -Base.@ccallable function trixi_get_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid +Base.@ccallable function trixi_load_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid simstate = load_simstate(simstate_handle) - trixi_get_cell_averages_jl(data, simstate) + trixi_load_cell_averages_jl(data, simstate) return nothing end -trixi_get_cell_averages_cfptr() = @cfunction(trixi_get_cell_averages, Cvoid, (Ptr{Cdouble}, Cint,)) \ No newline at end of file +trixi_load_cell_averages_cfptr() = @cfunction(trixi_load_cell_averages, Cvoid, (Ptr{Cdouble}, Cint,)) \ No newline at end of file diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index 6de8dcea..aa759c28 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -65,7 +65,7 @@ function trixi_nvariables_jl(simstate) return Trixi.nvariables(equations) end -function trixi_get_cell_averages_jl(data_, simstate) +function trixi_load_cell_averages_jl(data_, simstate) mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(simstate.semi) nelements = Trixi.nelements(solver, cache) diff --git a/examples/trixi_controller_data.c b/examples/trixi_controller_data.c index 3f120faf..7436a599 100644 --- a/examples/trixi_controller_data.c +++ b/examples/trixi_controller_data.c @@ -43,7 +43,7 @@ int main ( int argc, char *argv[] ) { double* data = malloc( sizeof(double) * nelements * nvariables ); // get averaged cell values for each variable - trixi_get_cell_averages(data, handle); + trixi_load_cell_averages(data, handle); // compute temperature const double gas_constant = 0.287; diff --git a/examples/trixi_controller_data.f90 b/examples/trixi_controller_data.f90 index 284d6fe0..7b5c7a38 100644 --- a/examples/trixi_controller_data.f90 +++ b/examples/trixi_controller_data.f90 @@ -57,7 +57,7 @@ program simple_trixi_controller_f allocate ( data(0:nelements*nvariables) ) ! get averaged cell values for each variable - call trixi_get_cell_averages(c_loc(data),handle); + call trixi_load_cell_averages(c_loc(data),handle); ! compute temperature gas_constant = 0.287; diff --git a/src/trixi.c b/src/trixi.c index 0e0b1046..bfd348e2 100644 --- a/src/trixi.c +++ b/src/trixi.c @@ -30,7 +30,7 @@ enum { TRIXI_FTPR_NDIMS, TRIXI_FTPR_NELEMENTS, TRIXI_FTPR_NVARIABLES, - TRIXI_FTPR_GET_CELL_AVERAGES, + TRIXI_FTPR_LOAD_CELL_AVERAGES, // The last one is for the array size TRIXI_NUM_FPTRS @@ -48,7 +48,7 @@ static const char* trixi_function_pointer_names[] = { [TRIXI_FTPR_NDIMS] = "trixi_ndims_cfptr", [TRIXI_FTPR_NELEMENTS] = "trixi_nelements_cfptr", [TRIXI_FTPR_NVARIABLES] = "trixi_nvariables_cfptr", - [TRIXI_FTPR_GET_CELL_AVERAGES] = "trixi_get_cell_averages_cfptr", + [TRIXI_FTPR_LOAD_CELL_AVERAGES] = "trixi_load_cell_averages_cfptr", }; // Default depot path *relative* to the project directory @@ -300,23 +300,25 @@ int trixi_nvariables(int handle) { /** - * @anchor trixi_get_cell_averages_api_c + * @anchor trixi_load_cell_averages_api_c * * @brief Return cell averaged values * - * Cell averaged values for each cell and each variable are stored in a contiguous array. + * 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. * The given array has to be of correct size and memory has to be allocated beforehand. * * @param[in] handle simulation handle - * @param[out] data cell averaged values for all cells and all variables + * @param[out] data cell averaged values for all cells and all primitive variables */ -void trixi_get_cell_averages(double * data, int handle) { +void trixi_load_cell_averages(double * data, int handle) { // Get function pointer - void (*get_cell_averages)(double *, int) = trixi_function_pointers[TRIXI_FTPR_GET_CELL_AVERAGES]; + void (*load_cell_averages)(double *, int) = trixi_function_pointers[TRIXI_FTPR_LOAD_CELL_AVERAGES]; // Call function - get_cell_averages(data, handle); + load_cell_averages(data, handle); } diff --git a/src/trixi.f90 b/src/trixi.f90 index a1ca416e..f83803ac 100644 --- a/src/trixi.f90 +++ b/src/trixi.f90 @@ -154,15 +154,15 @@ integer(c_int) function trixi_nvariables(handle) bind(c) end function !> - !! @fn LibTrixi::trixi_get_cell_averages::trixi_get_cell_averages(data, handle) + !! @fn LibTrixi::trixi_load_cell_averages::trixi_load_cell_averages(data, handle) !! !! @brief Return cell averaged values !! !! @param[in] handle simulation handle !! @param[out] data cell averaged values for all cells and all variables !! - !! @see @ref trixi_get_cell_averages_api_c "trixi_get_cell_averages (C API)" - subroutine trixi_get_cell_averages(data, handle) bind(c) + !! @see @ref trixi_load_cell_averages_api_c "trixi_load_cell_averages (C API)" + subroutine trixi_load_cell_averages(data, handle) bind(c) use, intrinsic :: iso_c_binding, only: c_int, c_ptr type(c_ptr), value :: data integer(c_int), value, intent(in) :: handle diff --git a/src/trixi.h b/src/trixi.h index d30bc523..fcb6a577 100644 --- a/src/trixi.h +++ b/src/trixi.h @@ -22,7 +22,7 @@ int trixi_nelements(int handle); // Return number of elements (cells) int trixi_nvariables(int handle); // Return number of (conservative) variables // Data -void trixi_get_cell_averages(double * data, int handle); +void trixi_load_cell_averages(double * data, int handle); // Misc double trixi_calculate_dt(int handle); From 874f710a0c6bc309e19e186386b597ede1290e1e Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Fri, 14 Jul 2023 16:16:54 +0200 Subject: [PATCH 10/25] Update LibTrixi.jl/src/api_c.jl Co-authored-by: Michael Schlottke-Lakemper --- LibTrixi.jl/src/api_c.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index 5941b6d3..3a9afe57 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -149,7 +149,7 @@ trixi_ndims_cfptr() = @cfunction(trixi_ndims, Cint, (Cint,)) """ trixi_nelements(simstate_handle::Cint)::Cint -Return number of elements (cells) +Return number of elements (cells). """ function trixi_nelements end From 9725e2a9452252a48b8296b37a76f62bd510ed1a Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Fri, 14 Jul 2023 16:17:06 +0200 Subject: [PATCH 11/25] Update LibTrixi.jl/src/api_c.jl Co-authored-by: Michael Schlottke-Lakemper --- LibTrixi.jl/src/api_c.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index 3a9afe57..a7d9a845 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -134,7 +134,7 @@ trixi_step_cfptr() = @cfunction(trixi_step, Cvoid, (Cint,)) """ trixi_ndims(simstate_handle::Cint)::Cint -Return number of spatial dimensions +Return number of spatial dimensions. """ function trixi_ndims end From cddfd8ba5419cd881e3fa9912faf15de925c06ac Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Fri, 14 Jul 2023 16:17:30 +0200 Subject: [PATCH 12/25] Update examples/trixi_controller_data.c Co-authored-by: Michael Schlottke-Lakemper --- examples/trixi_controller_data.c | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/trixi_controller_data.c b/examples/trixi_controller_data.c index 7436a599..87f705f6 100644 --- a/examples/trixi_controller_data.c +++ b/examples/trixi_controller_data.c @@ -49,7 +49,6 @@ int main ( int argc, char *argv[] ) { 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]) ); } From ab7bd6efe037a88c75a8b67235611d064b2b0e89 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Fri, 14 Jul 2023 16:19:25 +0200 Subject: [PATCH 13/25] Update LibTrixi.jl/src/api_jl.jl Co-authored-by: Michael Schlottke-Lakemper --- LibTrixi.jl/src/api_jl.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index aa759c28..90f5920f 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -66,7 +66,6 @@ function trixi_nvariables_jl(simstate) end function trixi_load_cell_averages_jl(data_, simstate) - mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(simstate.semi) nelements = Trixi.nelements(solver, cache) nvariables = Trixi.nvariables(equations) From d027822f0bed6a59cc67bbbbd92fcaeccb45d298 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 16:24:50 +0200 Subject: [PATCH 14/25] small fixes --- LibTrixi.jl/src/api_c.jl | 2 +- examples/trixi_controller_data.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index a7d9a845..876f3bdb 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -164,7 +164,7 @@ trixi_nelements_cfptr() = @cfunction(trixi_nelements, Cint, (Cint,)) """ trixi_nvariables(simstate_handle::Cint)::Cint -Return number of elements (cells) +Return number of variables. """ function trixi_nvariables end diff --git a/examples/trixi_controller_data.f90 b/examples/trixi_controller_data.f90 index 7b5c7a38..cbc71a9d 100644 --- a/examples/trixi_controller_data.f90 +++ b/examples/trixi_controller_data.f90 @@ -57,7 +57,7 @@ program simple_trixi_controller_f allocate ( data(0:nelements*nvariables) ) ! get averaged cell values for each variable - call trixi_load_cell_averages(c_loc(data),handle); + call trixi_load_cell_averages(c_loc(data), handle); ! compute temperature gas_constant = 0.287; From 4dc4cb5c418220ec25e0bcffb4999fd3da66e33a Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Fri, 14 Jul 2023 17:18:16 +0200 Subject: [PATCH 15/25] add all required Trixi symbols to using statement --- LibTrixi.jl/src/LibTrixi.jl | 3 ++- LibTrixi.jl/src/api_jl.jl | 36 +++++++++++++++++++----------------- 2 files changed, 21 insertions(+), 18 deletions(-) diff --git a/LibTrixi.jl/src/LibTrixi.jl b/LibTrixi.jl/src/LibTrixi.jl index 5c01200b..a8efec25 100644 --- a/LibTrixi.jl/src/LibTrixi.jl +++ b/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, + wrap_array, eachelement, cons2prim, get_node_vars, eachnode export trixi_initialize_simulation, trixi_initialize_simulation_cfptr, diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index 90f5920f..41fce557 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -51,39 +51,41 @@ function trixi_step_jl(simstate) end function trixi_ndims_jl(simstate) - mesh, _, _, _ = Trixi.mesh_equations_solver_cache(simstate.semi) + mesh, _, _, _ = mesh_equations_solver_cache(simstate.semi) return ndims(mesh) end function trixi_nelements_jl(simstate) - _, _, solver, cache = Trixi.mesh_equations_solver_cache(simstate.semi) - return Trixi.nelements(solver, cache) + _, _, solver, cache = mesh_equations_solver_cache(simstate.semi) + return nelements(solver, cache) end function trixi_nvariables_jl(simstate) - _, equations, _, api_jl = Trixi.mesh_equations_solver_cache(simstate.semi) - return Trixi.nvariables(equations) + _, equations, _, _ = mesh_equations_solver_cache(simstate.semi) + return nvariables(equations) end function trixi_load_cell_averages_jl(data_, simstate) - mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(simstate.semi) - nelements = Trixi.nelements(solver, cache) - nvariables = Trixi.nvariables(equations) + mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi) + nels = nelements(solver, cache) + nvars = nvariables(equations) # convert C to julia - data = unsafe_wrap(Array, data_, nelements*nvariables) + data = unsafe_wrap(Array, data_, nels*nvars) u_ode = simstate.integrator.u - u = Trixi.wrap_array(u_ode, mesh, equations, solver, cache) + u = wrap_array(u_ode, mesh, equations, solver, cache) - for element in Trixi.eachelement(solver, cache) + # temporary storage for mean value on current element for all variables + u_mean = zero(get_node_vars(u, equations, solver, 1, 1, 1)) - # temporary storage for mean value on current element for all variables - u_mean = zero(Trixi.get_node_vars(u, equations, solver, 1, 1, element)) + for element in eachelement(solver, cache) + + u_mean = zero(u_mean) # compute mean value using nodal dg values and quadrature - for j in Trixi.eachnode(solver), i in Trixi.eachnode(solver) - u_node_prim = Trixi.cons2prim(Trixi.get_node_vars(u, equations, solver, i, j, element), equations) + for j in eachnode(solver), i in eachnode(solver) + u_node_prim = cons2prim(get_node_vars(u, equations, solver, i, j, element), equations) u_mean += u_node_prim * solver.basis.weights[i] * solver.basis.weights[j] end @@ -91,8 +93,8 @@ function trixi_load_cell_averages_jl(data_, simstate) u_mean = u_mean / 2^ndims(mesh) # copy to provided array - for ivar = 0:nvariables-1 - data[element + ivar * nelements] = u_mean[ivar+1] + for ivar = 0:nvars-1 + data[element + ivar * nels] = u_mean[ivar+1] end end From 41c3bc4f9bf6d5b982fa160d92769d2022e3bf03 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Thu, 3 Aug 2023 21:13:52 +0200 Subject: [PATCH 16/25] changed data interface to real(c_double) --- examples/trixi_controller_data.f90 | 7 ++++--- src/trixi.f90 | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/trixi_controller_data.f90 b/examples/trixi_controller_data.f90 index cbc71a9d..a86ea797 100644 --- a/examples/trixi_controller_data.f90 +++ b/examples/trixi_controller_data.f90 @@ -1,13 +1,14 @@ program simple_trixi_controller_f use LibTrixi use, intrinsic :: iso_fortran_env, only: error_unit - use, intrinsic :: iso_c_binding, only: c_int, c_null_char, c_double, c_loc + use, intrinsic :: iso_c_binding, only: c_int, c_double implicit none integer(c_int) :: handle, nelements, nvariables, i character(len=256) :: argument - real(c_double), allocatable, target :: data(:) + integer, parameter :: dp = selected_real_kind(12) + real(dp), dimension(:), pointer :: data real(c_double) :: gas_constant @@ -57,7 +58,7 @@ program simple_trixi_controller_f allocate ( data(0:nelements*nvariables) ) ! get averaged cell values for each variable - call trixi_load_cell_averages(c_loc(data), handle); + call trixi_load_cell_averages(data, handle); ! compute temperature gas_constant = 0.287; diff --git a/src/trixi.f90 b/src/trixi.f90 index c4d4090e..f88c3933 100644 --- a/src/trixi.f90 +++ b/src/trixi.f90 @@ -163,8 +163,8 @@ integer(c_int) function trixi_nvariables(handle) bind(c) !! !! @see @ref trixi_load_cell_averages_api_c "trixi_load_cell_averages (C API)" subroutine trixi_load_cell_averages(data, handle) bind(c) - use, intrinsic :: iso_c_binding, only: c_int, c_ptr - type(c_ptr), value :: data + use, intrinsic :: iso_c_binding, only: c_int, c_double + real(c_double), dimension(*), intent(in) :: data integer(c_int), value, intent(in) :: handle end subroutine From a39b3c90efaeddd0c8f86a20fbaebc0f30affe67 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Aug 2023 10:03:12 +0200 Subject: [PATCH 17/25] make computation of averages dimension agnostic --- LibTrixi.jl/src/LibTrixi.jl | 2 +- LibTrixi.jl/src/api_jl.jl | 33 +++++++++++++++++++++------------ 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/LibTrixi.jl/src/LibTrixi.jl b/LibTrixi.jl/src/LibTrixi.jl index f4e73895..5bdb6791 100644 --- a/LibTrixi.jl/src/LibTrixi.jl +++ b/LibTrixi.jl/src/LibTrixi.jl @@ -2,7 +2,7 @@ module LibTrixi using OrdinaryDiffEq: OrdinaryDiffEq, step!, check_error, DiscreteCallback using Trixi: Trixi, summary_callback, mesh_equations_solver_cache, nelements, nvariables, - wrap_array, eachelement, cons2prim, get_node_vars, eachnode + nnodes, wrap_array, eachelement, cons2prim, get_node_vars, eachnode using MPI: MPI, run_init_hooks, set_default_error_handler_return export trixi_initialize_simulation, diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index 41fce557..cd4e86e6 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -67,34 +67,43 @@ end function trixi_load_cell_averages_jl(data_, simstate) mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi) - nels = nelements(solver, cache) - nvars = nvariables(equations) + n_elements = nelements(solver, cache) + n_variables = nvariables(equations) + n_nodes = nnodes(solver) + n_dims = ndims(mesh) # convert C to julia - data = unsafe_wrap(Array, data_, nels*nvars) + data = unsafe_wrap(Array, data_, n_elements*n_variables) 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 = zero(get_node_vars(u, equations, solver, 1, 1, 1)) + u_mean = get_node_vars(u, equations, solver, node_cis[1], 1) for element in eachelement(solver, cache) - u_mean = zero(u_mean) - # compute mean value using nodal dg values and quadrature - for j in eachnode(solver), i in eachnode(solver) - u_node_prim = cons2prim(get_node_vars(u, equations, solver, i, j, element), equations) - u_mean += u_node_prim * solver.basis.weights[i] * solver.basis.weights[j] + 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^ndims(mesh) + u_mean = u_mean / 2^n_dims # copy to provided array - for ivar = 0:nvars-1 - data[element + ivar * nels] = u_mean[ivar+1] + # 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 From 0f853f6e22917f37cf5d575ba6be58ee1f76499e Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Mon, 14 Aug 2023 10:06:32 +0200 Subject: [PATCH 18/25] Update LibTrixi.jl/src/api_c.jl Co-authored-by: Michael Schlottke-Lakemper --- LibTrixi.jl/src/api_c.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index 876f3bdb..5853d33f 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -177,7 +177,7 @@ trixi_nvariables_cfptr() = @cfunction(trixi_nvariables, Cint, (Cint,)) """ - trixi_load_cell_averages(simstate_handle::Cint)::Cint + trixi_load_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid Return cell averaged values From 63152971efe9b51eebef4b02fccd4c175f0687b0 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Mon, 14 Aug 2023 10:07:11 +0200 Subject: [PATCH 19/25] Update LibTrixi.jl/src/api_c.jl Co-authored-by: Michael Schlottke-Lakemper --- LibTrixi.jl/src/api_c.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index 5853d33f..dc3b826b 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -179,7 +179,7 @@ trixi_nvariables_cfptr() = @cfunction(trixi_nvariables, Cint, (Cint,)) """ trixi_load_cell_averages(data::Ptr{Cdouble}, simstate_handle::Cint)::Cvoid -Return cell averaged values +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 From 48aae60225235ca41d2c62325561cd19d139c6b7 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Mon, 14 Aug 2023 10:08:28 +0200 Subject: [PATCH 20/25] Update LibTrixi.jl/src/api_c.jl Co-authored-by: Michael Schlottke-Lakemper --- LibTrixi.jl/src/api_c.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index dc3b826b..f2c600cc 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -183,7 +183,8 @@ 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. +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 From ad0d308c997baff1047814636d39cfa718419ff3 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Aug 2023 10:32:12 +0200 Subject: [PATCH 21/25] doc consistent to suggestion --- src/trixi.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/trixi.c b/src/trixi.c index bfd348e2..2b41a4ad 100644 --- a/src/trixi.c +++ b/src/trixi.c @@ -306,7 +306,8 @@ int trixi_nvariables(int handle) { * * 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. + * the other variables subsequently (structure-of-arrays layout). + * * The given array has to be of correct size and memory has to be allocated beforehand. * * @param[in] handle simulation handle From 5d63a5320673d8ccecb6e1800233ae7f9a0f978d Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Aug 2023 10:32:46 +0200 Subject: [PATCH 22/25] move conversion of C to julia array from api_jl to api_c --- LibTrixi.jl/src/api_c.jl | 7 ++++++- LibTrixi.jl/src/api_jl.jl | 5 +---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index f2c600cc..a9bbbcbb 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -191,7 +191,12 @@ 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) - trixi_load_cell_averages_jl(data, simstate) + + # 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 diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index cd4e86e6..76f34b8b 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -65,16 +65,13 @@ function trixi_nvariables_jl(simstate) return nvariables(equations) end -function trixi_load_cell_averages_jl(data_, simstate) +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) - # convert C to julia - data = unsafe_wrap(Array, data_, n_elements*n_variables) - u_ode = simstate.integrator.u u = wrap_array(u_ode, mesh, equations, solver, cache) From 8a3b79397f2fea7b1227a758cca13e9f6113d2ec Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Aug 2023 10:45:57 +0200 Subject: [PATCH 23/25] switch to 1-based Fortran array --- examples/trixi_controller_data.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/trixi_controller_data.f90 b/examples/trixi_controller_data.f90 index a86ea797..a28ca003 100644 --- a/examples/trixi_controller_data.f90 +++ b/examples/trixi_controller_data.f90 @@ -55,7 +55,7 @@ program simple_trixi_controller_f write(*, '(a)') "" ! allocate memory - allocate ( data(0:nelements*nvariables) ) + allocate ( data(nelements*nvariables) ) ! get averaged cell values for each variable call trixi_load_cell_averages(data, handle); @@ -63,7 +63,7 @@ program simple_trixi_controller_f ! compute temperature gas_constant = 0.287; - do i = 0,nelements-1 + do i = 1,nelements print "('T[cell ', i4, '] = ', e14.8)", i, data(i+3*nelements)/(gas_constant*data(i)) end do From d71c53c30c0c8e971ef993fc13a84d5d42cb6d23 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Aug 2023 10:50:41 +0200 Subject: [PATCH 24/25] removed example, libelixir_p4est2d_dgsem_euler_sedov.jl can be used --- .../examples/libelixir_tree2d_dgsem_euler.jl | 74 ------------------- 1 file changed, 74 deletions(-) delete mode 100644 LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl diff --git a/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl b/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl deleted file mode 100644 index 92677e73..00000000 --- a/LibTrixi.jl/examples/libelixir_tree2d_dgsem_euler.jl +++ /dev/null @@ -1,74 +0,0 @@ - -using LibTrixi -using Trixi -using OrdinaryDiffEq - -# The function to create the simulation state needs to be named `init_simstate` -function init_simstate() - - ############################################################################### - # semidiscretization of the compressible Euler equations - - equations = CompressibleEulerEquations2D(1.4) - - initial_condition = initial_condition_convergence_test - - source_terms = source_terms_convergence_test - - boundary_condition = BoundaryConditionDirichlet(initial_condition) - - boundary_conditions = (x_neg=boundary_condition, - x_pos=boundary_condition, - y_neg=boundary_condition, - y_pos=boundary_condition,) - - solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) - - coordinates_min = (0.0, 0.0) - coordinates_max = (2.0, 2.0) - - mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=2, - n_cells_max=10_000) - - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms=source_terms, - boundary_conditions=boundary_conditions) - - - ############################################################################### - # ODE solvers, callbacks etc. - - tspan = (0.0, 1.0) - ode = semidiscretize(semi, tspan) - - summary_callback = SummaryCallback() - - analysis_interval = 100 - #analysis_callback = AnalysisCallback(semi, interval=analysis_interval) - - alive_callback = AliveCallback(analysis_interval=analysis_interval) - - stepsize_callback = StepsizeCallback(cfl=0.8) - - callbacks = CallbackSet(summary_callback, - alive_callback, - stepsize_callback) - - - ############################################################################### - # create the time integrator - - integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false), - dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep=false, callback=callbacks); - - - ############################################################################### - # Create simulation state - - simstate = SimulationState(semi, integrator) - - return simstate - -end From 1af00986ba38ff4ec6530e48bfcd22d1ed57dab5 Mon Sep 17 00:00:00 2001 From: Benedict Geihe Date: Mon, 14 Aug 2023 11:48:04 +0200 Subject: [PATCH 25/25] add new examples to CI testing --- .github/workflows/ci.yml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 142f5d2e..0e8bb708 100644 --- a/.github/workflows/ci.yml +++ b/.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: |