[Jaydi_21]

I'm working with the 2d shallow water equations, but my question is in a general sense. I'm not sure if I understood how a TVD scheme is implemented. Let's say we have a quadrilateral cell with two boundaries in the horizontal (x) direction and two boundaries in the vertical direction (y). This cell has neighbors, so I understand that TVD schemes intend to calculate the values of the variables at the face of the cells and now solve the Riemann problem for the interpolated values at each side of the cell. Is this correct?

Second, if my first assumption is correct. How do you calculate the values of the variables on both sides of the face? I've seen that you may use the cells you have before and next to the cell you are dealing with, but that only allows you to find ONE gradient. Am I suppose to use the same gradient for both faces of the cell.


For example, if I want to solve the RP at the yellow face, so I use the blue cell and the red cell, together with the one in the middle (black) and that allows me to compute a gradient. But only one gradient, then I can use the black cell and this gradient to calculate the value at the left side of the yellow face. But what happens with the right side of this face? Am I supposed to repeat the procedure with another two neighbors of the blue cell?

https://www.dropbox.com/s/mgamo1mux4...ptura.JPG?dl=0

[t.teschner]

I think there are some confusion about (TVD) schemes in general, maybe an explanation from a different angle may be helpful.

First of all, TVD (Total Variation Diminishing) is just doing what its name says (although this may not be helpful if you still are unsure how they work). Let's put into an practical example: Say you have two cells C0 and C1 which are connected by a face f. In the finite volume method, we want to obtain the values on the faces inside a cell from which we then reconstruct our variables at the cell center (if we use a cell center scheme, which I assume here). So we have to get our values somehow onto the face f from information that we have stored at C0's and C1's center. There are two easy approaches, one is using a first order upwind or a second order central scheme. Say we have for some variable phi a value of 0.5 at C0's center and 1.0 and C1's center. If we were to use the first order scheme (not upwind (yet), just first order for now), we simply copy the values that we have at the cell centers to the cell face. So this means we now have two values at the face, once 0.5 and once 1.0. You could use a Riemann solver here as you pointed out which just consolidates both values into one, based on some "magical math", lets say (to keep the discussion simple). Lets say the Riemann solver gives you a value of 0.70, so you replace your value at the face on both side by this value. This approach works well if you can use a Riemann solver (which you only can for hyperbolic systems, as you require the eigenvalue of the system of equations you are solving), but what if the Riemann solver approach is not available to you? Well, you can make your first order scheme an upwind scheme by simply checking from which direction the flow is coming and then selecting the value which is upwind of that flow direction. So if the flow is from C0 to C1, we select the value of 0.5 at f, and 1.0 if the flow is the reverse.

The second approach is to use a central scheme, where you simply average the values. Let's assume the cells are equal in size and f is directly in the middle of both cell centers, then you get a value of (0.5 + 1.0) / 2 = 0.75. Simple enough. The problem here is that the central scheme does not take into account where the flow is coming from. This will not allow to simulate flows with high Reynolds numbers (actually, Peclet numbers, but they are somewhat related and the Reynolds number is easier to visualise). So, you might say then lets use an upwind approach again, but make it second order.

At this point is where it all comes together, TVD schemes, gradients and a very famous theorem from Godunov which states that there are no schemes of second order or higher, which do not generate new extrema. What does that mean? Well, lets look at gradients for a moment. We compute these, for example with the Gauss theorem by looking at all the fluxes going in and out of a cell and dividing the cell by its volume. So each cell will have their own gradient in fact, so it is not shared by two cells. Once you have both gradients in cell C0 and C1, you construct the value at the face f by adding the value at the cell center plus the gradient at the cell center multiplied by the distance to the face. so for cell C0 you have phi_f,L = phi_C0 + grad_C0 * r_C0_f and for C1 you have similarly phi_f,R = phi_C1 + grad_C1 * r_C1_f, where r is the distance from the cell center of C0 or C1 to f. If we were to upwind those values, we would again look from where the flow is coming from and then make a decision which value to take. Lets say the flow is from C0 to C1, so we use the information from C0 to compute our value at the face. Lets say r = 1.0 for simplicity and a value of grad_C0 = 0.2. Then we have for the value of phi at the face f 0.5 + 0.2 * 1 = 0.7. So far so good, but now lets look what happens when grad_C1 becomes larger, say 0.6. Then you would have 0.5 + 0.6 * 1.0 = 1.1. This is a problem! Your value at C0 for phi is 0.5 and at C1 is 1.0. Yet, at the face you have generate a value of 1.1, so it is not between 0.5 and 1.0, which means we have created a new extrema. This is exactly why we use TVD schemes, they are used to circumvent Godunov's theorem by forcing the solution to not created new extrema, they limit your reconstructed values. So in our case, a TVD scheme would also provide a limiter to our second order upwind scheme and limit the result so that we are still bounded. If the limiter would see that we create a new extrema, as we did with 1.1, then it would simply limit that to 1.0. That's it, that's all a TVD scheme (rather, its limiter) is doing for you.

Hope this helped.

[FMDenaro]

Just to add that the Godunov theorem stated that there are "no linear scheme higher than first order", indeed non linear schemes (used also for linear equations) are the basis of TVD idea.

[Jaydi_21]

Originally Posted by t.teschner

I think there are some confusion about (TVD) schemes in general, maybe an explanation from a different angle may be helpful.

First of all, TVD (Total Variation Diminishing) is just doing what its name says (although this may not be helpful if you still are unsure how they work). Let's put into an practical example: Say you have two cells C0 and C1 which are connected by a face f. In the finite volume method, we want to obtain the values on the faces inside a cell from which we then reconstruct our variables at the cell center (if we use a cell center scheme, which I assume here). So we have to get our values somehow onto the face f from information that we have stored at C0's and C1's center. There are two easy approaches, one is using a first order upwind or a second order central scheme. Say we have for some variable phi a value of 0.5 at C0's center and 1.0 and C1's center. If we were to use the first order scheme (not upwind (yet), just first order for now), we simply copy the values that we have at the cell centers to the cell face. So this means we now have two values at the face, once 0.5 and once 1.0. You could use a Riemann solver here as you pointed out which just consolidates both values into one, based on some "magical math", lets say (to keep the discussion simple). Lets say the Riemann solver gives you a value of 0.70, so you replace your value at the face on both side by this value. This approach works well if you can use a Riemann solver (which you only can for hyperbolic systems, as you require the eigenvalue of the system of equations you are solving), but what if the Riemann solver approach is not available to you? Well, you can make your first order scheme an upwind scheme by simply checking from which direction the flow is coming and then selecting the value which is upwind of that flow direction. So if the flow is from C0 to C1, we select the value of 0.5 at f, and 1.0 if the flow is the reverse.

The second approach is to use a central scheme, where you simply average the values. Let's assume the cells are equal in size and f is directly in the middle of both cell centers, then you get a value of (0.5 + 1.0) / 2 = 0.75. Simple enough. The problem here is that the central scheme does not take into account where the flow is coming from. This will not allow to simulate flows with high Reynolds numbers (actually, Peclet numbers, but they are somewhat related and the Reynolds number is easier to visualise). So, you might say then lets use an upwind approach again, but make it second order.

At this point is where it all comes together, TVD schemes, gradients and a very famous theorem from Godunov which states that there are no schemes of second order or higher, which do not generate new extrema. What does that mean? Well, lets look at gradients for a moment. We compute these, for example with the Gauss theorem by looking at all the fluxes going in and out of a cell and dividing the cell by its volume. So each cell will have their own gradient in fact, so it is not shared by two cells. Once you have both gradients in cell C0 and C1, you construct the value at the face f by adding the value at the cell center plus the gradient at the cell center multiplied by the distance to the face. so for cell C0 you have phi_f,L = phi_C0 + grad_C0 * r_C0_f and for C1 you have similarly phi_f,R = phi_C1 + grad_C1 * r_C1_f, where r is the distance from the cell center of C0 or C1 to f. If we were to upwind those values, we would again look from where the flow is coming from and then make a decision which value to take. Lets say the flow is from C0 to C1, so we use the information from C0 to compute our value at the face. Lets say r = 1.0 for simplicity and a value of grad_C0 = 0.2. Then we have for the value of phi at the face f 0.5 + 0.2 * 1 = 0.7. So far so good, but now lets look what happens when grad_C1 becomes larger, say 0.6. Then you would have 0.5 + 0.6 * 1.0 = 1.1. This is a problem! Your value at C0 for phi is 0.5 and at C1 is 1.0. Yet, at the face you have generate a value of 1.1, so it is not between 0.5 and 1.0, which means we have created a new extrema. This is exactly why we use TVD schemes, they are used to circumvent Godunov's theorem by forcing the solution to not created new extrema, they limit your reconstructed values. So in our case, a TVD scheme would also provide a limiter to our second order upwind scheme and limit the result so that we are still bounded. If the limiter would see that we create a new extrema, as we did with 1.1, then it would simply limit that to 1.0. That's it, that's all a TVD scheme (rather, its limiter) is doing for you.

Hope this helped.

Thanks for your good explanation. So I think I understood...I find the gradient for each cell and then I reconstruct the variables at both sides of a given face (with the cell center information). This means that I need to use the neighbors of each cell, to find a gradient for that cell. Finally, I use the Riemann solver, which is a Godunov upwind scheme, with the reconstructed variables.
One more question, does it matter the direction of the gradient? Because I have a two dimensional hyperbolic system in a non uniform cartesian grid. I have quadrilateral cells and what I did is to calculate the values (from IDW interpolation) at the middle of each face and with this values I calculated a gradient for each cell (a different gradient for the two relevant directions, let's say "x" and "y"), but the middle of each face on both sides of a given cell, do not necessarily make a straight line. And the resulting scheme is not converging...what do you think? (It's a hierarchical grid, quadtree type).

[Jaydi_21]

Originally Posted by Jaydi_21

Thanks for your good explanation. So I think I understood...I find the gradient for each cell and then I reconstruct the variables at both sides of a given face (with the cell center information). This means that I need to use the neighbors of each cell, to find a gradient for that cell. Finally, I use the Riemann solver, which is a Godunov upwind scheme, with the reconstructed variables.
One more question, does it matter the direction of the gradient? Because I have a two dimensional hyperbolic system in a non uniform cartesian grid. I have quadrilateral cells and what I did is to calculate the values (from IDW interpolation) at the middle of each face and with this values I calculated a gradient for each cell (a different gradient for the two relevant directions, let's say "x" and "y"), but the middle of each face on both sides of a given cell, do not necessarily make a straight line. And the resulting scheme is not converging...what do you think? (It's a hierarchical grid, quadtree type).

I don't know if my question now is clear, so I want to add this. Since I'm working with a 2d set of equations, I'm not sure if the way I'm calculating the gradients is correct. I have a portion of my mesh...in the figure.

https://www.dropbox.com/s/7c3snrucx507v4t/mesh.JPG?dl=0

For example for the cell with the red "x", I used the interpolated values at the middle of the blue faces, to calculate a gradient for the y-direction and similarly, the green faces for the x-direction. As you can see, the blue faces have a gradient in the x-direction (I took them into account also), so I calculated two "gradients", one for each direction. To calculate this gradient, I used the distance, but you can see they don´t have the same direction and then, with this gradient, I calculated again the values at the middle of each cell (using only the distance along the axis I was considering). So I used some sort of directional derivative. But now the scheme does not converge so I'm wondering how many possibilities to make this I have...A LOT.

Should I calculate a gradient using all the neighbors of a given cell, or a gradient for each face, using the neighbors of this face and its extreme points? Once I have the gradient, should I calculate the value at the middle point of the face using the distance("r" from the cell center to the face center) or only the component in the x- or y-direction? In the second case, then, for the example of the figure, the value at both green faces on the left side of the cell would be the same.

I'm kind of confused...as you may see!

PS. I limited the reconstructed values so that it didn't create new local extrema, so I believe that the problem is the way I'm calculating the gradients and the direction of them.

[FMDenaro]

I suggest to think always in terms of the integral of the flux along the surface. This way you can consider to summ the contributions from each piece of the surface.



I don't know if you are aware of this well known work, you could start from this approach.


https://foam-extend.fsb.hr/wp-conten...k_PhD_1996.pdf

[TurbJet]

Originally Posted by FMDenaro

Just to add that the Godunov theorem stated that there are "no linear scheme higher than first order", indeed non linear schemes (used also for linear equations) are the basis of TVD idea.

This is where confuses me: How do ENO/WENO-like schemes count? Apparently they can achieve really high order. But, do they count as TVD? Do they count as linear scheme?

[FMDenaro]

Originally Posted by TurbJet

This is where confuses me: How do ENO/WENO-like schemes count? Apparently they can achieve really high order. But, do they count as TVD? Do they count as linear scheme?

They are non-linear schemes even if you solve a linear problem

[TurbJet]

Originally Posted by FMDenaro

They are non-linear schemes even if you solve a linear problem

I see. But what distinguishes linear/nonlinear scheme?

[FMDenaro]

Originally Posted by TurbJet

I see. But what distinguishes linear/nonlinear scheme?

Exactly as it for ODE/PDE, the coefficients depend on the solution itself.

As example, just consider the action of the limiter that introduces a coefficient that depends on the local field.

[TurbJet]

Originally Posted by FMDenaro

Exactly as it for ODE/PDE, the coefficients depend on the solution itself.

As example, just consider the action of the limiter that introduces a coefficient that depends on the local field.

Oh, I see. Thx.