[ANN] Trixi.jl v0.3: SciML integration and a new modular approach for easy extension

Sounds great! As @ranocha stated before, let us know if you have questions about using Trixi, either here or by creating an issue in GitHub.

Our current status is that we have a working implementation of your EC mortar flux correction in #247, but it needs to be ported to the newly released Trixi (it is still based on our earlier, monolithic approach). Once that works again, the next step would be to extend it to ES mortars, which shouldn’t be too hard (I think).

1 Like

Our current status is that we have a working implementation of your EC mortar flux correction in #247, but it needs to be ported to the newly released Trixi (it is still based on our earlier, monolithic approach). Once that works again, the next step would be to extend it to ES mortars, which shouldn’t be too hard (I think).

That was fast! And agreed, I don’t think ES mortars (assuming it’s from some entropy dissipative flux penalization) would be very difficult.

We now natively support visualizing results with Plots.jl both during and after a simulation. With the current Trixi.jl release (v0.3.6), plotting a solution from the REPL is as easy as loading Plots and throwing the ODESolution object from the call to DiffEqBase.solve! at it:

julia> trixi_include(default_example())
[...]

julia> using Plots

julia> plot(sol)

For in-situ visualization, use a VisualizationCallback to create a new plot of the current solution in regular time step intervals:
trixi-in-situ-visualization

The features and how to use them are described in detail in the official docs. We are looking for feedback on how to further improve the visualization/postprocessing capabilities in Trixi, so please feel welcome to leave a message here or create an issue.

Thanks @rvelts once more for these valuable suggestions!

11 Likes

That’s amazing.

5 Likes

That’s really cool!! I think it will boost the use of trixi (at least mine!).

On a side note, I cannot run your examples, like /examples/1d/ elixir_advection_amr.jl
For one, it says LinearScalarAdvectionEquation1D not defined

1 Like

Hm, that should not happen. If you create an issue on GitHub (with a little more info on how to reproduce your error), I will look into it. In general, all elixirs are run at least once during CI testing, where all tests are passing, thus this should not be too difficult to figure out.

Vielen Dank for making Trixi.jl available to my students!

My students would like to define the velocity vector by separately solving the eikonal equation and interpolating the velocity vector to the mesh that Trixi.jl employs.

How do you suggest to proceed?

Vielen dank und grusse.

1 Like

Dear @ziolai,
thank you very much! That sounds like an interesting application and we would like to know more about it - since we also need a bit more information about the problem to give a good answer. Would you mind opening an issue in our repository?
Some first hints: I guess you want to solve some linear advection equation. If so, you can setup a new equation type as described in our tutorials, e.g., this one on a conservative equation or this one on a non-conservative equation with variable coefficients. Then, you could solve the eikonal equation in whatever way you like and interpolate it to our solution nodes. You could for example wrap the interpolation object as callable function and pass that when creating the semidiscretization of Trixi.jl. Alternatively, you oculd initialize the variable coefficients with some default values and overwrite them after setting up the data, e.g., after setting up the ODE problem (via semidiscretize). If you do that, you need o get the node coordinates, e.g., via

semi = SemidiscretizationHyperbolic(...)

For example, you can get the following on a TreeMesh in 2D:

julia> using Trixi

julia> trixi_include(default_example())

julia> semi.cache.elements.node_coordinates # undocumented internal state - be cautious
2×4×4×256 Array{Float64, 4}:
[...]

Here, the first dimension of the array is the spatial dimension (2D in this case), the second and third dimensions are the tensor-product nodes of our Gauss-Lobatto-Legendre nodes in a quad element, and the last dimension of the array is the element index. If you want to have a plain vector of structs, you can wrap this, e.g., as

julia> reinterpret(SVector{2, Float64}, vec(semi.cache.elements.node_coordinates))
4096-element reinterpret(SVector{2, Float64}, ::Vector{Float64}):
 [-1.0, -1.0]
 [-0.9654508497187474, -1.0]
 [-0.9095491502812526, -1.0]
 [-0.875, -1.0]
[...]

Dear Hendrik,

Thank you so much for providing a wealth of information!

A group of four master students is looking into the mondeling of macroscopic pedestrian flow models based on a 2013 paper of Twarogowska, Goatin en Duvigneau. This assignment is part of compulsory mathematical modelling course that is part of their curriculum.

Please allow me to discuss the two options that you suggest with the students and get back to you next week.

Happy Easter weekend, Domenico.

2 Likes

Hello @ranocha ! I am interested in the package. Do you know if it can be used to solve the systems of partial differential equations that appear in General Relativity ?
Kind regards!

1 Like

Thanks for your interest! I don’t have experience with numerical solutions of the equations of general relativity. Could you please open an issue and provide a few more details there, e.g., equations desrcribing a minimal example that could be interesting for you - say with reduced dimensions etc. and not full 3+1D?

1 Like

Please advice on how to replace flux(u) = u^3 by e.g. flux(u,x) = x * u (or flux(u) = f(x) g(u)) in Tutorial 5: “Adding New Conservation Law”.

In a next step, we will be concerned with replacing f(x) with interpolation.

In a next-next step, we will concerned with loading the data to interpolate from file and the extension from 1D to 2D.

Thx! D.

https://trixi-framework.github.io/Trixi.jl/stable/tutorials/adding_nonconservative_equation/

Danke!

So to make sure that my students and I understand this well: there is no way to adapt the conservative formulation, i.e., Tutorial 5?

We desire to solve \partial u / \partial t + \partial [ f(x) g(u) ] / \partial x = 0, where f(x) is given by interpolation. Are you suggestion to interpolate using cubic splines such that the interpolation that defines f((x) can be differentiated?

Alternatively, would it be feasible to solve the eikonal equation using Trixi, first separately and later coupled with the transport equation. This would leverage the interpolation to functionality already in place in Trixi.

Edit: even with the differentiation of f(x), the flux remains explicitly dependent upon x. Are you suggesting to employ Tutorial 6 to make the conservative part of the flux x-dependent?

Thx!

The current version of Trixi.jl does not handle variable coefficients natively (see https://github.com/trixi-framework/Trixi.jl/issues/358). Thus, you need to introduce an auxiliary conserved variable as described in the tutorial on a non-conservative equation I linked above. You can of course use this setup for a conservative system. That’s what we do for our acoustics equations.

Trixi.jl will initialize the conserved variables (including the auxiliary variable for f(x)) at the discretization nodes. It doesn’t need to be smooth to run at all, but if your input is smooth, it can only help to preserve that (e.g., using splines as you wrote).

Right!

Suppose that we modify Tutorial 6 to solve the equation u_t + [ a(x,t) u ]_x = 0 instead of the current form u_t + a(x,t) u_x = 0 (i.e. placing the function a(x,t) inside of the x-derivative). How do you modify

function flux_nonconservative(u_mine, u_other, orientation,
                              equations::NonconservativeLinearAdvectionEquation)
    _, advection_velocity = u_mine
    scalar, _             = u_other

    return SVector(advection_velocity * scalar, zero(scalar))
end

accordingly?

Does

function flux_nonconservative(u_mine, u_other, orientation,
                              equations::NonconservativeLinearAdvectionEquation)
    scalar, advection_velocity  = u_other

    return SVector(advection_velocity * scalar, zero(scalar))
end

suffice? Thx.

If you don’t have a nonconservative equation, you do not need any flux_nonconservative - only Trixi.flux as in 7 Adding a new scalar conservation law · Trixi.jl

Or, alternatively, (staying with Tutorial 6), and do

flux(u, orientation, equation::NonconservativeLinearAdvectionEquation) = SVector(u[1]*u[2], zero(scalar) )

have_nonconservative_terms(::NonconservativeLinearAdvectionEquation) = Val(false)

?

have_nonconservative_terms(::NonconservativeLinearAdvectionEquation) = Val(false) is the default.
flux(u, orientation, equation::NonconservativeLinearAdvectionEquation) = SVector(u[1]*u[2], zero(scalar) ) must be something like flux(u, orientation, equation::NonconservativeLinearAdvectionEquation) = SVector(u[1]*u[2], zero(eltype(u))) (since scalar is not defined…).

Ganz sch"on! So far, so good!

Next, we would like to replace
advection_velocity = 2 + cos(x[1])
by
advection_velocity = interpolate(x[1],data-from-eikonal-solve)

Is is feasible to define data-from-eikonal-solve as a data member of the struct NonconservativeLinearAdvectionEquation as examples show in defining gamma as a data member of the CompressibleEulerEquations2D?

Edit: I have the interpolation in place. All good.