diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a4af483b..b8cfac1a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -243,10 +243,12 @@ jobs: # all controllers ../build/examples/trixi_controller_simple_c . libelixir_tree1d_advection_basic.jl ../build/examples/trixi_controller_simple_f . libelixir_tree1d_advection_basic.jl - ../build/examples/trixi_controller_data_c . libelixir_t8code2d_advection_amr.jl - ../build/examples/trixi_controller_data_f . libelixir_t8code2d_advection_amr.jl - ../build/examples/trixi_controller_t8code_c . libelixir_t8code2d_advection_amr.jl - ../build/examples/trixi_controller_t8code_f . libelixir_t8code2d_advection_amr.jl + ../build/examples/trixi_controller_data_load_c . libelixir_t8code2d_euler_tracer_amr.jl + ../build/examples/trixi_controller_data_load_f . libelixir_t8code2d_euler_tracer_amr.jl + ../build/examples/trixi_controller_data_store_c . libelixir_t8code2d_euler_tracer_amr.jl + ../build/examples/trixi_controller_data_store_f . libelixir_t8code2d_euler_tracer_amr.jl + ../build/examples/trixi_controller_t8code_c . libelixir_t8code2d_euler_tracer_amr.jl + ../build/examples/trixi_controller_t8code_f . libelixir_t8code2d_euler_tracer_amr.jl ../build/examples/trixi_controller_baroclinic_c . libelixir_t8code3d_euler_baroclinic_instability.jl ../build/examples/trixi_controller_baroclinic_f . libelixir_t8code3d_euler_baroclinic_instability.jl mpirun -n 2 ../build/examples/trixi_controller_mpi_c . libelixir_p4est2d_euler_sedov.jl @@ -275,10 +277,14 @@ jobs: "../build/examples/trixi_controller_simple_c ." \ "../build/examples/trixi_controller_simple_f" \ "../build/examples/trixi_controller_simple_f ." \ - "../build/examples/trixi_controller_data_c" \ - "../build/examples/trixi_controller_data_c ." \ - "../build/examples/trixi_controller_data_f" \ - "../build/examples/trixi_controller_data_f ." \ + "../build/examples/trixi_controller_data_load_c" \ + "../build/examples/trixi_controller_data_load_c ." \ + "../build/examples/trixi_controller_data_load_f" \ + "../build/examples/trixi_controller_data_load_f ." \ + "../build/examples/trixi_controller_data_store_c" \ + "../build/examples/trixi_controller_data_store_c ." \ + "../build/examples/trixi_controller_data_store_f" \ + "../build/examples/trixi_controller_data_store_f ." \ "../build/examples/trixi_controller_t8code_c" \ "../build/examples/trixi_controller_t8code_c ." \ "../build/examples/trixi_controller_t8code_f" \ diff --git a/LibTrixi.jl/Project.toml b/LibTrixi.jl/Project.toml index 99c61ccb..f11affa5 100644 --- a/LibTrixi.jl/Project.toml +++ b/LibTrixi.jl/Project.toml @@ -13,7 +13,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" MPI = "0.20.13" OrdinaryDiffEq = "6.53.2" Pkg = "1.8" -Trixi = "0.9.12, 0.10, 0.11, 0.12, 0.13" +Trixi = "0.11.12, 0.12, 0.13" julia = "1.8" [preferences.OrdinaryDiffEq] diff --git a/LibTrixi.jl/examples/libelixir_t8code2d_advection_amr.jl b/LibTrixi.jl/examples/libelixir_t8code2d_advection_amr.jl deleted file mode 100644 index a30080b7..00000000 --- a/LibTrixi.jl/examples/libelixir_t8code2d_advection_amr.jl +++ /dev/null @@ -1,83 +0,0 @@ -using LibTrixi -using OrdinaryDiffEq -using Trixi - -# The function to create the simulation state needs to be named `init_simstate` -function init_simstate() - - ############################################################################### - # semidiscretization of the linear advection equation - - advection_velocity = (0.2, -0.7) - equations = LinearScalarAdvectionEquation2D(advection_velocity) - - # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux - solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) - - coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) - coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) - - mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - - trees_per_dimension = (2, 2) - - mesh = T8codeMesh(trees_per_dimension, polydeg=3, - mapping=mapping, - initial_refinement_level=1) - - # A semidiscretization collects data structures and functions for the spatial discretization - semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) - - - ############################################################################### - # ODE solvers, callbacks etc. - - # Create ODE problem with time span from 0.0 to 0.2 - ode = semidiscretize(semi, (0.0, 0.2)); - - # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup - # and resets the timers - summary_callback = SummaryCallback() - - # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results - analysis_interval = 100 - analysis_callback = AnalysisCallback(semi, interval=analysis_interval) - - alive_callback = AliveCallback(analysis_interval=analysis_interval) - - # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step - stepsize_callback = StepsizeCallback(cfl=0.5) - - # 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) - amr_callback = AMRCallback(semi, amr_controller, - interval=10, - adapt_initial_condition=true, - adapt_initial_condition_only_refine=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, - stepsize_callback) - - - ############################################################################### - # create the time integrator - - # OrdinaryDiffEq's `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/LibTrixi.jl/examples/libelixir_t8code2d_euler_tracer_amr.jl b/LibTrixi.jl/examples/libelixir_t8code2d_euler_tracer_amr.jl new file mode 100644 index 00000000..2d927d05 --- /dev/null +++ b/LibTrixi.jl/examples/libelixir_t8code2d_euler_tracer_amr.jl @@ -0,0 +1,109 @@ +using LibTrixi +using OrdinaryDiffEq +using Trixi + +# The function to create the simulation state needs to be named `init_simstate` +function init_simstate() + + ############################################################################### + # initial condition: density wave + tracer blob + + function initial_condition_wave_blob(x, t, + equations::PassiveTracerEquations) + # Initial condition for flow equations: density wave + v1 = 0.1 + v2 = 0.2 + rho = 1 + 0.5 * sinpi(2 * (x[1] + x[2] - t * (v1 + v2))) + rho_v1 = rho * v1 + rho_v2 = rho * v2 + p = 20 + rho_e = p / (equations.flow_equations.gamma - 1) + 0.5 * rho * (v1^2 + v2^2) + + # Initial condition for tracers: blob in fraction of density + tracer = 0.2 * exp(-20 * (x[1] + 0.45)^2 - 10 * (x[2] - 0.15)^2) + + return SVector(rho, rho_v1, rho_v2, rho_e, rho * tracer) + end + + ############################################################################### + # semidiscretization of the compressible Euler equations + + gamma = 1.4 + flow_equations = CompressibleEulerEquations2D(gamma) + equations = PassiveTracerEquations(flow_equations, n_tracers = 1) + + # Create DG solver with polynomial degree = 5, Ranocha flux, and derived tracer flux + solver = DGSEM(polydeg = 5, surface_flux = FluxTracerEquationsCentral(flux_ranocha)) + + coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) + coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y)) + trees_per_dimension = (4, 4) # initial resolution (without refinement) + mesh = T8codeMesh(trees_per_dimension, polydeg = 1, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 1) + + # Create spatial discretization + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_wave_blob, + solver) + + ############################################################################### + # ODE solvers, callbacks etc. + + # Create ODE problem with time span from 0.0 to 2.0 + ode = semidiscretize(semi, (0.0, 2.0)); + + # SummaryCallback prints a summary of the simulation setup and recorded performance data + summary_callback = SummaryCallback() + + # AnalysisCallback analyses the solution at regular intervals and prints the results + analysis_interval = 100 + analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + + # AliveCallback prints a one line summary at regular intervals + alive_callback = AliveCallback(analysis_interval=analysis_interval) + + # StepsizeCallback handles the re-calculation of the maximum Δt after each time step + stepsize_callback = StepsizeCallback(cfl = 1.0) + + # AMRCallback triggers adaptive mesh refinement + @inline function first_tracer(u, equations::PassiveTracerEquations) + return Trixi.tracers(u, equations)[1] + end + amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first_tracer), + base_level=1, + med_level=2, med_threshold=0.05, + max_level=2, max_threshold=0.05) + amr_callback = AMRCallback(semi, amr_controller, + interval=50, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + + # SaveSolutionCallback writes the solution at regular intervals + save_solution = SaveSolutionCallback(interval=50, + save_initial_solution=true, + save_final_solution=true, + solution_variables = cons2prim) + + # CallbackSet collects all callbacks + callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + save_solution, + stepsize_callback) + + ############################################################################### + # create the time integrator + + # OrdinaryDiffEq's `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/LibTrixi.jl/src/LibTrixi.jl b/LibTrixi.jl/src/LibTrixi.jl index 0a84dfce..0c1ae24d 100644 --- a/LibTrixi.jl/src/LibTrixi.jl +++ b/LibTrixi.jl/src/LibTrixi.jl @@ -52,15 +52,24 @@ export trixi_load_node_reference_coordinates, export trixi_load_node_weights, trixi_load_node_weights_cfptr, trixi_load_node_weights_jl -export trixi_load_primitive_vars, - trixi_load_primitive_vars_cfptr, - trixi_load_primitive_vars_jl -export trixi_load_element_averaged_primitive_vars, - trixi_load_element_averaged_primitive_vars_cfptr, - trixi_load_element_averaged_primitive_vars_jl +export trixi_load_conservative_var, + trixi_load_conservative_var_cfptr, + trixi_load_conservative_var_jl +export trixi_load_primitive_var, + trixi_load_primitive_var_cfptr, + trixi_load_primitive_var_jl +export trixi_load_element_averaged_primitive_var, + trixi_load_element_averaged_primitive_var_cfptr, + trixi_load_element_averaged_primitive_var_jl +export trixi_store_conservative_var, + trixi_store_conservative_var_cfptr, + trixi_store_conservative_var_jl export trixi_register_data, trixi_register_data_cfptr, trixi_register_data_jl +export trixi_get_conservative_vars_pointer, + trixi_get_conservative_vars_pointer_cfptr, + trixi_get_conservative_vars_pointer_jl export trixi_version_library, trixi_version_library_cfptr, trixi_version_library_jl diff --git a/LibTrixi.jl/src/api_c.jl b/LibTrixi.jl/src/api_c.jl index 68bb408d..d5629197 100644 --- a/LibTrixi.jl/src/api_c.jl +++ b/LibTrixi.jl/src/api_c.jl @@ -433,8 +433,32 @@ trixi_load_node_weights_cfptr() = """ - trixi_load_primitive_vars(simstate_handle::Cint, variable_id::Cint, - data::Ptr{Cdouble})::Cvoid + trixi_load_conservative_var(simstate_handle::Cint, variable_id::Cint, + data::Ptr{Cdouble})::Cvoid + +Load conservative variable. +""" +function trixi_load_conservative_var end + +Base.@ccallable function trixi_load_conservative_var(simstate_handle::Cint, + variable_id::Cint, + data::Ptr{Cdouble})::Cvoid + simstate = load_simstate(simstate_handle) + + # convert C to Julia array + size = trixi_ndofs_jl(simstate) + data_jl = unsafe_wrap(Array, data, size) + + trixi_load_conservative_var_jl(simstate, variable_id, data_jl) +end + +trixi_load_conservative_var_cfptr() = + @cfunction(trixi_load_conservative_var, Cvoid, (Cint, Cint, Ptr{Cdouble})) + + +""" + trixi_load_primitive_var(simstate_handle::Cint, variable_id::Cint, + data::Ptr{Cdouble})::Cvoid Load primitive variable. @@ -443,22 +467,46 @@ are stored in the given array `data`. The given array has to be of correct size (ndofs) and memory has to be allocated beforehand. """ -function trixi_load_primitive_vars end +function trixi_load_primitive_var end -Base.@ccallable function trixi_load_primitive_vars(simstate_handle::Cint, variable_id::Cint, - data::Ptr{Cdouble})::Cvoid +Base.@ccallable function trixi_load_primitive_var(simstate_handle::Cint, variable_id::Cint, + data::Ptr{Cdouble})::Cvoid simstate = load_simstate(simstate_handle) # convert C to Julia array size = trixi_ndofs_jl(simstate) data_jl = unsafe_wrap(Array, data, size) - trixi_load_primitive_vars_jl(simstate, variable_id, data_jl) + trixi_load_primitive_var_jl(simstate, variable_id, data_jl) return nothing end -trixi_load_primitive_vars_cfptr() = - @cfunction(trixi_load_primitive_vars, Cvoid, (Cint, Cint, Ptr{Cdouble})) +trixi_load_primitive_var_cfptr() = + @cfunction(trixi_load_primitive_var, Cvoid, (Cint, Cint, Ptr{Cdouble})) + + +""" + trixi_store_conservative_var(simstate_handle::Cint, variable_id::Cint, + data::Ptr{Cdouble})::Cvoid + +Store conservative variable. +""" +function trixi_store_conservative_var end + +Base.@ccallable function trixi_store_conservative_var(simstate_handle::Cint, + variable_id::Cint, + data::Ptr{Cdouble})::Cvoid + simstate = load_simstate(simstate_handle) + + # convert C to Julia array + size = trixi_ndofs_jl(simstate) + data_jl = unsafe_wrap(Array, data, size) + + trixi_store_conservative_var_jl(simstate, variable_id, data_jl) +end + +trixi_store_conservative_var_cfptr() = + @cfunction(trixi_store_conservative_var, Cvoid, (Cint, Cint, Ptr{Cdouble})) """ @@ -509,8 +557,8 @@ trixi_get_simulation_time_cfptr() = @cfunction(trixi_get_simulation_time, Cdoubl """ - trixi_load_element_averaged_primitive_vars(simstate_handle::Cint, variable_id::Cint, - data::Ptr{Cdouble})::Cvoid + trixi_load_element_averaged_primitive_var(simstate_handle::Cint, variable_id::Cint, + data::Ptr{Cdouble})::Cvoid Load element averages for primitive variable. @@ -520,9 +568,9 @@ element are stored in the given array `data`. The given array has to be of correct size (nelements) and memory has to be allocated beforehand. """ -function trixi_load_element_averaged_primitive_vars end +function trixi_load_element_averaged_primitive_var end -Base.@ccallable function trixi_load_element_averaged_primitive_vars(simstate_handle::Cint, +Base.@ccallable function trixi_load_element_averaged_primitive_var(simstate_handle::Cint, variable_id::Cint, data::Ptr{Cdouble})::Cvoid simstate = load_simstate(simstate_handle) @@ -530,14 +578,34 @@ Base.@ccallable function trixi_load_element_averaged_primitive_vars(simstate_han size = trixi_nelements_jl(simstate) data_jl = unsafe_wrap(Array, data, size) - trixi_load_element_averaged_primitive_vars_jl(simstate, variable_id, data_jl) + trixi_load_element_averaged_primitive_var_jl(simstate, variable_id, data_jl) return nothing end -trixi_load_element_averaged_primitive_vars_cfptr() = - @cfunction(trixi_load_element_averaged_primitive_vars, Cvoid, (Cint, Cint, Ptr{Cdouble})) +trixi_load_element_averaged_primitive_var_cfptr() = + @cfunction(trixi_load_element_averaged_primitive_var, Cvoid, (Cint, Cint, Ptr{Cdouble})) + +""" + trixi_get_conservative_vars_pointer(simstate_handle::Cint)::Ptr{Cdouble} +Return pointer to internal data vector. +""" +function trixi_get_conservative_vars_pointer end + +Base.@ccallable function trixi_get_conservative_vars_pointer(simstate_handle::Cint)::Ptr{Cdouble} + simstate = load_simstate(simstate_handle) + return trixi_get_conservative_vars_pointer_jl(simstate) +end + +trixi_get_conservative_vars_pointer_cfptr() = + @cfunction(trixi_get_conservative_vars_pointer, Ptr{Cdouble}, (Cint,)) + + + +############################################################################################ +# t8code +############################################################################################ """ trixi_get_t8code_forest(simstate_handle::Cint)::Ptr{Trixi.t8_forest} @@ -557,6 +625,8 @@ end trixi_get_t8code_forest_cfptr() = @cfunction(trixi_get_t8code_forest, Ptr{Trixi.t8_forest}, (Cint,)) + + ############################################################################################ # Auxiliary ############################################################################################ diff --git a/LibTrixi.jl/src/api_jl.jl b/LibTrixi.jl/src/api_jl.jl index ee9c354f..f60968ee 100644 --- a/LibTrixi.jl/src/api_jl.jl +++ b/LibTrixi.jl/src/api_jl.jl @@ -140,7 +140,32 @@ function trixi_load_node_weights_jl(simstate, data) end -function trixi_load_primitive_vars_jl(simstate, variable_id, data) +function trixi_load_conservative_var_jl(simstate, variable_id, data) + mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi) + n_nodes_per_dim = nnodes(solver) + n_dims = ndims(mesh) + n_nodes = n_nodes_per_dim^n_dims + + 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_per_dim, n_dims)) + node_lis = LinearIndices(node_cis) + + for element in eachelement(solver, cache) + for node_ci in node_cis + node_vars = get_node_vars(u, equations, solver, node_ci, element) + node_index = (element-1) * n_nodes + node_lis[node_ci] + data[node_index] = node_vars[variable_id] + end + end + + return nothing +end + + +function trixi_load_primitive_var_jl(simstate, variable_id, data) mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi) n_nodes_per_dim = nnodes(solver) n_dims = ndims(mesh) @@ -165,7 +190,7 @@ function trixi_load_primitive_vars_jl(simstate, variable_id, data) end -function trixi_load_element_averaged_primitive_vars_jl(simstate, variable_id, data) +function trixi_load_element_averaged_primitive_var_jl(simstate, variable_id, data) mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi) n_nodes = nnodes(solver) n_dims = ndims(mesh) @@ -201,6 +226,30 @@ function trixi_load_element_averaged_primitive_vars_jl(simstate, variable_id, da end +function trixi_store_conservative_var_jl(simstate, variable_id, data) + mesh, equations, solver, cache = mesh_equations_solver_cache(simstate.semi) + n_nodes_per_dim = nnodes(solver) + n_dims = ndims(mesh) + n_nodes = n_nodes_per_dim^n_dims + + 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_per_dim, n_dims)) + node_lis = LinearIndices(node_cis) + + for element in eachelement(solver, cache) + for node_ci in node_cis + node_index = (element-1) * n_nodes + node_lis[node_ci] + u[variable_id, node_ci, element] = data[node_index] + end + end + + return nothing +end + + function trixi_register_data_jl(simstate, index, data) simstate.registry[index] = data if show_debug_output() @@ -210,6 +259,11 @@ function trixi_register_data_jl(simstate, index, data) end +function trixi_get_conservative_vars_pointer_jl(simstate) + return pointer(simstate.integrator.u) +end + + function trixi_get_simulation_time_jl(simstate) return simstate.integrator.t end diff --git a/LibTrixi.jl/test/Project.toml b/LibTrixi.jl/test/Project.toml index 19feeb5d..3174d971 100644 --- a/LibTrixi.jl/test/Project.toml +++ b/LibTrixi.jl/test/Project.toml @@ -5,7 +5,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] OrdinaryDiffEq = "6.53.2" -Trixi = "0.9.12, 0.10, 0.11, 0.12, 0.13" +Trixi = "0.11.12, 0.12, 0.13" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false diff --git a/LibTrixi.jl/test/test_interface.jl b/LibTrixi.jl/test/test_interface.jl index 7278cbf0..97dc6bd3 100644 --- a/LibTrixi.jl/test/test_interface.jl +++ b/LibTrixi.jl/test/test_interface.jl @@ -136,17 +136,38 @@ end # compare element averaged values data_c = zeros(nelements_c) - trixi_load_element_averaged_primitive_vars(handle, Int32(1), pointer(data_c)) + trixi_load_element_averaged_primitive_var(handle, Int32(1), pointer(data_c)) data_jl = zeros(nelements_jl) - trixi_load_element_averaged_primitive_vars_jl(simstate_jl, 1, data_jl) + trixi_load_element_averaged_primitive_var_jl(simstate_jl, 1, data_jl) + @test data_c == data_jl + + # compare conservative variable values on all dofs + data_c = zeros(ndofs_c) + trixi_load_conservative_var(handle, Int32(1), pointer(data_c)) + data_jl = zeros(ndofs_jl) + trixi_load_conservative_var_jl(simstate_jl, 1, data_jl) @test data_c == data_jl # compare primitive variable values on all dofs data_c = zeros(ndofs_c) - trixi_load_primitive_vars(handle, Int32(1), pointer(data_c)) + trixi_load_primitive_var(handle, Int32(1), pointer(data_c)) data_jl = zeros(ndofs_jl) - trixi_load_primitive_vars_jl(simstate_jl, 1, data_jl) + trixi_load_primitive_var_jl(simstate_jl, 1, data_jl) @test data_c == data_jl + + # write 1.0 to first variable and compare via raw access + data_c = fill(1.0, ndofs_c) + trixi_store_conservative_var(handle, Int32(1), pointer(data_c)) + data_ptr_c = trixi_get_conservative_vars_pointer(handle) + data_jl = unsafe_wrap(Array, data_ptr_c, ndofs_c) + @test all(data_jl .== 1.0) + + # write 2.0 to first variable and compare via raw access + data_jl = fill(2.0, ndofs_jl) + trixi_store_conservative_var_jl(simstate_jl, 1, data_jl) + data_ptr_jl = trixi_get_conservative_vars_pointer_jl(simstate_jl) + data_jl = unsafe_wrap(Array, data_ptr_jl, ndofs_jl) + @test all(data_jl .== 2.0) end diff --git a/LibTrixi.jl/test/test_t8code.jl b/LibTrixi.jl/test/test_t8code.jl index 2f581595..dc69b8e1 100644 --- a/LibTrixi.jl/test/test_t8code.jl +++ b/LibTrixi.jl/test/test_t8code.jl @@ -6,7 +6,7 @@ using Trixi libelixir = joinpath(dirname(pathof(LibTrixi)), - "../examples/libelixir_t8code2d_advection_amr.jl") + "../examples/libelixir_t8code2d_euler_tracer_amr.jl") # initialize a simulation via API, receive a handle handle = trixi_initialize_simulation(libelixir) diff --git a/README.md b/README.md index 47735339..771d1ae0 100644 --- a/README.md +++ b/README.md @@ -213,7 +213,7 @@ aspects on how to use the C and Fortran APIs of libtrixi: - `trixi_controller_simple.(c|f90)`: basic usage - `trixi_controller_mpi.(c|f90)`: usage in the presence of MPI -- `trixi_controller_data.(c|f90)`: simulation data access +- `trixi_controller_data_(load|store).(c|f90)`: simulation data access - `trixi_controller_t8code.(c|f90)`: interacting with t8code If you just want to test the Julia part of libtrixi, i.e., LibTrixi.jl, you can also run diff --git a/docs/src/index.md b/docs/src/index.md index 3c19903d..62f1714c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -216,7 +216,7 @@ examples demonstrate different aspects on how to use the C and Fortran APIs of l - `trixi_controller_simple.(c|f90)`: basic usage - `trixi_controller_mpi.(c|f90)`: usage in the presence of MPI -- `trixi_controller_data.(c|f90)`: simulation data access +- `trixi_controller_data_(load|store).(c|f90)`: simulation data access - `trixi_controller_t8code.(c|f90)`: interacting with t8code If you just want to test the Julia part of libtrixi, i.e., LibTrixi.jl, you can also run diff --git a/examples/.gitignore b/examples/.gitignore index 70e92df7..f95b5fed 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -1,8 +1,12 @@ -trixi_controller_simple_c -trixi_controller_simple_f +trixi_controller_baroclinic_c +trixi_controller_baroclinic_f +trixi_controller_data_load_c +trixi_controller_data_load_f +trixi_controller_data_store_c +trixi_controller_data_store_f trixi_controller_mpi_c trixi_controller_mpi_f -trixi_controller_data_c -trixi_controller_data_f +trixi_controller_simple_c +trixi_controller_simple_f trixi_controller_t8code_c trixi_controller_t8code_f diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 1cdf293b..508fb4bf 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -3,8 +3,10 @@ set ( EXAMPLES trixi_controller_simple.f90 trixi_controller_mpi.c trixi_controller_mpi.f90 - trixi_controller_data.c - trixi_controller_data.f90 + trixi_controller_data_load.c + trixi_controller_data_load.f90 + trixi_controller_data_store.c + trixi_controller_data_store.f90 trixi_controller_t8code.c trixi_controller_t8code.f90 trixi_controller_baroclinic.c diff --git a/examples/trixi_controller_baroclinic.c b/examples/trixi_controller_baroclinic.c index cec6cc0e..22cdefe3 100644 --- a/examples/trixi_controller_baroclinic.c +++ b/examples/trixi_controller_baroclinic.c @@ -128,10 +128,10 @@ int main ( int argc, char *argv[] ) { while ( !trixi_is_finished( handle ) ) { // Get current state - trixi_load_primitive_vars( handle, 1, u1 ); - trixi_load_primitive_vars( handle, 2, u2 ); - trixi_load_primitive_vars( handle, 3, u3 ); - trixi_load_primitive_vars( handle, 4, u4 ); + trixi_load_primitive_var( handle, 1, u1 ); + trixi_load_primitive_var( handle, 2, u2 ); + trixi_load_primitive_var( handle, 3, u3 ); + trixi_load_primitive_var( handle, 4, u4 ); // Compute source terms source_terms_baroclinic( nnodes, nodes, forest, diff --git a/examples/trixi_controller_baroclinic.f90 b/examples/trixi_controller_baroclinic.f90 index d0b551b3..f0a95ff0 100644 --- a/examples/trixi_controller_baroclinic.f90 +++ b/examples/trixi_controller_baroclinic.f90 @@ -155,10 +155,10 @@ program trixi_controller_baroclinic_f if ( trixi_is_finished(handle) ) exit ! Get current state - call trixi_load_primitive_vars( handle, 1, u1 ) - call trixi_load_primitive_vars( handle, 2, u2 ) - call trixi_load_primitive_vars( handle, 3, u3 ) - call trixi_load_primitive_vars( handle, 4, u4 ) + call trixi_load_primitive_var( handle, 1, u1 ) + call trixi_load_primitive_var( handle, 2, u2 ) + call trixi_load_primitive_var( handle, 3, u3 ) + call trixi_load_primitive_var( handle, 4, u4 ) ! Compute source terms call source_terms_baroclinic( nnodes, nodes, forest, ndofs, & diff --git a/examples/trixi_controller_data.c b/examples/trixi_controller_data_load.c similarity index 96% rename from examples/trixi_controller_data.c rename to examples/trixi_controller_data_load.c index 35e88cdd..9ed41f3d 100644 --- a/examples/trixi_controller_data.c +++ b/examples/trixi_controller_data_load.c @@ -49,7 +49,7 @@ int main ( int argc, char *argv[] ) { data = realloc( data, sizeof(double) * nelements ); // Get element averaged values for first variable - trixi_load_element_averaged_primitive_vars(handle, 1, data); + trixi_load_element_averaged_primitive_var(handle, 1, data); } } diff --git a/examples/trixi_controller_data.f90 b/examples/trixi_controller_data_load.f90 similarity index 96% rename from examples/trixi_controller_data.f90 rename to examples/trixi_controller_data_load.f90 index d5d09667..b041d502 100644 --- a/examples/trixi_controller_data.f90 +++ b/examples/trixi_controller_data_load.f90 @@ -1,4 +1,4 @@ -program trixi_controller_data_f +program trixi_controller_data_load_f use LibTrixi use, intrinsic :: iso_fortran_env, only: error_unit use, intrinsic :: iso_c_binding, only: c_int, c_double @@ -66,7 +66,7 @@ program trixi_controller_data_f allocate( data(nelements) ) ! get element averaged values for first variable - call trixi_load_element_averaged_primitive_vars(handle, 1, data) + call trixi_load_element_averaged_primitive_var(handle, 1, data) end if end do diff --git a/examples/trixi_controller_data_store.c b/examples/trixi_controller_data_store.c new file mode 100644 index 00000000..7b28e515 --- /dev/null +++ b/examples/trixi_controller_data_store.c @@ -0,0 +1,95 @@ +#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], 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; + double* rho = NULL; + double* rho_tracer = NULL; + + printf("\n*** Trixi controller *** Entering main loop\n"); + while ( !trixi_is_finished( handle ) ) { + + trixi_step( handle ); + steps++; + + if (steps % 100 == 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_conservative_vars_pointer(handle); + + for (int i = 0; i < ndofs; ++i) { + // Density comes first + const double rho = raw_data[5*i]; + + // Tracer comes last + const double rho_tracer = raw_data[5*i + 4]; + + // Apply 10% damping to tracer (fraction of density) + const double tracer = 0.9 * (rho_tracer / rho); + raw_data[5*i + 4] = tracer * rho; + } + } + + if (steps % 100 == 50) { + + // Get number of degrees of freedom + int ndofs = trixi_ndofsglobal( handle ); + + // Allocate memory + rho = realloc( rho, sizeof(double) * ndofs ); + rho_tracer = realloc( rho_tracer, sizeof(double) * ndofs ); + + // Get density and tracer + trixi_load_conservative_var(handle, 1, rho); + trixi_load_conservative_var(handle, 5, rho_tracer); + + for (int i = 0; i < ndofs; ++i) { + // Apply 5% amplification to tracer (fraction of density) + const double tracer = 1.05 * (rho_tracer[i] / rho[i]); + rho_tracer[i] = tracer * rho[i]; + } + + // Write back tracer + trixi_store_conservative_var(handle, 5, rho_tracer); + } + } + + // 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(); + + free(rho); + free(rho_tracer); + + return 0; +} diff --git a/examples/trixi_controller_data_store.f90 b/examples/trixi_controller_data_store.f90 new file mode 100644 index 00000000..f18b2a7c --- /dev/null +++ b/examples/trixi_controller_data_store.f90 @@ -0,0 +1,115 @@ +program trixi_controller_data_store_f + use LibTrixi + use, intrinsic :: iso_fortran_env, only: error_unit + use, intrinsic :: iso_c_binding, only: c_int, c_double, c_f_pointer, c_ptr + + implicit none + + integer(c_int) :: handle, ndofs, steps, i + character(len=256) :: argument + integer, parameter :: dp = selected_real_kind(12) + real(dp) :: tracer, rho_val, rho_tracer_val + real(dp), dimension(:), pointer :: rho => null(), rho_tracer => null() + type(c_ptr) :: raw_data_c + real(c_double), dimension(:), pointer :: raw_data + + + 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" + call exit(2) + 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" + call exit(2) + end if + + + ! Initialize Trixi + write(*, '(a)') "" + write(*, '(a)') "*** Trixi controller *** Initialize Trixi" + call get_command_argument(1, argument) + call trixi_initialize(argument) + + ! Set up the Trixi simulation + ! We get a handle to use subsequently + write(*, '(a)') "*** Trixi controller *** Set up Trixi simulation" + call get_command_argument(2, argument) + handle = trixi_initialize_simulation(argument) + + ! Main loop + steps = 0 + write(*, '(a)') "*** Trixi controller *** Entering main loop" + + do + ! Exit loop once simulation is completed + if ( trixi_is_finished(handle) ) exit + + call trixi_step(handle) + steps = steps + 1 + + if (modulo(steps, 100) == 0) then + ! get number of degrees of freedom + ndofs = trixi_ndofsglobal(handle) + + ! Get a pointer to Trixi's internal simulation data + raw_data_c = trixi_get_conservative_vars_pointer(handle) + call c_f_pointer(raw_data_c, raw_data, [ndofs]) + + do i = 1,ndofs + ! density comes first + rho_val = raw_data(5*(i-1) + 1 ) + + ! tracer comes last + rho_tracer_val = raw_data(5*i) + + ! Apply 10% damping to tracer (fraction of density) + tracer = 0.9 * (rho_tracer_val / rho_val) + raw_data(5*i) = tracer * rho_val + end do + end if + + if (modulo(steps, 100) == 50) then + ! get number of degrees of freedom + ndofs = trixi_ndofsglobal(handle) + + ! allocate memory + if ( associated(rho) ) deallocate(rho) + allocate( rho(ndofs) ) + if ( associated(rho_tracer) ) deallocate(rho_tracer) + allocate( rho_tracer(ndofs) ) + + ! get density and tracer + call trixi_load_conservative_var(handle, 1, rho) + call trixi_load_conservative_var(handle, 5, rho_tracer) + + do i = 1,ndofs + ! apply 5% amplification to tracer (fraction of density) + tracer = 1.05 * (rho_tracer(i) / rho(i)) + rho_tracer(i) = tracer * rho(i) + end do + + ! write back tracer + call trixi_store_conservative_var(handle, 5, rho_tracer) + end if + end do + + ! Finalize Trixi simulation + write(*, '(a)') "" + write(*, '(a)') "*** Trixi controller *** Finalize Trixi simulation" + call trixi_finalize_simulation(handle) + + ! Finalize Trixi + write(*, '(a)') "" + write(*, '(a)') "*** Trixi controller *** Finalize Trixi" + call trixi_finalize() + + deallocate(rho) + nullify(rho) + deallocate(rho_tracer) + nullify(rho_tracer) +end program diff --git a/src/api.c b/src/api.c index da2b17d2..74f0d5f8 100644 --- a/src/api.c +++ b/src/api.c @@ -22,9 +22,12 @@ enum { TRIXI_FPTR_NNODES, TRIXI_FPTR_LOAD_NODE_REFERENCE_COORDINATES, TRIXI_FPTR_LOAD_NODE_WEIGHTS, - TRIXI_FTPR_LOAD_PRIMITIVE_VARS, - TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VARS, + TRIXI_FPTR_LOAD_CONSERVATIVE_VAR, + TRIXI_FTPR_LOAD_PRIMITIVE_VAR, + TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VAR, + TRIXI_FPTR_STORE_CONSERVATIVE_VAR, TRIXI_FTPR_REGISTER_DATA, + TRIXI_FPTR_GET_CONSERVATIVE_VARS_POINTER, TRIXI_FTPR_VERSION_LIBRARY, TRIXI_FTPR_VERSION_LIBRARY_MAJOR, TRIXI_FTPR_VERSION_LIBRARY_MINOR, @@ -60,9 +63,12 @@ static const char* trixi_function_pointer_names[] = { [TRIXI_FPTR_NNODES] = "trixi_nnodes_cfptr", [TRIXI_FPTR_LOAD_NODE_REFERENCE_COORDINATES] = "trixi_load_node_reference_coordinates_cfptr", [TRIXI_FPTR_LOAD_NODE_WEIGHTS] = "trixi_load_node_weights_cfptr", - [TRIXI_FTPR_LOAD_PRIMITIVE_VARS] = "trixi_load_primitive_vars_cfptr", - [TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VARS] = "trixi_load_element_averaged_primitive_vars_cfptr", + [TRIXI_FPTR_LOAD_CONSERVATIVE_VAR] = "trixi_load_conservative_var_cfptr", + [TRIXI_FTPR_LOAD_PRIMITIVE_VAR] = "trixi_load_primitive_var_cfptr", + [TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VAR] = "trixi_load_element_averaged_primitive_var_cfptr", + [TRIXI_FPTR_STORE_CONSERVATIVE_VAR] = "trixi_store_conservative_var_cfptr", [TRIXI_FTPR_REGISTER_DATA] = "trixi_register_data_cfptr", + [TRIXI_FPTR_GET_CONSERVATIVE_VARS_POINTER] = "trixi_get_conservative_vars_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", @@ -644,7 +650,33 @@ void trixi_load_node_weights(int handle, double* node_weights) { /** - * @anchor trixi_load_primitive_vars_api_c + * @anchor trixi_load_conservative_var_api_c + * + * @brief Load conservative variable. + * + * The values for the conservative variable at position `variable_id` at every degree of + * freedom are stored in the given array `data`. + * + * The given array has to be of correct size (ndofs) and memory has to be allocated + * beforehand. + * + * @param[in] handle simulation handle + * @param[in] variable_id index of variable + * @param[out] data values for all degrees of freedom + */ +void trixi_load_conservative_var(int handle, int variable_id, double * data) { + + // Get function pointer + void (*load_conservative_var)(int, int, double *) = + trixi_function_pointers[TRIXI_FPTR_LOAD_CONSERVATIVE_VAR]; + + // Call function + return load_conservative_var(handle, variable_id, data); +} + + +/** + * @anchor trixi_load_primitive_var_api_c * * @brief Load primitive variable * @@ -658,19 +690,19 @@ void trixi_load_node_weights(int handle, double* node_weights) { * @param[in] variable_id index of variable * @param[out] data values for all degrees of freedom */ -void trixi_load_primitive_vars(int handle, int variable_id, double * data) { +void trixi_load_primitive_var(int handle, int variable_id, double * data) { // Get function pointer - void (*load_primitive_vars)(int, int, double *) = - trixi_function_pointers[TRIXI_FTPR_LOAD_PRIMITIVE_VARS]; + void (*load_primitive_var)(int, int, double *) = + trixi_function_pointers[TRIXI_FTPR_LOAD_PRIMITIVE_VAR]; // Call function - load_primitive_vars(handle, variable_id, data); + load_primitive_var(handle, variable_id, data); } /** - * @anchor trixi_load_element_averaged_primitive_vars_api_c + * @anchor trixi_load_element_averaged_primitive_var_api_c * * @brief Load element averages for primitive variable * @@ -684,14 +716,39 @@ void trixi_load_primitive_vars(int handle, int variable_id, double * data) { * @param[in] variable_id index of variable * @param[out] data element averaged values for all elements */ -void trixi_load_element_averaged_primitive_vars(int handle, int variable_id, double * data) { +void trixi_load_element_averaged_primitive_var(int handle, int variable_id, double * data) { // Get function pointer - void (*load_element_averaged_primitive_vars)(int, int, double *) = - trixi_function_pointers[TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VARS]; + void (*load_element_averaged_primitive_var)(int, int, double *) = + trixi_function_pointers[TRIXI_FTPR_LOAD_ELEMENT_AVERAGED_PRIMITIVE_VAR]; // Call function - load_element_averaged_primitive_vars(handle, variable_id, data); + load_element_averaged_primitive_var(handle, variable_id, data); +} + + +/** + * @anchor trixi_store_conservative_var_api_c + * + * @brief Store conservative variable. + * + * The values for the conservative variable at position `variable_id` at every degree of + * freedom are read from the given array `data` and written to Trixi.jl's internal storage. * + * + * The given array has to be of correct size (ndofs). + * + * @param[in] handle simulation handle + * @param[in] variable_id index of variable + * @param[in] data values for all degrees of freedom + */ +void trixi_store_conservative_var(int handle, int variable_id, double * data) { + + // Get function pointer + void (*store_conservative_var)(int, int, double *) = + trixi_function_pointers[TRIXI_FPTR_STORE_CONSERVATIVE_VAR]; + + // Call function + return store_conservative_var(handle, variable_id, data); } @@ -726,6 +783,30 @@ void trixi_register_data(int handle, int index, int size, const double * data) { } +/** + * @anchor trixi_get_conservative_vars_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_conservative_vars_pointer(int handle) { + + // Get function pointer + double * (*get_conservative_vars_pointer)(int) = + trixi_function_pointers[TRIXI_FPTR_GET_CONSERVATIVE_VARS_POINTER]; + + // Call function + return get_conservative_vars_pointer(handle); +} + + /** * @anchor trixi_get_simulation_time_api_c * diff --git a/src/api.f90 b/src/api.f90 index 223ece95..1b7b2933 100644 --- a/src/api.f90 +++ b/src/api.f90 @@ -375,7 +375,28 @@ subroutine trixi_load_node_weights(handle, node_weights) bind(c) end subroutine !> - !! @fn LibTrixi::trixi_load_primitive_vars::trixi_load_primitive_vars(handle, variable_id, data) + !! @anchor trixi_load_conservative_var_api_c + !! + !! @brief Load conservative variable. + !! + !! The values for the conservative variable at position `variable_id` at every degree of + !! freedom are stored in the given array `data`. + !! + !! The given array has to be of correct size (ndofs) and memory has to be allocated + !! beforehand. + !! + !! @param[in] handle simulation handle + !! @param[in] variable_id index of variable + !! @param[out] data values for all degrees of freedom + subroutine trixi_load_conservative_var(handle, variable_id, data) bind(c) + use, intrinsic :: iso_c_binding, only: c_int, c_double + integer(c_int), value, intent(in) :: handle + integer(c_int), value, intent(in) :: variable_id + real(c_double), dimension(*), intent(out) :: data + end subroutine + + !> + !! @fn LibTrixi::trixi_load_primitive_var::trixi_load_primitive_var(handle, variable_id, data) !! !! @brief Load primitive variable !! @@ -383,8 +404,8 @@ subroutine trixi_load_node_weights(handle, node_weights) bind(c) !! @param[in] variable_id index of variable !! @param[out] data primitive variable values for all degrees of freedom !! - !! @see @ref trixi_load_primitive_vars_api_c "trixi_load_primitive_vars (C API)" - subroutine trixi_load_primitive_vars(handle, variable_id, data) bind(c) + !! @see @ref trixi_load_primitive_var_api_c "trixi_load_primitive_var (C API)" + subroutine trixi_load_primitive_var(handle, variable_id, data) bind(c) use, intrinsic :: iso_c_binding, only: c_int, c_double integer(c_int), value, intent(in) :: handle integer(c_int), value, intent(in) :: variable_id @@ -407,7 +428,7 @@ real(c_double) function trixi_get_simulation_time(handle) bind(c) end function !> - !! @fn LibTrixi::trixi_load_element_averaged_primitive_vars::trixi_load_element_averaged_primitive_vars(handle, variable_id, data) + !! @fn LibTrixi::trixi_load_element_averaged_primitive_var::trixi_load_element_averaged_primitive_var(handle, variable_id, data) !! !! @brief Load element averages for primitive variable !! @@ -415,14 +436,35 @@ real(c_double) function trixi_get_simulation_time(handle) bind(c) !! @param[in] variable_id index of variable !! @param[out] data averaged values for all elements !! - !! @see @ref trixi_load_element_averaged_primitive_vars_api_c "trixi_load_element_averaged_primitive_vars (C API)" - subroutine trixi_load_element_averaged_primitive_vars(handle, variable_id, data) bind(c) + !! @see @ref trixi_load_element_averaged_primitive_var_api_c "trixi_load_element_averaged_primitive_var (C API)" + subroutine trixi_load_element_averaged_primitive_var(handle, variable_id, data) bind(c) use, intrinsic :: iso_c_binding, only: c_int, c_double integer(c_int), value, intent(in) :: handle integer(c_int), value, intent(in) :: variable_id real(c_double), dimension(*), intent(out) :: data end subroutine + !> + !! @anchor trixi_store_conservative_var_api_c + !! + !! @brief Store conservative variable. + !! + !! The values for the conservative variable at position `variable_id` at every degree of + !! freedom are read from the given array `data` and written to Trixi.jl's internal + !! storage. + !! + !! The given array has to be of correct size (ndofs). + !! + !! @param[in] handle simulation handle + !! @param[in] variable_id index of variable + !! @param[in] data values for all degrees of freedom + subroutine trixi_store_conservative_var(handle, variable_id, data) bind(c) + use, intrinsic :: iso_c_binding, only: c_int, c_double + integer(c_int), value, intent(in) :: handle + integer(c_int), value, intent(in) :: variable_id + real(c_double), dimension(*), intent(in) :: data + end subroutine + !> !! @fn LibTrixi::trixi_register_data::trixi_register_data(handle, variable_id, data) !! @@ -442,6 +484,24 @@ subroutine trixi_register_data(handle, index, size, data) bind(c) real(c_double), dimension(*), intent(in) :: data end subroutine + !> + !! @anchor trixi_get_conservative_vars_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_conservative_vars_pointer(handle) bind(c) + use, intrinsic :: iso_c_binding, only: c_int, c_ptr + integer(c_int), value, intent(in) :: handle + end function + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/trixi.h b/src/trixi.h index b8d3f433..7abeacc1 100644 --- a/src/trixi.h +++ b/src/trixi.h @@ -37,9 +37,12 @@ double trixi_calculate_dt(int handle); double trixi_get_simulation_time(int handle); void trixi_load_node_reference_coordinates(int handle, double* node_coords); void trixi_load_node_weights(int handle, double* node_weights); -void trixi_load_primitive_vars(int handle, int variable_id, double * data); -void trixi_load_element_averaged_primitive_vars(int handle, int variable_id, double * data); +void trixi_load_conservative_var(int handle, int variable_id, double * data); +void trixi_load_primitive_var(int handle, int variable_id, double * data); +void trixi_load_element_averaged_primitive_var(int handle, int variable_id, double * data); +void trixi_store_conservative_var(int handle, int variable_id, double * data); void trixi_register_data(int handle, int index, int size, const double * data); +double * trixi_get_conservative_vars_pointer(int handle); // T8code #if !defined(T8_H) && !defined(T8_FOREST_GENERAL_H) diff --git a/test/c/simulation.cpp b/test/c/simulation.cpp index 36fbccbb..56320f85 100644 --- a/test/c/simulation.cpp +++ b/test/c/simulation.cpp @@ -100,11 +100,21 @@ TEST(CInterfaceTest, SimulationRun) { } EXPECT_NEAR(integral, 0.4, 1e-17); - // Check primitive variable values on all dofs + // Check conservative variable values std::vector rho(ndofs); + std::vector rho_energy(ndofs); + trixi_load_conservative_var(handle, 1, rho.data()); + trixi_load_conservative_var(handle, 4, rho_energy.data()); + // check memory borders + EXPECT_DOUBLE_EQ(rho[0], 1.0); + EXPECT_DOUBLE_EQ(rho[ndofs-1], 1.0); + EXPECT_DOUBLE_EQ(rho_energy[0], 2.5e-5); + EXPECT_DOUBLE_EQ(rho_energy[ndofs-1], 2.5e-5); + + // Check primitive variable values std::vector energy(ndofs); - trixi_load_primitive_vars(handle, 1, rho.data()); - trixi_load_primitive_vars(handle, 4, energy.data()); + trixi_load_primitive_var(handle, 1, rho.data()); + trixi_load_primitive_var(handle, 4, energy.data()); // check memory borders EXPECT_DOUBLE_EQ(rho[0], 1.0); EXPECT_DOUBLE_EQ(rho[ndofs-1], 1.0); @@ -116,10 +126,10 @@ TEST(CInterfaceTest, SimulationRun) { std::vector v1_averages(nelements); std::vector v2_averages(nelements); std::vector e_averages(nelements); - trixi_load_element_averaged_primitive_vars(handle, 1, rho_averages.data()); - trixi_load_element_averaged_primitive_vars(handle, 2, v1_averages.data()); - trixi_load_element_averaged_primitive_vars(handle, 3, v2_averages.data()); - trixi_load_element_averaged_primitive_vars(handle, 4, e_averages.data()); + trixi_load_element_averaged_primitive_var(handle, 1, rho_averages.data()); + trixi_load_element_averaged_primitive_var(handle, 2, v1_averages.data()); + trixi_load_element_averaged_primitive_var(handle, 3, v2_averages.data()); + trixi_load_element_averaged_primitive_var(handle, 4, e_averages.data()); if (nranks == 1) { // check memory borders (densities at the beginning, energies at the end) EXPECT_DOUBLE_EQ(rho_averages[0], 1.0); @@ -166,6 +176,15 @@ TEST(CInterfaceTest, SimulationRun) { FAIL() << "Test cannot be run with " << nranks << " ranks."; } + // Check storing of conservative variables + rho[0] = 42.0; + rho[ndofs-1] = 23.0; + trixi_store_conservative_var(handle, 1, rho.data()); + + double * raw_data = trixi_get_conservative_vars_pointer(handle); + EXPECT_DOUBLE_EQ(rho[0], raw_data[0]); + EXPECT_DOUBLE_EQ(rho[ndofs-1], raw_data[4*(ndofs-1)]); + // Finalize Trixi simulation trixi_finalize_simulation(handle); diff --git a/test/c/t8code.cpp b/test/c/t8code.cpp index 9339b9b1..18bf50e4 100644 --- a/test/c/t8code.cpp +++ b/test/c/t8code.cpp @@ -9,7 +9,7 @@ const char * julia_project_path = JULIA_PROJECT_PATH; // Example libexlixir const char * libelixir_path = - "../../../LibTrixi.jl/examples/libelixir_t8code2d_advection_amr.jl"; + "../../../LibTrixi.jl/examples/libelixir_t8code2d_euler_tracer_amr.jl"; TEST(CInterfaceTest, T8code) { diff --git a/test/fortran/simulationRun_suite.f90 b/test/fortran/simulationRun_suite.f90 index c9261f02..ed40707c 100644 --- a/test/fortran/simulationRun_suite.f90 +++ b/test/fortran/simulationRun_suite.f90 @@ -1,6 +1,7 @@ module simulationRun_suite use LibTrixi use testdrive, only : new_unittest, unittest_type, error_type, check + use, intrinsic :: iso_c_binding, only: c_double, c_f_pointer, c_ptr implicit none private @@ -29,6 +30,8 @@ subroutine test_simulationRun(error) integer, parameter :: dp = selected_real_kind(15) real(dp) :: dt, time, integral real(dp), dimension(:), allocatable :: data, weights + type(c_ptr) :: raw_data_c + real(c_double), dimension(:), pointer :: raw_data ! Initialize Trixi call trixi_initialize(julia_project_path) @@ -99,10 +102,16 @@ subroutine test_simulationRun(error) call check(error, integral, 0.4_dp) deallocate(data) - ! Check primitive variable values + ! Check conservative variable values size = ndofs allocate(data(size)) - call trixi_load_primitive_vars(handle, 1, data) + call trixi_load_conservative_var(handle, 1, data) + call check(error, data(1), 1.0_dp) + call check(error, data(3200), 1.0_dp) + call check(error, data(size), 1.0_dp) + + ! Check primitive variable values + call trixi_load_primitive_var(handle, 1, data) call check(error, data(1), 1.0_dp) call check(error, data(3200), 1.0_dp) call check(error, data(size), 1.0_dp) @@ -111,10 +120,21 @@ subroutine test_simulationRun(error) ! Check element averaged values size = nelements allocate(data(size)) - call trixi_load_element_averaged_primitive_vars(handle, 1, data) + call trixi_load_element_averaged_primitive_var(handle, 1, data) call check(error, data(1), 1.0_dp) call check(error, data(94), 0.99833232379996562_dp) call check(error, data(size), 1.0_dp) + + ! Check storing of conservative variables + data(1) = 42.0 + data(ndofs) = 23.0 + call trixi_store_conservative_var(handle, 1, data) + + raw_data_c = trixi_get_conservative_vars_pointer(handle) + call c_f_pointer(raw_data_c, raw_data, [ndofs]) + call check(error, data(1), raw_data(1)) + call check(error, data(ndofs), raw_data(4*ndofs - 3)) + deallocate(data) ! Finalize Trixi simulation diff --git a/test/fortran/t8code_suite.f90 b/test/fortran/t8code_suite.f90 index 24d71545..e95afd84 100644 --- a/test/fortran/t8code_suite.f90 +++ b/test/fortran/t8code_suite.f90 @@ -10,7 +10,7 @@ module t8code_suite character(len=*), parameter, public :: julia_project_path = JULIA_PROJECT_PATH character(len=*), parameter, public :: libelixir_path = & - "../../../LibTrixi.jl/examples/libelixir_t8code2d_advection_amr.jl" + "../../../LibTrixi.jl/examples/libelixir_t8code2d_euler_tracer_amr.jl" contains