[naffrancois]

Dear all,

I am trying to find out a way to interpolate a field between two grids, where one grid has its nodes moved and/or more nodes than the other one. Is there a conservative and possibly simple way to achieve this ? I can think of a sort of nearest neighbour approach + piecewise/linear interpolation but this cannot be conservative.

How is it done in the frame of unsteady mesh adaptation ?

Thank you for your help

[FMDenaro]

Originally Posted by naffrancois

Dear all,

I am trying to find out a way to interpolate a field between two grids, where one grid has its nodes moved and/or more nodes than the other one. Is there a conservative and possibly simple way to achieve this ? I can think of a sort of nearest neighbour approach + piecewise/linear interpolation but this cannot be conservative.

How is it done in the frame of unsteady mesh adaptation ?

Thank you for your help

I am aware of some proposals in literqature but I personally never worked on conservative interpolation.
Maybe you can think to intepolate the flux functions at the new edges? This way you could compute the sum of the fluxes in the new FV.

[aerosayan]

Would following literature on overset/chimera mesh based solvers help?


They're conservative in nature and interpolate data between grids. Interpolation error is a function of cell size differences between the two grids.

[naffrancois]

Originally Posted by FMDenaro

I am aware of some proposals in literqature but I personally never worked on conservative interpolation.
Maybe you can think to intepolate the flux functions at the new edges? This way you could compute the sum of the fluxes in the new FV.

Thank you FMDenaro pour your answer. Unfortunately I am not sure to understand your suggestion but this is maybe because I did not explain well my doubts.

Let's imagine I have two domains, two squares of same dimensions to keep things "simple", meshed with triangles. The two triangulations do not share anything. I have a solution field on mesh 1 and want to interpolate it on mesh 2 keeping total mass and energy constant (having Euler equations in mind). I am wondering if this is even possible, at least much more complicated than I hoped. I remember a paper where I think they did something similar, but I cannot find it anymore. They considered all possible element intersections, not so easy though...

[FMDenaro]

Originally Posted by naffrancois

Thank you FMDenaro pour your answer. Unfortunately I am not sure to understand your suggestion but this is maybe because I did not explain well my doubts.

Let's imagine I have two domains, two squares of same dimensions to keep things "simple", meshed with triangles. The two triangulations do not share anything. I have a solution field on mesh 1 and want to interpolate it on mesh 2 keeping total mass and energy constant (having Euler equations in mind). I am wondering if this is even possible, at least much more complicated than I hoped. I remember a paper where I think they did something similar, but I cannot find it anymore. They considered all possible element intersections, not so easy though...

I mean that instead of interpolating the solution maybe one can interpolate the flux functions and then compute the solution.
But I think you can find the literature on conservative interpolation on google.

[naffrancois]

Originally Posted by aerosayan

Would following literature on overset/chimera mesh based solvers help?


They're conservative in nature and interpolate data between grids. Interpolation error is a function of cell size differences between the two grids.

It may, I did not think about that, thanks for the idea.

[sbaffini]

First of all, a formal conservation statement is needed. Assuming a finite volume discretization at second order, with constant cell values, the conservation statement for a scalar variable is something like: volume integral of new interpolated variables over old grid cells must be equal.

Unfortunately, there is no way to achieve conservation by relying on interpolation. Chimera and overset are not, indeed, conservative when done by interpolation alone.

The only way I know is to actually split the new grid with the faces from the old one and compute the new variables by proper volume average.

If the two meshes cover the exact same domain it is still complex for general grids, but doable. If you allow the domain to change I have doubts it is doable without further conditions.

[naffrancois]

Originally Posted by sbaffini

The only way I know is to actually split the new grid with the faces from the old one and compute the new variables by proper volume average.

I see... that's rather complicated indeed for general tets grids. I will do a literature review anyway with the ideas you all gave me. Although I believe now there's no simple magic solution here

[aerosayan]

Originally Posted by sbaffini

First of all, a formal conservation statement is needed. Assuming a finite volume discretization at second order, with constant cell values, the conservation statement for a scalar variable is something like: volume integral of new interpolated variables over old grid cells must be equal.

Unfortunately, there is no way to achieve conservation by relying on interpolation. Chimera and overset are not, indeed, conservative when done by interpolation alone.

The only way I know is to actually split the new grid with the faces from the old one and compute the new variables by proper volume average.

If the two meshes cover the exact same domain it is still complex for general grids, but doable. If you allow the domain to change I have doubts it is doable without further conditions.

Splitting the cells does allow to conserve values, but a somewhat more complex but accurate solution might be to create an empty buffer region to be filled with triangle cells between the two grids' interface. We fill it with advancing front triangulation, which connects cells from the two grids perfectly. Then your solver just treats it as a normal partitioned mesh. This obviously is meant for static grids, and doesn't work very well for moving grids. For moving grids, your buffer region needs to be large, and you can create perfect connectivity between the two grids in each iteration, as the grids move.

[sbaffini]

Originally Posted by aerosayan

Splitting the cells does allow to conserve values, but a somewhat more complex but accurate solution might be to create an empty buffer region to be filled with triangle cells between the two grids' interface. We fill it with advancing front triangulation, which connects cells from the two grids perfectly. Then your solver just treats it as a normal partitioned mesh. This obviously is meant for static grids, and doesn't work very well for moving grids. For moving grids, your buffer region needs to be large, and you can create perfect connectivity between the two grids in each iteration, as the grids move.

Yes, that's more or less how it is done in chimera/overset conservative frameworks (but with care), which for full overlapping grids, I think, just kind of becomes splitting the old mesh with the new one and redo the volume integrals.

I don't exclude that there might be some interpolation frameworks where you add conservation as additional constraint but, in the end, that constraint would just need the same information of a full grid split. That is, interpolation is just the last step, which in the second order framework I used as example is just trivial. Higher order interpolations might require to conserve additional volume moments.

I purposedly left out the treatment of specific variables here and just mentioned scalars. The additional complexity is that, of course, you need to conserve the actually conserved variables, which might not be your independent variables. Then you need to properly switch between them before and after the interpolation.

[sbaffini]

Also, you only implicitly explain your task, by mentioning unsteady moving mesh.

If, instead, you just need to interpolate without local conservation, you could indeed proceed by a simple interpolation and later fix conservation globally, redistributing the imbalance according to, say, a local mesh metric. But I don't have a reference, you are on your own here. It is typically done with immersed boundaries

[naffrancois]

Originally Posted by sbaffini

Also, you only implicitly explain your task, by mentioning unsteady moving mesh.

You are right, it is indeed about mesh adaptation of a 3d unstructured solver for unsteady flows. I am at step 0, trying to evaluate the amount of work involved depending on the different approaches found here and there.

This mesh interpolation question arises with libraries such as MMG (https://www.mmgtools.org/) which is from my understanding a remeshing strategy based on a size map. Apparently it is used within the SU2 framework, I do not know for now though if they use it for unsteady flows or exclusively for steady computations. If only for steady flows, conservation is not an issue as the solution is restarted anyway from scratch (or any crude interpolation) on the adapted mesh. I could check though if they apply it for unsteady flow and if so how they manage solution interpolation.

Otherwise I am also interested in building my own AMR-like adaptation method, in that case I think conservation can be enforced in a more natural way thanks to the relation between the adapted meshes. But this approach for sure needs a big amount of modifications to my code, that's why I was considering first, probably naively, remeshing strategies through existing libraries.

[arjun]

Originally Posted by naffrancois

You are right, it is indeed about mesh adaptation of a 3d unstructured solver for unsteady flows. I am at step 0, trying to evaluate the amount of work involved depending on the different approaches found here and there.

This mesh interpolation question arises with libraries such as MMG (https://www.mmgtools.org/) which is from my understanding a remeshing strategy based on a size map. Apparently it is used within the SU2 framework, I do not know for now though if they use it for unsteady flows or exclusively for steady computations. If only for steady flows, conservation is not an issue as the solution is restarted anyway from scratch (or any crude interpolation) on the adapted mesh. I could check though if they apply it for unsteady flow and if so how they manage solution interpolation.

Otherwise I am also interested in building my own AMR-like adaptation method, in that case I think conservation can be enforced in a more natural way thanks to the relation between the adapted meshes. But this approach for sure needs a big amount of modifications to my code, that's why I was considering first, probably naively, remeshing strategies through existing libraries.

There is a very simple way to handle all this and this is the best way to handle this.


Don't worry about anything in interpolation. After the interpolation is done. Simple solve continuity one time. This enforces divergence free condition on velocity field which is what you need here.


If you can not do this for some reason then second best option is to just skip first momentum equation solve in the time step this way continuity gets the velocity field before momentum.

Both these options look the same but they are not. In second option you might see pressure spikes that no one will be explaining to you why happens. (I know but i will leave it as a puzzle in case you reach that point).


First options is best option if you can do so.

[FMDenaro]

I am not sure but I remember some papers of T. Barth about adaptive unstructured grids with conservation properties. Sorry but I don't remember exactly the sources ...

[naffrancois]

Originally Posted by arjun

There is a very simple way to handle all this and this is the best way to handle this.


Don't worry about anything in interpolation. After the interpolation is done. Simple solve continuity one time. This enforces divergence free condition on velocity field which is what you need here.


If you can not do this for some reason then second best option is to just skip first momentum equation solve in the time step this way continuity gets the velocity field before momentum.

Both these options look the same but they are not. In second option you might see pressure spikes that no one will be explaining to you why happens. (I know but i will leave it as a puzzle in case you reach that point).


First options is best option if you can do so.

Thank you Arjun for your explanation. Unfortunately I am not sure this procedure applies to my cases of interest as I am working primarily with highly compressible flows. Tell me if I understand it wrong. I think the interpolation has to be done carefully to maintain TVD property and hopefully conservation of mass and energy.

By the way, these pressure spikes you mention, aren't they due to spurious oscilations caused by a polynomial interpolation ?

[naffrancois]

Originally Posted by FMDenaro

I am not sure but I remember some papers of T. Barth about adaptive unstructured grids with conservation properties. Sorry but I don't remember exactly the sources ...

Thank you, I will definitely look in his papers

[LuckyTran]

There needs to be a formal statement of conservativeness so some lit review is warranted.

Let's say you want the interpolated solution which satisfies the local governing equations on the new mesh, this is trying to find the inverse to a dual mapping problem.

Given two different size meshes and , and and the respective solutions of the discretized equations on each mesh , you have a transformation matrix which maps the input mesh coordinate vectors to the solution vectors to and to which I will label and which is the thingy you hope to learn.
and

For any given mesh you know the (linear) discretization matrix of the linearized system of equations and with:
because both are approximation solutions to the same global conservation equation. I don't think we need to consider non-linear L's because CFD doesn't work that way.

Therefore you have
and the transformation matrix you're looking for is

where denotes the pseudo-inverses of their respective matrices.

Note that the solution of this system is at least as expensive as solving so you have to determine for some given conservativeness criteria whether it is worth it or not That is, it is well within the realm of expectation that you can do a very fast 0th or linear interpolation and just crank the solution on the new mesh faster than you can find this mapping. Also note that you wouldn't actually solve for the pseudo-inverses, you would solve for A in the same way you solve for V without trying to invert L. For local mesh adaptation, you operate on only subsets of this full matrix system.

[naffrancois]

Originally Posted by LuckyTran

There needs to be a formal statement of conservativeness so some lit review is warranted.

Let's say you want the interpolated solution which satisfies the local governing equations on the new mesh, this is trying to find the inverse to a dual mapping problem.

Given two different size meshes and , and and the respective solutions of the discretized equations on each mesh , you have a transformation matrix which maps the input mesh coordinate vectors to the solution vectors to and to which I will label and which is the thingy you hope to learn.
and

For any given mesh you know the (linear) discretization matrix of the linearized system of equations and with:
because both are approximation solutions to the same global conservation equation. I don't think we need to consider non-linear L's because CFD doesn't work that way.

Therefore you have
and the transformation matrix you're looking for is

where denotes the pseudo-inverses of their respective matrices.

Note that the solution of this system is at least as expensive as solving so you have to determine for some given conservativeness criteria whether it is worth it or not That is, it is well within the realm of expectation that you can do a very fast 0th or linear interpolation and just crank the solution on the new mesh faster than you can find this mapping. Also note that you wouldn't actually solve for the pseudo-inverses, you would solve for A in the same way you solve for V without trying to invert L. For local mesh adaptation, you operate on only subsets of this full matrix system.

Thank you LuckyTran for this different view, it sheds some light on the potential costs involved and trade off necessary to make the adaptation strategy viable. I am getting more convinced that I should opt for a local subdivision approach (something in line with AMR) with a grid hierarchy rather than a full remeshing approach (mesh movement+enrichment).

[naffrancois]

Many thanks to you all for your inputs. In the end I applied the method described in the following paper:

https://doi.org/10.1002/nme.2951

It describes an algorithm to construct a mesh as the intersection between two meshes, e.g. a mesh containing both mesh 1 and mesh 2. From there, maximum principle and global conservation can be achieved in a straightforward way.