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

Package Announcements Package announcements
, , , ,
#41

Yes, I don’t see a problem there.

#42

Thanks for all the help! (I’m one of the students of @ziolai). With your and his help, we managed to implement the linear advection equation for a 2D-problem (with a space-dependent flux function) for a space-dependent fixed advection velocity field (given in the initial condition):

u_t +\nabla \cdot (\mathbf{a}(\mathbf{x}) \; u(\mathbf{x}, t)) = 0

However, we want to make the velocity field dependent of u. This means that we need to solve the eikonal equation, given the density field, for each time step in the time integration again, like

u_t +\nabla \cdot (\mathbf{a}(\mathbf{x}, u) \; u(\mathbf{x}, t)) = 0

where a(\mathbf{x}, u) is the advection velocity field computed by solving the eikonal equation for each time step using the values of u (\mathbf{x}, t). Do you have any suggestions for how we can make the velocity field dependent of u like in this equation?

#43

In short:

ode_baseline = semidiscretize(semi, tspan) # as usual

function modified_rhs!(du_ode, u_ode, semi, t)
  # solve eikonal equation given u_ode (vector of numerical solution values)
  # update coefficients
  Trixi.rhs!(du_ode, u_ode, semi, t)
end
ode = ODEProblem(modified_rhs!, ode_baseline.u0_ode, ode_baseline.tspan, ode_baseline.p #= = semi =#)

# usual solve setup etc.
2 Likes
#44

Thank you very much for your reaction! Two questions about your answer:

  • Is it correct that the variable u_ode is of the same form as ode_baseline.u0?
  • If we implement it like this, it seems we can’t use the Trixi.PlotData2D functions anymore, because the modified_rhs! is not of the type Trixi.rhs!. Is there any solution to this problem? Can we adjust the function or the Trixi.PlotData2D functionality such that it also accepts the modified function?
#45

Yes :+1:

You should be able to call pd = PlotData2D(sol.u[end], sol.prob.p) as in Trixi.jl/types.jl at d1f98e854f3a84e59c4c9e38a12a95349b2b8be5 · trixi-framework/Trixi.jl · GitHub. But we should probably just drop the restriction to Trixi.rhs! in Trixi.jl/types.jl at d1f98e854f3a84e59c4c9e38a12a95349b2b8be5 · trixi-framework/Trixi.jl · GitHub. I made a PR to fix this. It would be great if you could let me know whether that works in your case.

1 Like
#46

We just made a new release of Trixi.jl including the PR I linked above. Could you please check whether that fixes the plotting problem you encountered?

#47

Thank you very much for your help! The new release indeed seems to fix our plotting problem.

2 Likes
#48

image

I performed the implementation as you suggested here and it seems to work for the biggest part. For a given velocity function \mathbf{a}(\mathbf{x}, u) we get results as shown in the plot above (on the left we see the density field u, on the right the x-component of the velocity field a_x.
What I’m wondering however is why the color map of this velocity field is so “discrete”: only at the positions of the grid centers the value changes to the value of the updated velocity field. Does this mean that my implementation is wrong or is there some problem with the visualization of the updated velocity field?

Edit: We used the following velocity field function \mathbf{a}(\mathbf{x}, u) to test the implementation:

function advection_vel(x, u)
    a_x = max(2*u^2 + x[1]^2, 0.1) # define the x-component of the velocity
    a_y = max(2*u^2 + x[2]^2, 0.1) # define the y-component of the velocity
    return a_x, a_y
end
#49

It’s hard to help since I don’t know some details. Some questions/remarks:

  • How do you set the advection velocity in the end? Did you verify that you get the correct values at the grid nodes?
  • If you update the advection velocity based on the solution itself, you end up with a nonlinear problem - and can get all the problems this means (possible discontinuities etc.).
  • We probably interpolate to a uniform grid for plotting - did you try changing the grid resolution of the numerical solution and/or the plotting resolution (via keyword arguments of PlotData2D, see the docstring at Trixi.jl/types.jl at d9d469fc9c32aba129e7ce7976f340b2afa6b4aa · trixi-framework/Trixi.jl · GitHub)?
1 Like
#50

I checked the values again and it seems that they are incorrect - the value of a_x should be at least 0.1 everywhere on the domain (because of the definition of \mathbf{a}(\mathbf{x}, u), but it is zero on most of the grid nodes for some reason. So, something seems to be wrong with the implementation of the modified_rhs!. This function is defined as follows now:

function modified_rhs!(du_ode, u_ode, semi, t)
    u_sol = reinterpret(SVector{3, Float64}, vec(u_ode)) # reinterpret the solution values (u, a_x, a_y)
    us = [elm[1] for elm in u_sol]
    x = reinterpret(SVector{2, Float64}, vec(semi.cache.elements.node_coordinates))
    a = [advection_vel(x[i], us[i]) for i=1:length(x)] # advection_vel is the function defined in my previous message
    ax = [elm[1] for elm in a]
    ay = [elm[2] for elm in a]
    u_sol = [SVector(us[i], ax[i], ay[i]) for i in 1:length(us)]
    u_sol = reinterpret(SVector{1, Float64}, u_sol)
    u_odenew = [u[1] for u in u_sol] # this is necessary to get u_odenew in the right format (Vector{Float64}}
    Trixi.rhs!(du_ode, u_odenew, semi, t)
end

Is maybe something wrong in using reinterpret(SVector{2, Float64}, vec(semi.cache.elements.node_coordinates)) to get the coordinates of the grid nodes?

#51

I guess you create a modified_ode and solve it to get a modified_sol? Note that you do not update your u_ode in-place. Thus, the new advection velocity is only visible in Trixi.rhs! but what you get in modified_sol.u[end] (the final solution) as advection velocity is basically garbage. You should plot your u_odenew to see whether things work as they should.

1 Like
#52

Ah wow, thank you very much for pointing that out! I was indeed comparing the wrong values, forgot to update the advection velocity in the modified_sol. I’ve implemented a function very similar to the modifed_rhs! now which updates the advection velocity after the solution is returned, and these results seem to be correct (however we will probably still need to look into the problems that the nonlinear form brings with itself).

The plotting problem is by the way also fixed after updating the modified_sol, the plots of the updated solution looks much nicer, as you can see below (again, the density u is on the left hand side and the x-component of the velocity field on the left hand side).

image

2 Likes