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

Package Announcements
, , , ,
62

How about splitting second order in time into two coupled first order in time equations? Thx.

63

It’s second-order in space, not time, but the basic idea behind Incorporate parabolic terms into `main` by jlchan · Pull Request #1149 · trixi-framework/Trixi.jl · GitHub. We just need more time to polish it.

64

Expand second order in space as div ( grad() ) ?

65

Yes

66

Is there currently or plans in the future for multi material support?

67

Do you mean multiphase, or multicomponent? Trixi already supports multiple components for the Euler equation (limited to an ideal gas EOS) [ref]

68

I meant multicomponent so thanks for the references!

69

Hey! Trixi uses DifferentialEquations.jl to solve the system of (semidiscretized) ODEs. We were wondering if Trixi also supports solving a PDE for steady state, using one of the Steady State Solvers which DifferentialEquations.jl offers (see this page). Thanks in advance!

EDIT: just tried this and this seems to work. Now, suppose the steady state solver works and we want to solve the following transport equation with source term for steady state:

\partial_t s + \partial_x(s \, q_1) + \partial_y (s \, q_2) = \frac{1}{f^2} + s(\partial_x q_1 + \partial_y q_2)

How would you advise on approximating the partial derivatives \partial_x q_1, \partial_y q_2 efficiently and how can we add the right hand side term as a source term in the Trixi environment?

70

What is required to extend InviscidBurgersEquation1D() to a 2D equivalent?

EDIT: a guide in defining UserDefinedEquation() would be valuable to have.

71

Depends on your q1, q2. If they are known analytically, just implement them like that. Otherwise, fake variable coefficients as in our tutorial 8 Adding a non-conservative equation · Trixi.jl.

72

The tutorial 7 Adding a new scalar conservation law · Trixi.jl describes the 1D setting. If you go to 2D, you need to set the number of spatial dimensions in the subtyping appropriately. Then, you need to define

Trixi.flux(u, orientation, equations)

where orientation == 1 for the x direction and orientation == 1 for the y direction. See for example Trixi.jl/elixir_kpp.jl at main · trixi-framework/Trixi.jl · GitHub (linked from the 1D tutorial).

Do you have specific suggestions for improving this part of our docs/tutorials? If so, could you please create a pull request (or open an issue)?

73

Thx! Will go over it and get back to you. Cheers.

74

Thank you very much for your fast reaction and for your useful answers! Our q1, q2 are not known analytically, we have to estimate them from the data. If we try implement the 2D Inviscid Burgers equations

\frac{\partial q_1}{\partial t} + q_1\frac{\partial q_1}{\partial x} + q_2 \frac{\partial q_1}{\partial y} = 0
\frac{\partial q_2}{\partial t} + q_1\frac{\partial q_2}{\partial x} + q_2 \frac{\partial q_2}{\partial y} = 0

following the implementation in tutorial 8, we get the following structure:

struct BurgersEquations2D <: AbstractEquations{2, 2} #(spatial dimensions, variables)
end

The conservative flux part in these equations is 0:

flux(u, orientation, equation::BurgersEquations2D) = zero(u)

and the non-conservative flux is defined as follows.

# This "nonconservative numerical flux" implements the nonconservative terms.
# In general, nonconservative terms can be written in the form
#   g(u) ∂ₓ h(u)
# Thus, a discrete difference approximation of this nonconservative term needs
# - `u mine`:  the value of `u` at the current position (for g(u))
# - `u_other`: the values of `u` in a neighborhood of the current position (for ∂ₓ h(u))

function flux_nonconservative(u_mine, u_other, orientation, equations::BurgersEquations2D)
    if orientation == 1
        q1_mine,  q2_mine  = u_mine
        q1_other, q2_other = u_other
        f = SVector(q1_mine * q1_other, q1_mine * q2_other)
    else orientation == 2 
        q1_mine,  q2_mine  = u_mine
        q1_other, q2_other = u_other
        f = SVector(q2_mine * q1_other, q2_mine * q2_other)
    end
    return f
end

We are not really certain about the definition of the non-conservative flux here: did we implement this in a correct way? Where in this code are the partial derivatives \partial_x(q_1) and \partial_y(q_2) estimated? Is this already done in the definition of u_other? Or should we do this ourselves in the definition of the non-conservative flux (and what would you recommend to do so)?

EDIT: we found out how to implement Dirichlet Boundary conditions, so removed that part of the question.

75

Looks good to me (but you should try it out with a toy example to be sure).

Not directly in this code, but it will be called in loops described in the tutorial:

du_m(D, u) = sum(D[m, l] * flux_nonconservative(u[m], u[l], 1, equations)) # orientation 1: x

Here, D is the derivative matrix. Thus, the derivative at index m is obtained by summing over l with fluxes where the first argument is fixed to u[m] and the second one varies with u[l]. Thus, the way you have written the non-conservative flux ensures that you take the derivatives of q1_other and q2_other.

76

Hey there, I am also one of ziolai’s student. We currently have insufficient information about the TreeMesh. I have two questions about it and I hope you guys can help me out/ point me to the right direction.

  1. How do we extract/calculate the average solution over a quadrilateral element?
  2. Once we know a new average over a quadrilateral element, how do we update the internal nodes?

Some more background information:
For each timestep we currently modify our rhs before letting trixi do it’s magic. For this we first obtain the average density of crowds over each rectangular element by computing a 2D spline over the whole domain and then integrating over each rectangular element (our first mistake). Then we apply an algorithm to determine the direction of the crowd for each element. And lastly, we update the direction for each node used in our current mesh by interpolation and extrapolation (our second mistake).

We constructed our mesh as

TreeMesh((0, 0), (64, 64), initial_refinement_level=6, n_cells_max=10000, periodicity = false)

which constructs a 64x64 uniform square mesh. However each square element contains 16 internal nodes and we do not know how to extract and import data to each of the internal nodes.

77

See for example

Note that this uses undocumented internal API, so be aware…

You want to change the position of the integration nodes? You shouldn’t do that…

78

Thanks for your feedback. So I can calculate the average values now.

I absolutely do not want to change the position of the integration nodes. I want to update the values of the solution over a quadrilateral element, so I also need to update the values on the integration nodes.

79

Ah, I see. Just scroll a bit further down in the code I linked above, to

This sets the nodal variables in the vector u at index (i, j) in the given element to theta * u_node + (1-theta) * u_mean, a convex combination of the current node variables and the mean value. I guess that should help you in your case?

1 Like
80

Yes, this is exactly what I needed. Thanks!

81

Is there a way to obtain the spatial partial derivatives of a computed field in post-processing stage?