Small offset when reading gradients (fluid-fluid coupling specific)

[MakisH]

As presented in the 14th OpenFOAM Conference, when reading gradients we end up with a minor increase in the respective quantity.

We observed this in fluid-fluid coupling (see #60) when coupling {pressure} with {velocity, pressure gradient} and using a coarse mesh. However, this problem also appears in CHT cases (reading heat flux), but we hadn't noticed until now since we are always using very fine meshes there.

In the following figure you can see this error when coupling gradients. When coupling only velocity and pressure, this problem does not appear, although we still get an oscillation around the interface.

We currently suspect a technical problem in the way we read and write gradients, or some other numerical issue.

openfoam-adapter/FF/PressureGradient.C

Lines 21 to 75 in 857f18c

	void preciceAdapter::FF::PressureGradient::write(double * buffer)
	{
	int bufferIndex = 0;
	// For every boundary patch of the interface
	for (uint j = 0; j < patchIDs_.size(); j++)
	{
	int patchID = patchIDs_.at(j);
	// Get the pressure gradient boundary patch
	scalarField gradientPatch
	=
	refCast<fixedValueFvPatchScalarField>
	(
	p_->boundaryFieldRef()[patchID]
	).snGrad();
	// For every cell of the patch
	forAll(gradientPatch, i)
	{
	// Copy the pressure gradient into the buffer
	buffer[bufferIndex++]
	=
	-gradientPatch[i];
	}
	}
	}
	void preciceAdapter::FF::PressureGradient::read(double * buffer)
	{
	int bufferIndex = 0;
	// For every boundary patch of the interface
	for (uint j = 0; j < patchIDs_.size(); j++)
	{
	int patchID = patchIDs_.at(j);
	// Get the pressure gradient boundary patch
	scalarField & gradientPatch
	=
	refCast<fixedGradientFvPatchScalarField>
	(
	p_->boundaryFieldRef()[patchID]
	).gradient();
	// For every cell of the patch
	forAll(gradientPatch, i)
	{
	// Set the pressure gradient as the buffer value
	gradientPatch[i]
	=
	buffer[bufferIndex++];
	}
	}
	}

You can try this in our pipe-pipe tutorial.

[uekerman]

We should test whether we can reproduce this behavior with a simple heat transfer case, so coupling laplacianFoam to laplacianFoam; maybe even steady-state.

[davidscn]

I set up a heat transfer case coupling deal.II and a custom laplacianFoam considering non-zero RHS in order to investigate this. The example is our partitioned-heat example, which has a quadratic solution in space and a linear solution in time. OpenFOAM uses 5 x 5 FV elements and deal.II discretizes the domain using quadratic shape functions which solves the problem more or less exact.

Doing a single explicitly coupled time step in time reveals the following curvatures:

and the deal.II side reports an L2 error of e = 0.00604907. OpenFOAM acts here as Neumann participant (the right part of the domain) reading Heat-Flux (the Temperature gradients) and Temperature values.

There is clearly a kink at the coupling interface. However, there is no permanent offset of the Temperature curvature. The reason is probably the fixedGradient boundary condition, which computes respective face values using linear reconstruction in space. (In the picture above, it might be an additional representation issue with paraView, but I double checked the quantitative values in OpenFOAM itself.)

When writing gradients, we use surface normal gradients. I'm not yet completely sure because I have not fully investigated this case, but the problem might be similar. However, the surface normal gradient scheme can be specified in the fvSchemes dict as opposed to the fixed gradient and there might be an option I have not yet seen.

In order to investigate this further I boiled our partitioned heat example from second order in space to first order in space down (RHS = beta) and the example looks as follows:

and the deal.II side reports an L2 error of e = 9.82347e-14. As we can see, we don't have any kinks present, i.e. a linear error on the OpenFOAM side at the interface can be confirmed.

In addition, the observation would be consistent with the issue occurring only for coarse meshes as the linear error goes down with very meshes in interface normal direction.

Regarding the Fluid-Fluid coupling above:
I took a closer look at out partitioned pipe example. I think the general issue is different: taking the setup and running it in an explicit coupling scheme is not even stable. I think there is an issue on how and which boundary conditions are applied. The pressure across the interface goes just crazy during the simulation

Tweaking the boundary conditions a bit and running an explicit coupling gives me

Not sure if this is still some linear error; the overall curvature is more or less linear in the region so that I would not expect an issue at this point.

[MakisH]

Thank you very much for looking so carefully on this! Just a question:

In order to investigate this further I boiled our partitioned heat example from second order in space to first order in space down (RHS = beta) and the example looks as follows:

What exactly did you change here? I guess something in the fvSchemes? Anything else?

OpenFOAM acts here as Neumann participant (the right part of the domain) reading Heat-Flux (the Temperature gradients) and Temperature values.

I assume here you mean that deal.II reads Temperature values, right?

[davidscn]

What exactly did you change here? I guess something in the fvSchemes? Anything else?

No it is indeed a change within the OpenFOAM solver. @uekerman thought the partitioned-heat solver with OpenFOAM might be valuable for different reasons. We still need to discuss where to put it..

I assume here you mean that deal.II reads Temperature values, right?

..reading Heat-flux and writing Temperature values. So, yes deal.II reads them.

Will look a bit more in the partitioned pipe as there are some things I want to try.

[uekerman]

What if OpenFOAM is the Dirichlet participant. Do we then also see a kink?

[davidscn]

After more investigations (see #189 #188 #187 ) I can now finally run the partitioned-heat OpenFOAM-OpenFOAM case and visualize the correct results:

[davidscn]

To sum this up: I cannot reproduce the original issue here with a partitioned-heat example and I still think that you cannot simply use fixedValue (not even an inlet) and fixedGradient boundary conditions in order to tackle the FF coupling. There are more specialized BCs like InletOutlet or fixedFluxPressure in OpenFOAM which might be more appropriate for this type of coupling.

[davidscn]

@MakisH what I observed for the CHT cases now (considering the fixes) looks as follows:

I only observe this for intermediate time steps, not for the final time step. For me, it looks like a thermal boundary layer/ the transient behavior of the simulation and I can also observe it using the default (fine) mesh and coarser meshes. Is this what you observed or was it something different?

[MakisH]

Yes, this is what I observed as well. Do you mean that this has a physical explanation?

[uekerman]

I can also observe it using the default (fine) mesh and coarser meshes

And the quantitative result is then rather independent of the mesh?

[davidscn]

After some further investigations I'm a bit stuck here. Trying to visualize the results with paraFoam or vtk exports results always in a small 'gap' between the participants. Grouping the participants in paraView together and plotting them afterwards results in a segfault since both participants contain different data fields which seems to cause trouble for paraView:

The green curve belongs to the fluid part and the curve does not fully reach the interface located at y=0.25. However, I am a bit puzzled how to bypass the problem.

[MakisH]

To sum this up: I cannot reproduce the original issue here with a partitioned-heat example and I still think that you cannot simply use fixedValue (not even an inlet) and fixedGradient boundary conditions in order to tackle the FF coupling. There are more specialized BCs like InletOutlet or fixedFluxPressure in OpenFOAM which might be more appropriate for this type of coupling.

I also noticed that the previous coupling of ATHLET-OpenFOAM (by GRS) used modified boundary conditions based on inletOutlet (for the temperature), normalInletOutletVelocity (velocity), fixedMean (velocity), fixedValue/buoyantPressure (pressure). See presentation at the OpenFOAM Workshop 2014, slide 8.

After some further investigations I'm a bit stuck here. Trying to visualize the results with paraFoam or vtk exports results always in a small 'gap' between the participants. Grouping the participants in paraView together and plotting them afterwards results in a segfault since both participants contain different data fields which seems to cause trouble for paraView:

I guess the gap corresponds to the distance between the interface cell center and the corresponding face?

Why "different data fields"? Both are named T, both are (I guess) cell data, I do not understand what the incompatibility is that could lead to a segfault. What does the message say, exactly? (although this is often not helpful)

The green curve belongs to the fluid part and the curve does not fully reach the interface located at y=0.25. However, I am a bit puzzled how to bypass the problem.

Maybe some conversion between point and cell data would help?

[davidscn]

I also noticed that the previous coupling of ATHLET-OpenFOAM (by GRS) used modified boundary conditions based on inletOutlet (for the temperature), normalInletOutletVelocity (velocity), fixedMean (velocity), fixedValue/buoyantPressure (pressure). See presentation at the OpenFOAM Workshop 2014, slide 8.

What do you mean by previous? There was a Fluid-Fluid coupling before?

Why "different data fields"? Both are named T, both are (I guess) cell data, I do not understand what the incompatibility is that could lead to a segfault. What does the message say, exactly? (although this is often not helpful)

With different I mean quantities such as U, which are apparently not part of the solid. The output is some lengthy and cryptic paraView stack trace.

Maybe some conversion between point and cell data would help?

The problem might be that the stored gradient is not resolved correctly? How can I exclude that a potential misalignment does not stem from paraView. Visualizing these coupled cases is a bit of a challenge..

[MakisH]

What do you mean by previous? There was a Fluid-Fluid coupling before?

Only between ATHLET and OpenFOAM, using an in-house system, not preCICE. Reproducing something similar with preCICE is the project I am working on.

With different I mean quantities such as U, which are apparently not part of the solid. The output is some lengthy and cryptic paraView stack trace.

In ParaView, you can disable specific fields. Maybe this helps?

The problem might be that the stored gradient is not resolved correctly? How can I exclude that a potential misalignment does not stem from paraView. Visualizing these coupled cases is a bit of a challenge..

It definitely is!

One approach I can think of would be inspecting specific points in a spreadsheet: https://docs.paraview.org/en/latest/UsersGuide/selectingData.html

Another approach would be converting first to VTK and comparing: do you get the same effect?

[davidscn]

Only between ATHLET and OpenFOAM, using an in-house system, not preCICE. Reproducing something similar with preCICE is the project I am working on.

Ah ok, but all Fluid-Fluid OpenFoam progress is collected in the tutorials, right?

In ParaView, you can disable specific fields. Maybe this helps?

This helps, indeed, and results in no further offsets. The two lines are essentially connected.

Another approach would be converting first to VTK and comparing: do you get the same effect?

foamToVTK has the same issue as before, doesn't really help and confuses a bit more since the implementations among the vendors seems to be different.

[MakisH]

Ah ok, but all Fluid-Fluid OpenFoam progress is collected in the tutorials, right?

yes

This helps, indeed, and results in no further offsets. The two lines are essentially connected.

Could you please still upload a picture here for documentation purposes?

[davidscn]

Generated with paraFoam by loading the temperature field exclusively.

[MakisH]

To sum this up: I cannot reproduce the original issue here with a partitioned-heat example and I still think that you cannot simply use fixedValue (not even an inlet) and fixedGradient boundary conditions in order to tackle the FF coupling. There are more specialized BCs like InletOutlet or fixedFluxPressure in OpenFOAM which might be more appropriate for this type of coupling.

I also noticed that the previous coupling of ATHLET-OpenFOAM (by GRS) used modified boundary conditions based on inletOutlet (for the temperature), normalInletOutletVelocity (velocity), fixedMean (velocity), fixedValue/buoyantPressure (pressure). See presentation at the OpenFOAM Workshop 2014, slide 8.

Just to document:

technically, InletOutlet seems to work with the adapter for pressure and velocity. For the gradients, however, we need to change the class to something that can work with both InletOutletand fixedGradient, otherwise it crashes e.g. in PressureGradient.C. As a next step, I would look into the implementation of InletOutlet and try to understand how it is implemented type-wise.

So, currently the fluid-fluid module technically works for fixedValue and fixedGradient, but these may be terrible selections in the first place and may explain why only these FF cases show these offsets.

[MakisH]

To document the state: we are working on this together with @thesamriel, in the context of his Master's thesis.

To sum this up: I cannot reproduce the original issue here with a partitioned-heat example and I still think that you cannot simply use fixedValue (not even an inlet) and fixedGradient boundary conditions in order to tackle the FF coupling. There are more specialized BCs like InletOutlet or fixedFluxPressure in OpenFOAM which might be more appropriate for this type of coupling.

Indeed, what @thesamriel found already, is that fixedFluxExtrapolatedPressure boundary conditions on the side that reads velocity and writes pressure (without exchanging pressure gradient) is enough to solve this issue. More details once we integrate the fix into the tutorials.

[raynoldtan]

Dear all,

I am trying to do a fluid fluid coupling involving an interFoam solver on the left and right domain respectively( See attached picture). The case presented here is the wave flume case inside the tutorial case of openFoam interFoam/laminar/wave folder.

In this case, inside precise dict, the left side is writing velocity and alpha while reading in pressure and purge. the right side is writing pressure and prgh and reading in velocity and alpha. Alpha is the interphase fraction. The issue I am facing is some sort of wave reflection on the left domain, more specifically, it is an issue of back flow. In this example, I have simply used here fixed gradient for the velocity outlet on the left domain and fixed value for the velocity inlet of the right domain. The simulation was able to run for a while until wave reflections start building up on the left domain outlet. Please advise on what type of boundary conditions can be used on the left domain outlet to prevent this issue? I have also tried inletOutlet on the left domain outlet but it only allow the simulation to run longer, but wave reflection eventually build up on the left domain outlet.

![Screenshot 2022-08-16 at 10 41 43](https://user-image

s.githubuse

rcontent.com/107682091/184850858-07e12aee-0781-4614-9362-85e175a3871b.png)

[raynoldtan]

I am attaching the pictures here again. The first two are the velocity field. The third is the alpha field and the last one is the prgh field

Any advice is much appreciated !

.

[MakisH]

Dear @raynoldtan,

even though this is also related to fluid-fluid coupling, there can be several reasons while this is happening, and it does not seem to be directly related to the specific gradients reading issue. Please describe your case/support request thoroughly in the preCICE forum. There are other people in our community that have coupled interFoam (at least for FSI) and may also be able to help.

[raynoldtan]

@MakisH

HI Makis, we have identified the issue as something related to the outlet boundary of the left domain. Subsequently, we were able to get the simulation running and obtain better results through the use of a wave absorption Boundary condition from olaflow (openfoam solver catered for ocean engineering problems). Although there is still some wave reflection at the left domain outlet even after using the wave Absorption boundary. This of course could not be completely removed unless a complete wave absorbing boundary condition can be devised. The complexity in this issues lies in the fact that this is a multiphase flow.

[MakisH]

Closed by #281
The fix was mainly the fixedFluxExtrapolatedPressure boundary condition.