[FMDenaro]

The key is to think always in terms of a totally unstructured grid. You need to extend the topology of the grid, that is each centroid of the faces must be associated to a proper simplex (tetrahedron) containing it. Using the four linear shape function you have all you need.

[arjun]

Originally Posted by sbaffini

Now, to me, there seems that, at least at some point back in the 90's, someone in Fluent was quite sure about this. I still have to try this, and I don't expect great differences (nor convenience in terms of computational costs), but it has a certain consistency. The whole algorithm now requires only face values: those for convection and those in the green gauss procedure for the cross diffusion gradients; both are computed via reconstruction gradients, which are the limited version of the base ones (which in contrast are not anymore needed).

By reading this part i get the impression that you are saying that they are using reconstruction gradients for diffusion terms.

If this is what you meant then I would point out that the quotes from Sanjay Mathur et al. are not saying this.

This is how it shall work:

1. Compute gradients by GG or least square means. Call them unlimited gradients.
2. Limit these gradients - call them reconstruction gradients.
3. Use these gradients and reconstruct face values and compute another gradients by GG. These are the gradients used in diffusion terms. (they are unlimited).

Here are these quote again:

Originally Posted by sbaffini

- The gradients used ( step (3)) in the diffusion terms are computed using the reconstructed face values (step 1 grad) - Mathur and Murthy, Pressure Boundary Conditions for Incompressible Flow Using Unstructured Meshes, Numerical Heat Transfer B, 1997

- The cell derivatives (step 3) used for the secondary diffusion terms in Eq. 12 are computed by again applying the divergence theorem over the control volume and using the averaged reconstructed values (using grad from step 1) at the faces - Mathur and Murthy, A Pressure Based Method for Unstructured Meshes, Numerical Heat Transfer B, 1997

- The cell derivatives needed to evaluate in the secondary diffusion term of equation 23 are computed by applying the Gauss theorem again (using limited grads of (1)) over the cell control volume using the reconstructed face values as follows. - Kim, Mathur, Murthy and Choudhury, A Reynolds-Averaged Navier-Stokes Solver Using Unstructured Mesh-Based Finite-Volume Scheme, AIAA 98-0231

- The gradients used in the diffusion term in Eq. 7 are computed by applying the divergence theorem with the average of reconstructed face values.(using grads of (1))- Mathur and Murthy, All Speed Flows on Unstructured Meshes Using a Pressure Correction Approach, AIAA 99-3365

- When evaluating diffusive fluxes the gradient at each face is computed by first considering an average gradient composed of gradients within the two cells on either side of the face where the cell gradients and are evaluated using the divergence theorem (see Eq. 8), but are different from the reconstruction gradients appearing in Eq. 7 in that they are computed using an average of the reconstructed values and (as opposed to cell-average values), and are not limited. - Weiss, Maruszewski and Smith, Implicit Solution of the Navier-Stokes Equations on Unstructured Meshes, AIAA 97-2103

For your concern about what happens in shock or discontinuity, I would say that once you have gradients of effective viscosity you are no longer bound to using average, you could apply other schemes like upwinding or some other scheme that is used in shock for velocity (what is good enough for goose good enough for gander). The exact same thing is done for density which is normally just used as average or first order upwind.

This all need efforts and cost so people usually avoid it. Seeing all this discussion i would implement this and give user option in wildkatze to switch to different schemes. Who knows to someone it might be useful.

[sbaffini]

@Filippo: Thanks, that's the idea. But, in general, I only have access to the two values across a face, L and R.

Besides linear/harmonic interpolation along the face normal coordinate, or taking a simple/harmonic average, I can at most interpolate, still along the connecting line, in the point nearest to the face center, so that I can minimize the error. But not much more.

@arjun: yes, I probably expressed more obscurely than I tought, but that's what I meant. So you also confirm that there is a third gradient, built via GG and using as face values those reconstructed via limited gradients. This third gradient is NOT itself limited again, but only used in the secondary diffusion term.

@michael: at moment I'm just SDed. Do I still qualify?

[arjun]

Originally Posted by sbaffini

@arjun: yes, I probably expressed more obscurely than I tought, but that's what I meant. So you also confirm that there is a third gradient, built via GG and using as face values those reconstructed via limited gradients. This third gradient is NOT itself limited again, but only used in the secondary diffusion term.

Yes. Not only there but all the other things that need gradient use this gradient. (ie turbulence production etc etc)

I will note down one more thing that many commercial codes when using least square gradients just make a copy of the unlimited gradient and use this (they avoid one more GG application). Then the same gradients are limited and used for reconstruction.

In theory it makes sense but in practice it could be problems in prism type cells so I in wildkatze use extra step of GG and use that as mentioned.

[FMDenaro]

Originally Posted by sbaffini

@Filippo: Thanks, that's the idea. But, in general, I only have access to the two values across a face, L and R.

Besides linear/harmonic interpolation along the face normal coordinate, or taking a simple/harmonic average, I can at most interpolate, still along the connecting line, in the point nearest to the face center, so that I can minimize the error. But not much more.

@arjun: yes, I probably expressed more obscurely than I tought, but that's what I meant. So you also confirm that there is a third gradient, built via GG and using as face values those reconstructed via limited gradients. This third gradient is NOT itself limited again, but only used in the secondary diffusion term.

@michael: at moment I'm just SDed. Do I still qualify?

Paolo, during the run of the flux computation you have L, R and the node on the face "f". What you need is a subroutine to run before any flow computation that, for each face "f", associates and stores the 4 nodes (two are always L and R) that constitute the tetrahedron. I don't understand why you have only access to L and R...

[FMDenaro]

And note that the same procedure can be used for computing easily the derivative df/dn ...

[arjun]

Originally Posted by FMDenaro

Paolo, during the run of the flux computation you have L, R and the node on the face "f". What you need is a subroutine to run before any flow computation that, for each face "f", associates and stores the 4 nodes (two are always L and R) that constitute the tetrahedron. I don't understand why you have only access to L and R...

I believe CFX does what you are mentioning. They work with shape functions.

[FMDenaro]

Originally Posted by arjun

I believe CFX does what you are mentioning. They work with shape functions.

Well that seems to me a straightforward way to work 😉

[sbaffini]

Originally Posted by arjun

Yes. Not only there but all the other things that need gradient use this gradient. (ie turbulence production etc etc)

I will note down one more thing that many commercial codes when using least square gradients just make a copy of the unlimited gradient and use this (they avoid one more GG application). Then the same gradients are limited and used for reconstruction.

In theory it makes sense but in practice it could be problems in prism type cells so I in wildkatze use extra step of GG and use that as mentioned.

At the moment I think I can't really go with another round of GG. I already have three of them to compute, iteratively, the unlimited gradient. But good to know what I should actually do.

[sbaffini]

Yeah, actually I could build such shape functions. What I meant is that I would like to stick just to information available in cells L and R. Otherwise, storing 4 more scalars per face would blow up my memory. There is a tradeoff with respect to gradients but, honestly, I also don't want to compute them for properties.

[FMDenaro]

Originally Posted by sbaffini

Yeah, actually I could build such shape functions. What I meant is that I would like to stick just to information available in cells L and R. Otherwise, storing 4 more scalars per face would blow up my memory. There is a tradeoff with respect to gradients but, honestly, I also don't want to compute them for properties.

You need to store just 2 nodes, not 4. And the shape functions can be computed out of the flow computation

[arjun]

I made some studies about this diffusion term (Poisson Problem) . I got some improvements over the second order method without increasing the cost too much.

For the test problem 1x1x1 m^3 problem with function sin(PI x) . sin( PI y) . sin(PI z)
the following are the results on tetrahydral meshes.

Mesh Size ___ starccm ___ NewAlog
36051 ___ 6.42E-06 ___ 4.00E-06
181461 ___ 1.17E-06 ___ 7.00E-07
455037 ___ 5.71E-07 ___ 4.55E-07
934437 ___ 2.00E-07 ___ 1.40E-07

The results are compared against starccm, the table shows L2 error norm against exact value for computed solution.

The good thing about it is that on uniform meshes the solution shall be much more accurate than 2nd order method, on skewed meshes the worst case is same as second order method.

It is just primary calculation so take it with a grain of salt. So far this is what i got.

[sbaffini]

You refer to the shape function approach you were mentioning in previous posts, right?

[arjun]

Originally Posted by sbaffini

You refer to the shape function approach you were mentioning in previous posts, right?

No.

Basically the idea is to adapt 4rth order finite difference method to unstructured grid. This is why on uniform meshes the results shall be exactly as the 4rth order method will produce.

The locally refined Cartesian mesh type solvers will benefit the most because majority of those meshes are "boxes" where it shal be close to 4rth order.

I will make more tests when i get time and will outline what exactly is done. Nothing complicated is done and your second order descretization is used as it is.