WENO scheme and transport equation


As given in the documentation, the transport equation solved is given by


In order to solve the advection term using WENO scheme on 5 point stencil, to my understanding explicit weno5Z_upwind is present in the code.
I would like to know how this term is discetized? The advection term which reads as :
V.grad(Sj) can be re-formulated as
div ( V Sj) - Sj div (V)

Ques 1: As the code is incompressible, is te last term neglected in the code?
Ques 2: Does the explicit WENO upwind formulation refers to WENO dsicretization of Sj but upwind velocity?

Ques 3: If I wish to add more terms on RHS of equation which are dependent of Sj in addition to already present diffusion term, say grad(Sj) and add the discretization implicitly, then as per my understanding, the coefficients of the matrix needs to be modified. May I know which one for transport equation.

Thanking You


The advection term is re-formulated as you wrote (conservative form). The Sj*div(V) term is kept because in some cases divergence may not be strictly equal to computer precision. You can disable this term in the .nts file, in the modeling block of the transport equation, thanks to the following keyword:


To add implicit terms to the transport equation matrix, I advice you to read the following code :

  • src/lib/modeling/generic_transport_equation/solve.f90
  • src/lib/discretization/implicit_discretization/cell_scalar/discretize_cell_transport_equation.f90
    in which each term is discretized. You have to add here a call to a new routine you will create.

For instance, the diffusion term is discretized into

and 2 others routines that are into the src/lib/discretization/implicit_discretization/cell_scalar/diffusion_term
directory. Coefficients of the diffusion term are added to the matrix in the following routine:


You may copy this routine and adapt it in order to discretize the terms you want to add (same way you are used to do with Thétis code)

Be carefull to the stencil part if your stencil is not a 5 points one. We may talk of this point later.

See you

To answer your 2nd question, the WENO scheme solves:
div( u \cdot S_j )
in the Finite Volume approach with the flux
u \cdot S_j
being discretized at the cell’s faces.

This flux can be treated in several manners (with different flux types and limiters).
By default it is a low order Godunov flux, which is simply an upwind reconstruction scheme of S_j based on the entering/exiting velocity at the face (see src/lib/modeling/species_transport/variables.f90 and src/lib/discretization/explicit_discretization/type_FV_flux.f90). You can change this for higher order (LW, FLIC, or Force) from the NTS file.

The WENO reconstruction (and not interpolation) schemes (“3Z” and “5Z”) are ‘traditional’ ones (and needs to have enough ghost nodes to work correctly)…


As per your reply, I have gone through the subroutine described and to my understanding, I need to make changed in add_cell_diffusion_term_centered_o2 in order to implicitly add higher order terms in the coefficient matrix.
As I could see in the subroutine , the coefficient matrix has been filled as follows
( for example left node)
matrix(ia + i_left)= matrix(ia+i_left) -left
where left = alpha_u(i)*flux_coef%u(i,j,k)/dx(i)/alpha(i).

In the above expression, my query is, if it is second order scheme, then shouldn’t dx(i) be dx(i)**2 ?
Kindly clarify.

if I understood correctly, then WENO scheme solves:

div. (u.Sj) - Sj div(u) for advection of scalar term. Here, you refer to the fact that flux of u.Sj can be treated in several manners and based on velocity on the faces. ( which is taken to be at previous time step).The flux of u.Sj is defined using several option ( as you have mentioned) and in case of WENO (5z), which will require 5 point stencil and thus 3 ghost cell.

May you please describe difference between interpolation and reconstruction as you have highlighted in particular.

It is a second order centered scheme. The dx(i) you expect is hidden into flux_coef variable.
See you

So if I understand correctly, then f_coef is actually f_coef/dx and other dx comes as is written in the code.

flux is defined as D*grad (S) where D is the diffusion coefficient. The flux coefficients, coming from the discretization on faces, can be roughly written as :

flux_coef%u = D/dx
flux_coef%v = D/dy

  1. When using an explicit scheme for advecting a scalar, velocity at times n-1, n and n+1 are known. Then, we use 3rd order interpolation of the velocity field to integrate at the desired intermediate time (for example when using RK schemes or when the code used sub iterations to ensure stability)
  2. The term interpolation is used when dealing with node based values and reconstruction with cell-meaned values used to reconstruct the field at a particular position x:
    \phi_i=\phi(x_i) ===> \tilda{\phi} (x) = interpolation( \phi_i, x )
    \psi_i=\Delta x^{-1} \int_{x_{i-1/2}}^{x_{i+1/2}} \psi(x) ===> (\tilda{\psi} (x) = reconstruction( \psi_i, x )

The interpolation of mean values is different from the reconstruction using mean values when using anything different than a centered 2nd order scheme.