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.
2 Likes
Thank you very much for your reaction! Two questions about your answer:
- Is it correct that the variable
u_odeis of the same form asode_baseline.u0? - If we implement it like this, it seems we can’t use the
Trixi.PlotData2Dfunctions anymore, because themodified_rhs!is not of the typeTrixi.rhs!. Is there any solution to this problem? Can we adjust the function or theTrixi.PlotData2Dfunctionality 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 https://github.com/trixi-framework/Trixi.jl/blob/d1f98e854f3a84e59c4c9e38a12a95349b2b8be5/src/visualization/types.jl#L286. But we should probably just drop the restriction to Trixi.rhs! in https://github.com/trixi-framework/Trixi.jl/blob/d1f98e854f3a84e59c4c9e38a12a95349b2b8be5/src/visualization/types.jl#L7. 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
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.
2 Likes
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)?
1 Like
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.
1 Like
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).
2 Likes
Hey there! I’m trying to follow this tutorial for HOHQmesh (14 Unstructured meshes with HOHQMesh.jl · Trixi.jl), specifically the ice cream mesh example of section 9: Unstructured meshes with HOHQMesh.jl. Right now I’m simply trying to run their example. I’ve copied the control file into a seperate julia jupyter notebook, and have the mesh generating code they provide:
using HOHQMesh
control_file = joinpath(“out”, “ice_cream_straight_sides.control”)
output = generate_mesh(control_file)
which outputs the following when run:
" \r\n *******************\r\n 2D Mesh Statistics:\r\n *******************\r\n Total time = 1.5625000000000000E-002\r\n Number of nodes = 296\r\n Number of Edges = 548\r\n Number of Elements = 252\r\n \r\n Mesh Quality:\r\n Me" ⋯ 747 bytes ⋯ "000000 90.00000000 90.00000000\r\n Maximum Angle 90.00200200 135.35199230 94.62391514 90.00000000 135.00000000 90.00000000\r\n Area Sign 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000\r\n "
This creates two files I believe, the .mesh and .tec files. In the tutorial it says the .tec file is used to plot the mesh grid, but it doesn’t particularly indicate how to plot it. I attempted to plot this ice_cream_straight_sides.tec file using Plots with the following:
using Plots
plot(ice_cream_straight_sides.tec)
But it outputs:
UndefVarError: ice_cream_straight_sides not defined
Stacktrace:
[1] top-level scope
@ In[15]:5
[2] eval
@ .\boot.jl:373 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1196
I’m very new to Julia so I’m not sure what went wrong or how to proceed, I’d greatly appreciate any help to solve this!
additional info: I’m running this on jupyter notebooks with windows 11.
Hey! The output for the generate mesh will display in a better format if you add a println(output). For example,
using HOHQMesh
control_file = joinpath(“out”, “ice_cream_straight_sides.control”)
output = generate_mesh(control_file); println(output)
Currently, the .tec file must be visualized in an external program like ParaView (https://www.paraview.org/). In this program you can open the ice_cream_straight_sides.tec and choose the “TecPlot read in” to draw the mesh grid. We are working on an interactive tool (Add interactive HQMTool and testing by DavidAKopriva · Pull Request #15 · trixi-framework/HOHQMesh.jl · GitHub) that can generate and draw the mesh project directly in a REPL session. For this interactive tool visualization we are using Makie.
Besides the visualization tooling that Andrew mentioned, your plot command isn’t valid as written. You’re seeing a missing variable error because Julia is looking for a nonexistent variable called ice_cream_straight_sides containing a field tec, but filepaths need to be enclosed in double-quotes like ”path/to/file.ext”.
Thanks so much for the help @awinters, I’ll try paraview!
Another question I have now: Is it possible to generate a regular mesh with no interior boundaries (so no ice cream cone inside, just the background grid) by adapting the ice cream control file? Perhaps by deleting the {MODEL} block of code?
Ah ok i see - thank you @stillyslalom, I will keep this in mind!
I forgot to mention, once you load the .tec file into ParaView you will want to plot it using the Wireframe attribute.
Yes, if you want a mesh of square Cartesian elements you can delete the {MODEL} block from the control file and only pass the {RUN_PARAMETERS}. This will generate the mesh but will also throw a spurious warning about an uninitialized dictionary. We have since fixed the HOHQMesh Fortran source to avoid throwing this unnecessary warning, but have not created a new release.
For reference, here is the output I get when removing the model block and generating a Cartesian box mesh:
------------------------------------------------------------------
The following errors were found when constructing the project:
Exception Named: Warning error exception
{
poster = initWithDictionary
}
{
message = MODEL block is missing from control file
}
------------------------------------------------------------------
*******************
2D Mesh Statistics:
*******************
Total time = 1.5669000000000002E-002
Number of nodes = 289
Number of Edges = 544
Number of Elements = 256
Mesh Quality:
Measure Minimum Maximum Average Acceptable Low Acceptable High Reference
Signed Area 1.00000000 1.00000000 1.00000000 0.00000000 999.99900000 1.00000000
Aspect Ratio 1.00000000 1.00000000 1.00000000 1.00000000 999.99900000 1.00000000
Condition 1.00000000 1.00000000 1.00000000 1.00000000 4.00000000 1.00000000
Edge Ratio 1.00000000 1.00000000 1.00000000 1.00000000 4.00000000 1.00000000
Jacobian 1.00000000 1.00000000 1.00000000 0.00000000 999.99900000 1.00000000
Minimum Angle 90.00000000 90.00000000 90.00000000 40.00000000 90.00000000 90.00000000
Maximum Angle 90.00000000 90.00000000 90.00000000 90.00000000 135.00000000 90.00000000
Area Sign 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000
With the help of the Trixi community and our teacher, we managed to solve the eikonal equation to find a velocity field which was used to update a density field with the nonlinear 2D advection equation. An example of a solution can be seen her: crowding problem animation. Thanks for your help to get at this point!
As you can see in the animation, we solve the eikonal and advection equation on a simple, square grid. We would like to solve it on more complex meshes. For the advection equation, this is simple since we are using Trixi to solve this equation, which allows us to use complexer meshes. However, for the eikonal equation, we now use some Fast-Marching Algorithm package (FEFMM.jl) which only works for rectangular meshes. Therefore, it would be very nice if there existed some way to solve the eikonal equation using Trixi on the same mesh as on which we solve the advection equation. Does someone have an advice on how we could approach this problem?
- In this paper, the eikonal equation is transformed into some form better suited for a discontinuous Galerkin framework. Can this be useful in some way for solving the eikonal equation with
Trixi? - The eikonal equation can also be rewritten in an ODE system using the method of characteristics (as in this paper). Can this help us?
Do you think it is possible at all to solve the eikonal equation on the same mesh as the advection equation using Trixi? Or is it a better idea to consider other Eikonal solvers?
1 Like
Right now, we don’t have second-derivative terms in Trixi.jl yet (see e.g., viscous semidiscretizations · Issue #704 · trixi-framework/Trixi.jl · GitHub, Notes on parabolic terms · Issue #1147 · trixi-framework/Trixi.jl · GitHub) but there is some work in progress. Thus, we can’t discretize eq. (13) of the paper you mentioned at the moment.
I didn’t read the paper. There are a lot of efficient ODE solvers in Julia, look at the SciML ecosystem.