How to implement the TVD scheme correctly?

[mhvwerts]

When updating the PyFVTool-introduction-demo Notebook, I came across the question on how to correctly implement the TVD scheme in solving convective transport.

My present version of the Notebook is here (on the solvePDE_notebooks branch of my PyFVTool fork):

https://github.com/mhvwerts/PyFVTool/blob/solvePDE_notebooks/examples-notebooks/PyFVTool-introduction-demo.ipynb

The proposed TVD scheme is in cell number 52.

Two questions:

1. The pf.convectionTvdRHSTerm() function only provides a RHS term. What should be used as the convective matrix term? I guessed that it could be the upwind matrix term, provided by pf.convectionUpwindTerm(). Is this correct? Are there perhaps several possibilities?
2. There is an inner loop for the TVD scheme, which is looped over 5 times. Why 5? Does 5 always work fine? Can it be optimized? Perhaps there are some convergence criteria that may be applied to the inner loop. We do not have to do this, but we may mention it in the Notebook text.

Many thanks for your help. I will document the answers in the Notebook.

---

Excerpt from the present (proposed) TVD code

# Constant coefficient matrices:
Mconvuw = pf.convectionUpwindTerm(uf)   # upwind convection term

# Choose a flux limiter
FL = pf.fluxLimiter('SUPERBEE')

# Solver
dt = 0.001  # [s], time step
final_t = W/u

tt = 0.0
count = 0

phi_tvd = []
phi_uw = []
while (tt<final_t) and (count < 5000):
    tt += dt

    # inner loop for TVD scheme
    Mt, RHSt = pf.transientTerm(phi_old, dt, alfa) # this remains constant over the inner loop
    for jj in range(5):   # Why 5?
        # TVD RHS term
        RHSconv = pf.convectionTvdRHSTerm(uf, phi, FL)

        # Form the matrix equation
        # MW: It is not clear how to construct the Mconv term in the TVD scheme:
        #     convectionTvdRHSTerm only provides RHS. How to obtain Mconv?
        # M = Mconv + Mt + Mbc  # This is an error since Mconv undefined
        # M = Mt + Mbc          # There should be a Mconv term, I guess...
        M = Mconvuw + Mt + Mbc  # Use upwind matrix for Mconv in TVD scheme?
        RHS = RHSt + RHSbc + RHSconv # Combined with the RHS provided by convectionTvdRHSTerm.

        # Solve the PDE
        phi = pf.solveMatrixPDE(mesh1, M, RHS)

    # Store the TVD solution for later animation
    phi_tvd.append(phi)

    # Calculate the tarnsient term for the upwind scheme
    Mtuw, RHStuw = pf.transientTerm(phiuw_old, dt, alfa);

    # Form the system of equations
    Muw = Mconvuw + Mtuw + Mbc
    RHSuw = RHStuw + RHSbc

    # Solve the system of equations using the Upwind discretization    
    phiuw = pf.solveMatrixPDE(mesh1, Muw, RHSuw)

    # store the solution 
    phiuw_old = phiuw
    phi_old = phi
    phi_uw.append(phiuw)

    count += 1

[simulkade]

The TVD scheme I implemented here provides an 'explicit' correction term for the first degree upwind based on the values of phi in the previous time step. At the beginning, I had implemented a function convectionTVDTerm that returned both the matrix of coefficients and the RHS corrector. Soon I realized that it is wasteful, since the matrix of coefficient comes from upwind and does not change in the inner 'corrector' loop (if at all!). Now, we have the pf.convectionTvdRHSTerm() that returns only the corrector term. Your example is therefore correct.
Strictly speaking, we have to use a convergence criterion for the inner loop. However, I consulted some other codes (written for flow in porous media) and noticed that they just iterate a few times, usually 3. I must admit that the flow in porous media community uses these tricks all the time simply because they work for their specific problem (often multiphase flow in porous media). It makes a lot of sense there because there is another loop for solving the nonlinear equation and it adds another time-consuming step to the solver. However, it does not mean that this trick always works.

[mhvwerts]

Thanks for these explanations.

I have cleaned up the part on the different discretization schemes for convection, switching also to 'new style' solvePDE() which gives results identical to the original solveMatrixPDE() version. I also added the 'standard' central difference convection term to the comparison.

The document still needs some further editing, which I hope to finish soon.

P.S. I also tried 3 TVD iterations instead of 5. This seems to work fine as well in this case, and speeds up the calculation.

[simulkade]

I also tried 3 TVD iterations instead of 5. This seems to work fine as well in this case, and speeds up the calculation.

I guess ultimately it depends on the problem, and how much one cares about the numerical diffusion. This simple TVD scheme does not replace proper shock-capturing methods but does a pretty decent job at capturing some features in the solution by eliminating some of the numerical diffusion. I have a couple of beautiful examples of multiphase flow in porous media with analytical fractional flow solutions, e.g. this one. Hopefully, I can add them as advanced examples to pyfvtool.

[mhvwerts]

It would be excellent to have modeling examples of multiphase flow in porous media with PyFVTool. This may also help in getting perhaps some geophysics and hydrogeology colleagues and/or their students interested in the project.

I see that you programmed your example in Julia, and I am sorry about pushing Python. Julia is clearly the language of the future for scientific computing, but all the kids (present and future students) nowadays learn Python in school.

[mhvwerts]

So, thanks again for making the FVTool available in Python!

[simulkade]

It would be excellent to have modeling examples of multiphase flow in porous media with PyFVTool. This may also help in getting perhaps some geophysics and hydrogeology colleagues and/or their students interested in the project.

I see that you programmed your example in Julia, and I am sorry about pushing Python. Julia is clearly the language of the future for scientific computing, but all the kids (present and future students) nowadays learn Python in school.

I am actually quite happy with Python these days, especially because of the more mature ecosystem of packages. It is kind of fun to program in Julia but I must admit debugging can be a struggle, perhaps because I'm not good enough at it! And the packages are changing way too fast so for a casual user like me it was a bit difficult to keep up with the great developments happening in Julia community.