trixi-framework / Trixi.jl

Trixi.jl: Adaptive high-order numerical simulations of conservation laws in Julia
https://trixi-framework.github.io/Trixi.jl
MIT License
533 stars 108 forks source link

RFC: alternative computation of source terms based on indices #2010

Open benegee opened 3 months ago

benegee commented 3 months ago

In our current project ADAPTEX we use Trixi.jl/TrixiAtmo.jl via libtrixi. In a typical workflow, source terms (aka forcings) will be computed by the external application and stored as a vector, which will be accessible to Trixi.jl. In such a situation the typical callback

du_local = source_terms(u_local, x_local, t, equations)

is difficult to use. We would rather use the element and node indices to figure out which values to read from the supplied vector. https://github.com/trixi-framework/Trixi.jl/blob/d0d704dc5948c15cb3367bb1c2d3bb63235d4d28/src/solvers/dgsem_tree/dg_3d.jl#L1390-L1398

Would it be possible to implement an alternative callback, which could then be somehow selected. Maybe based on a trait similar to have_nonconservative_terms but for source_terms?

ranocha commented 3 months ago

I do not see a clean and nice way how to do this (but this may be my fault). Instead, I would create another ODE RHS in libtrixi that adds the required source terms via indices. It kind of has to live there since the conversion of the indices to the external data is something Trixi.jl does not know.

lchristm commented 3 months ago

I think the simplest way to do this would be to overload the calc_sources! function in the elixir that is used with libtrixi. E.g. define something like

struct Forcing
    data::Vector{Float64} # received from external application
    dof2index::Array{Int, 4} # maps (i,j,k,elem) to idx in data
end
function Trixi.calc_sources!(du, u, t, source_terms::Forcing,
                             equations::AbstractEquations{3}, dg::DG, cache)
    # do something
end

Note that overloading an internal function leads to stricter compat limits, i.e. such elixir is not guaranteed to work with any future version of Trixi.jl. But for initial experimentation this may be fine.