Yes, I don’t see a problem there.

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?

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.
```

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?

Yes

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.

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?

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

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
```

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)?

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?

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.

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).