[amine.hanini]

Hello,

I am trying to implement in 2D with matlab the least square method to calculate the gradient of my variables to reach a second order accuracy in space. Together with the minmod limiter to avoid oscillations and the appearance of new extrema. My problem here is that I have the feeling that my second order is not really second order and the results is very similar to the first order .
I appreciate any help.

[sbaffini]

Well, the LSQ method is actually only first order for general unstructured grids, and not even always if it doesn't employ any weighting. But on fully structured and cartesian grids you should achieve 2nd order indeed.

The limiter has nothing to do with the gradient accuracy.

[sbaffini]

If you give more details on what you're doing and how are you checking your accuracy, maybe we can give more help

[amine.hanini]

@Sbaffini, Thank you for your response.

Concerning the type of LSQ method, I am using the weighted averaging procedure to compute the values on the vertices of each triangular cell, this computation needs the information of the surrounding cells with common vertex. After the cell node variables are computed, the gradient on control volume can be calculated via some arithmetic average on each cell face. I am supposed to be more accurate with this node-based reconstruction and this bigger stencil than the cell-based one, but it appears to that there is no change between 1st and 2nd order and when I plot the two results, It's exactly the same.

P.S. To limit gradient on each control volume, I use the Barth-Jespersen slope limiter.

[amine.hanini]

untitled.png
LSQ1.jpg
LSQ2.jpg

In these attachement the method the reconstruction used in my matlab code and the results of a simple dam-break and a superposed water height level with the two orders.

[sbaffini]

Originally Posted by amine.hanini

@Sbaffini, Thank you for your response.

Concerning the type of LSQ method, I am using the weighted averaging procedure to compute the values on the vertices of each triangular cell, this computation needs the information of the surrounding cells with common vertex. After the cell node variables are computed, the gradient on control volume can be calculated via some arithmetic average on each cell face. I am supposed to be more accurate with this node-based reconstruction and this bigger stencil than the cell-based one, but it appears to that there is no change between 1st and 2nd order and when I plot the two results, It's exactly the same.

P.S. To limit gradient on each control volume, I use the Barth-Jespersen slope limiter.

So, you are using a Green-Gauss gradient with face values computed as "some average" of node values, which in turn are computed from a WLSQ
interpolation? A variant of this scheme is also available in Fluent as node-based green-gauss gradient.

I have no experience with this approach but, I expect "some average" to actually have a large influence on the final result.

Moreover, on general unstructured grids, you need higher than 2nd order accuracy on face values to actually hope to get higher than 1st order accuracy on gradients. You can see this by noting that the final green-gauss formula is a sum of face values multiplied by area and divided by volume, so if your face values are like

phi_f = phi_f_exact + O(dx^2)

the final gradient will be something like (not considering other error sources):

grad_phi = grad_phi_exact + O(Area*dx^2/Vol) = grad_phi_exact + O(dx)

If your grid has some regularity, it might happen that O(dx^2) terms on opposite faces would cancel each other, but this doesn't happen in general.

Gradient tests should always be done on a set of grids:

1) Uniform cartesian
2) Non-uniform cartesian
3) Prismatic (split the uniform cartesian cells along the diagonal of one plane only)
4) Unstructured (both regular and not, if possible)

and on a set of functions:

1) linear in 1D
2) linear in multi dimensions
3) mixed linear-higher order (e.g., x*y^2)
4) polynomial in multi-dimensions (e.g., x^2+y^2 + (x*y)^2)
5) full non linear ones (sin, cos, exp, etc.)

But in no case you should include the limiter in this exercise.

[sbaffini]

Originally Posted by amine.hanini

Attachment 74724
Attachment 74725
Attachment 74726

In these attachement the method the reconstruction used in my matlab code and the results of a simple dam-break and a superposed water height level with the two orders.

I just posted few seconds after you posted your last message. Few things:

1) Yes, your gradient method is indeed the node-based green-gauss one. Also note that it is mentioned as more robust and not more accurate.

2) You seem to be comparing the two methods on the final full solution (of a dam break problem!!!), which explains why you keep mentioning the limiter. While this might make sense, it has not if you limit the discussion to the gradient. It is highly probable that in your case the gradient has no such a great difference.

3) Are you using a specific code of someone else or something coded by you? In the last case, I suggest to test the gradient alone as mentioned in my last post (if this is what you actually want to test).

[amine.hanini]

Originally Posted by sbaffini

3) Are you using a specific code of someone else or something coded by you? In the last case, I suggest to test the gradient alone as mentioned in my last post (if this is what you actually want to test).

I am using my own matlab code. It's a 2D unstructured triangular code.

Yes, It will be very interesting if I could test this reconstruction on a simple test cases, but I can't follow your steps in your last post; in fact uniform cartesian, non-uniform one and prismatic can not be reached by my code, 1D either.

I can test my reconstruction method using no regular triangular grid.

[QUOTE]

1) linear in 1D
2) linear in multi dimensions
3) mixed linear-higher order (e.g., x*y^2)
4) polynomial in multi-dimensions (e.g., x^2+y^2 + (x*y)^2)
5) full non linear ones (sin, cos, exp, etc.)

[QUOTE]

1) not possible with my code

If I understand well, I have to change my system of equation by a simple linear one, first, like transport equation ?

[FMDenaro]

From a general theoretical point of view, a linear reconstruction of the function can never produce a second order accurate gradient. The only case that happens is for a uniform grid where you have some terms in the local truncation error disappearing.
Thus, the question is what do you want more by doing a LSR?

[sbaffini]

I assumed you were using a 3D code, but you can do almost all the things I mentioned in 2D as well.

You certainly have an initialization part in your code, and this certainly precedes the one where the gradient is computed for the first time. So, the idea to test the gradient alone involves:

1) Initializing the variables (those whose gradient will be computed) with some known function, taking care of correctly applying the boundary conditions to them (using the same known function in a dirichlet bc or neaumann bc... better to test both implementations actually).

2) Once their gradient is computed you can compute the error with respect to the known one (because you know the function used to initialize the fields, so you can manually take their analytical derivative in the two directions x and y), take some norm of it (I suggest both L2 and L_inf) and print it

3) Just stop the code because you don't need anything else from it

You should repeat the above 3 steps, carefully annotating the resulting error norms, for different classes of functions and different grids, and for each combination of the former using different refinement levels for the given type of grid.

For what concerns the possible functions to use as initialization, you could use:

1) k*x (linear in 1D), you choose k
2) k*x*y, m*x + n*y (linear in multidimensions), you choose k, m and n
.
.
.

etc. with the remaining ones I already mentioned in my previous post. That is, you can certainly initialize your fields with whatever you want.

For what concerns the types of grid, some of the cases I mentioned have indeed no sense in 2D. So you could just do:

1) Uniform cartesian
2) Non uniform cartesian
3) Regular triangular (obtained by splitting the quads of a cartesian along the diagonal)
4) Fully triangular with no special properties

If you can't do cartesian grids, well, just do the kind of grids you can actually use

[arkie87]

Originally Posted by sbaffini

The limiter has nothing to do with the gradient accuracy.

I thought that when the limiter is engaged, accuracy approaches first order for all limiters? Isn't that how limiters work?

[amine.hanini]

Originally Posted by FMDenaro

From a general theoretical point of view, a linear reconstruction of the function can never produce a second order accurate gradient. The only case that happens is for a uniform grid where you have some terms in the local truncation error disappearing.
Thus, the question is what do you want more by doing a LSR?

@FMDenaro, thank you for your reponse.

If I can't obtain second order, I must be a little more than 1st order, say 1.2 or 1.3, for instance, you don't think so ?

Can you guide me to any reference on that topic from a general theoritical point of view ?

And is there any other method to compute the gradient that allow me to be arround 2nd order with my configuration ?

Many thanks

[sbaffini]

Originally Posted by arkie87

I thought that when the limiter is engaged, accuracy approaches first order for all limiters? Isn't that how limiters work?

(Gradient) Limiters are a tool used, specifically, to avoid that variables reconstructed on faces trough gradients get out of bounds. So, you have a gradient and, after the limiter is applied to it, you also have the limited gradient.

You are right that, when the limiter is on, the limited gradient is of reduced accuracy (so that reconstructed values end up being first order) but it has no sense to test the accuracy of the limited gradient. Gradient accuracy is tested with unlimited gradients.

[arkie87]

Originally Posted by sbaffini

(Gradient) Limiters are a tool used, specifically, to avoid that variables reconstructed on faces trough gradients get out of bounds. So, you have a gradient and, after the limiter is applied to it, you also have the limited gradient.

You are right that, when the limiter is on, the limited gradient is of reduced accuracy (so that reconstructed values end up being first order) but it has no sense to test the accuracy of the limited gradient. Gradient accuracy is tested with unlimited gradients.

Makes sense.

Is it possible though that OP is experiencing first-order solution with 2nd order technique because his gradients are limited? I suppose it would have to be limited in every cell for that to occur...

[amine.hanini]

Originally Posted by arkie87

Makes sense.

Is it possible though that OP is experiencing first-order solution with 2nd order technique because his gradients are limited? I suppose it would have to be limited in every cell for that to occur...

Yes I think in every cell and at every time step I am limited and I go back to first order. Without limiting I have lot of overshootings

[sbaffini]

It depends from a lot of factors, that's why I suggested him testing the gradient separately if that is what he actually wants to do.

On a separate note, the method he is using is not 2nd order, no matter what.

[sbaffini]

But how do you know that your two gradient methods actually work and are doing what they are expected to do and you don't have bugs in them? Just out of curiosity

[arkie87]

Originally Posted by sbaffini

But how do you know that your two gradient methods actually work and are doing what they are expected to do and you don't have bugs in them? Just out of curiosity

Along the same lines, it might make sense to build a test case (such as the classic case of flow at 45° angle to the square Cartesian grid with bottom at 0K inlet and left at 100K and with thermal conductivity (diffusion) set to 0), and check if it converges as expected. If you select the test case properly, you can ensure you dont have flux limiter engaged for all the points in the domain, so you can eliminate that effect as well.

I actually did something similar (since i wanted to see effect on false/numerical/artificial diffusion), and discovered a bug in my code when i reversed the flow direction. So, I would recommend trying all combinations.

Perhaps sbaffini can suggest a better test case?

[sbaffini]

Originally Posted by arkie87

Along the same lines, it might make sense to build a test case (such as the classic case of flow at 45° angle to the square Cartesian grid with bottom at 0K inlet and left at 100K and with thermal conductivity (diffusion) set to 0), and check if it converges as expected. If you select the test case properly, you can ensure you dont have flux limiter engaged for all the points in the domain, so you can eliminate that effect as well.

I actually did something similar (since i wanted to see effect on false/numerical/artificial diffusion), and discovered a bug in my code when i reversed the flow direction. So, I would recommend trying all combinations.

Perhaps sbaffini can suggest a better test case?

I actually did it in posts n. 6 and 10.


What you suggest is good as well, but comes much later in code verification, it basically is a whole code test.

[FMDenaro]

Just to make clear the issues:


- The flux limiting concept is developed in the framework of monotonic-preserving property, it has nothing to do with the accuracy order. It is applied to prevent numerical oscillations around sharp or infinite gradients while the accuracy order is determined on smooth solutions


- the linear reconstruction, that is a polynomial of degree 1, has a second order accuracy for the approximation of the function but when you make the derivative of the polynomial you get a constant function in the interval between the two values. That is a first order accuracy for the derivative. You cannot improve the accuracy by doing a LSR or some type of averaging.


- Obviously, a quadratic reconstruction (that is a polynomial of degree 2) allows for a second order accurate gradient.


- The accuracy of the gradient enters by means of the diffusive flux therefore you need to test a case in which you have an exact solution from a diffusion problem.