Monotoinic upstream scheme for discretizing the convection term

[ThitiACHR]

Dear BOUT++ community

I am studying the transport model, where the expression of the model can be found at https://iopscience.iop.org/article/10.1088/1741-4326/ad2d39 (i.e., equation 1, 6, 10, A5, and A6). We will focus on the convection term of equations A5 and A6.

Regarding the numerical method in Appendix A.3 of the reference paper, they use the monotonic upstream scheme to discretize the convection term. I have seen that in a document of BOUT++, the convection term can be implemented through the function FDDY(V_n,n) where V_n is the particle pinch velocity that depends on profile evolution and n is the density.
The built-in numerical schemes to discretize this term are upwinding (U1) and central finite difference (C2 and C4) methods. I want to implement the monotonic upstream scheme for discretizing this term. May anybody suggest to me where I should start? So far, unwinding and central finite difference methods are not working properly.

Note that once I have seen the MUSCL scheme is available in BOUT++ version 4.2.3

BOUT-dev/src/mesh/index_derivs.cxx

Line 245 in 9ca15c3

//////////////////////// MUSCL scheme ///////////////////////

, but it disappears after 4.3.0 for some reasons.

[bendudson]

Hi @ThitiACHR thanks for trying BOUT++!
I'm not sure what the paper means by "monotonic upstream schemes": That could mean several things. They are only using 1st order Euler for time integration, so "monotonic upstream" might just mean 1st order upwinding (U1).
If you want a higher order method, then "W3" (WENO3) should also be available, as well as 2nd and 3rd-order upwinding (U2 and U3). The implementations are here in the code: https://github.com/boutproject/BOUT-dev/blob/master/src/mesh/index_derivs.cxx#L194

Many schemes, including MUSCL, are difficult to implement in a simple per-operator way, because they need to know characteristic wave speeds of the system for the approximate Riemann solve. The implementation was probably removed because it didn't work well and wasn't being used.

I think the most advanced methods for BOUT++ are implemented in Hermes-3 (but I'm biased). See shock tube tests: https://hermes3.readthedocs.io/en/latest/tests.html#sod-shock and writeup of the method: https://hermes3.readthedocs.io/en/latest/numerical_methods.html

[ThitiACHR]

Hi @bendudson, thank you for your quick response.

In this case, the pinch velocity (V_n) is not constant where it is position-dependent. Suppose we consider the system of PDEs in the y-direction for simplicity, the expression of this term in equation A6 is ddy(V_n*n). So, in the physics model, I should implement this term as FDDY(V_n, n) not VDDY(V_n, n), where the latter treats the derivative as V_n*ddy(n). It looks like the code you provided me is the derivative methods (U1, U2, U3, and W3) for VDDY not the FDDY. In the code of the FDDY, there are only U1, U2, C2, and C4 available. Please correct this if I understand incorrectly.

Why, if the time integrator is 1st order Euler, the monotonic upstream scheme might be 1st order unwinding. Is there any reason behind this?

I will take a look at the sod-shock test in Hermes3.

[ThitiACHR]

Hi @bendudson

1. Are there any ways to implement a new finite volume method in BOUT++? I want to implement something like the Lax-Wendroff method and Beam-Warimng method to study constant-coefficient advection equation, even though the finite volume method provided in BOUT++ is better in accuracy and can handle the numerical oscillation using the reconstruction with slope limiter.

   BOUT-dev/include/bout/fv_ops.hxx

   Line 192 in c193e73

   const Field3D Div_par(const Field3D& f_in, const Field3D& v_in,

My final goal is to implement a high-resolution method such as MUSCL for nonlinear problems like those in the aforementioned paper by Xiaobo Li et al., 2024. It would be better if I start with something simple that can be used as a starting point to implement higher-resolution schemes.

2. In the context of nonlinear conservation law, where the flux function: f(q) is a nonlinear function of quantity (q) when using the FV::Div_par operator to discretize the convection term, how can we obtain wave_speed_in or a, which is the input of the operator. I guess it is given by the derivative of the flux function with respect to q i.e. a = f'(q) where q is the density of conserved quantity.