View Issue Details

IDProjectCategoryView StatusLast Update
0003660OpenFOAMBugpublic2021-04-12 14:58
ReporterjanGaertnerAssigned Tohenry 
PrioritynormalSeverityminorReproducibilityalways
Status closedResolutionno change required 
PlatformUnixOSUbuntuOS Version18.04
Product Version 
Fixed in Version 
Summary0003660: Processor Boundary for Compressibility psi
DescriptionTransonic solvers use the compressibility and pressure to replace the density in the pressure equation and to derive the typical hyperbolic equation system. For this the flux phid is calculated by interpolating the compressibility psi from the EoS to the cell faces and then multiplying this with the volume or mass flux phi, see lines 54-58 in pEqn.H of rhoPimpleFoam OF v8.
The compressibility is calculated when thermo.correct() is called which then sets the cell internal field and boundary values according to the temperature and pressure, see hePsiThermo::calculate().
The problem is now, that when the compressibility psi is interpolated to the cell faces the neighbor cells are required for the standard linear interpolation that is used when fvc::interpolate() is called. At the processor boundaries this is the patchNeighbourField. However, this field is only updated once the boundary conditions are updated. This does not happen when the field entries are manually set, as it is the case for hePsiThermo::calculate().
In conlcusion, this leads to a mass flux error at the processor boundaries, as the compressibility values of the patchNeighbourField are not updated when the flux is solved!

This can be seen when a decomposed case is simulated and the continuity error is visualized for each cell.
TagsNo tags attached.

Activities

henry

2021-04-12 10:34

manager   ~0011978

> At the processor boundaries this is the patchNeighbourField. However, this field is only updated once the boundary conditions are updated.


template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::processorFvPatchField<Type>::patchNeighbourField() const
{
    if (debug && !this->ready())
    {
        FatalErrorInFunction
            << "On patch " << procPatch_.name()
            << " outstanding request."
            << abort(FatalError);
    }
    return *this;
}

i.e. the patchNeighbourField are the values on the patch which are calculated in hePsiThermo::calculate() from the corresponding patchNeighbourField values of the state variables. It is not clear why psi should be inconsistent on the boundaries at this point, could you elaborate?

henry

2021-04-12 10:46

manager   ~0011979

How large is the continuity error you see? For the tests I have run it is tiny and consistent with the drop in accuracy at processor patches due to the reduced effectiveness of the solver preconditioning.

janGaertner

2021-04-12 10:49

reporter   ~0011980

Hello Henry,

the patchNeighbourField can be thought of as ghost cells local to your current processor that are updated with the values from the bordering processor through MPI once the overloaded fvPatchField::updateCoeffs() is called which in return calls the processorFvPatchField::evaluate() function. Only when this function is called the patchNeighbouringField is updated!
It is correct, that the boundary values are updated in hePsiThermo::calculate() but the patchNeighbourField is not. Consider that for an update of the field an MPI call is required.

Typically the update function of the boundary conditions is called before a matrix system is solved. This is built in the fvMatrix functions, see the solve command etc. (It is not easy to see but it is there ;) )
However, for the case of the transonic option of rhoPimpleFoam the fluxes are explicitly calculated before they are added to the matrix system. Hence, the patchNeighbourField of psi is never updated.

henry

2021-04-12 10:49

manager   ~0011981

I added

    this->psi_.correctBoundaryConditions();

to Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()

and I get the same result as expected. Can you provide more details about the problem you see, how to reproduce it and what you propose to fix it as currently I am unable to reproduce it.

janGaertner

2021-04-12 10:51

reporter   ~0011982

To see the error you would need to have large gradients in the temperature field, such that psi is changing significantly over the field. If the temperature changes are small the error is small as well as the changes of psi, compared to the initial solution, are not large either.

henry

2021-04-12 10:53

manager   ~0011983

The small difference in continuity error at processor boundaries I see can be reduced by tightening the pressure solver tolerance, as expected.

janGaertner

2021-04-12 10:53

reporter   ~0011984

I will try to create a small test case where this can be seen. My current simulation is too large and has some either features that I cannot show here.

janGaertner

2021-04-12 10:54

reporter   ~0011985

Maybe to add, the error is not large in most cases. I do not expect different results or current results to be invalidated. However, it is an inconsistency that can be easily solved by adding the update function to the hePsiThermo::calculate() function.

henry

2021-04-12 11:01

manager   ~0011986

I don't see any inconsistency in the code and it is not clear what "update function" should be added to hePsiThermo::calculate(), please be more specific. I added this->psi_.correctBoundaryConditions();
to Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate() and see absolutely no difference in the results, as expected.

janGaertner

2021-04-12 11:12

reporter   ~0011987

Could you tell me why you expect no change?
Possibly I do not see where this happens in the code, but to reiterate the flux phid is calculated before the matrix and the interpolate function uses the patchNeighbourField at the processor boundaries. This patchNeighbourField however is not updated, at least I do not see where this would happen.
I actually did not check, but I can write out the patchNeighbourField as a list over the time steps and see if the values change. However, to see a change I would need to create some kind of temperature gradient in the field that develops over time, thus that there is a time change of psi over the processor boundary.

Could you also tell me what test case you are using? So that we can compare?

henry

2021-04-12 11:13

manager   ~0011988

Last edited: 2021-04-12 11:15

View 2 revisions

> At the processor boundaries this is the patchNeighbourField. However, this field is only updated once the boundary conditions are updated.

template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::processorFvPatchField<Type>::patchNeighbourField() const
{
    if (debug && !this->ready())
    {
        FatalErrorInFunction
            << "On patch " << procPatch_.name()
            << " outstanding request."
            << abort(FatalError);
    }
    return *this;
}

i.e. the patchNeighbourField are the values on the patch which are calculated in hePsiThermo::calculate() from the corresponding patchNeighbourField values of the state variables. It is not clear why psi should be inconsistent on the boundaries at this point.

henry

2021-04-12 11:14

manager   ~0011989

Could you let me know specifically what change to made to the hePsiThermo::calculate() and how it affected your results?

janGaertner

2021-04-12 13:56

reporter   ~0011990

I had added this->psi_.correctBoundaryConditions() to the calculations file.
But I now see that I had misunderstood the processorFvPatchField. I had the understanding, that there are two fields, the normal field at the patch and the patchNeighbourField. However, now I see that they are both the same field and that the update in hePsiThermo does then consider the newest temperature and pressure field.

The error I saw at the boundaries of processors I see now was not caused by psi, it actually has something to do with the interpolation scheme I chose... As I now tried to built a small test case I found that two things changed at once, where one was the change of the psi boundary and the other the interpolation.

To summarize, it was my mistake, as psi is updated with the corrected p and T fields the psi at the boundary is also correct.

Thank you for still taking the time to discuss the issue.

henry

2021-04-12 14:58

manager   ~0011991

User error

Issue History

Date Modified Username Field Change
2021-04-12 09:39 janGaertner New Issue
2021-04-12 10:34 henry Note Added: 0011978
2021-04-12 10:46 henry Note Added: 0011979
2021-04-12 10:49 janGaertner Note Added: 0011980
2021-04-12 10:49 henry Note Added: 0011981
2021-04-12 10:51 janGaertner Note Added: 0011982
2021-04-12 10:53 henry Note Added: 0011983
2021-04-12 10:53 janGaertner Note Added: 0011984
2021-04-12 10:54 janGaertner Note Added: 0011985
2021-04-12 11:01 henry Note Added: 0011986
2021-04-12 11:12 janGaertner Note Added: 0011987
2021-04-12 11:13 henry Note Added: 0011988
2021-04-12 11:14 henry Note Added: 0011989
2021-04-12 11:15 henry Note Edited: 0011988 View Revisions
2021-04-12 13:56 janGaertner Note Added: 0011990
2021-04-12 14:58 henry Assigned To => henry
2021-04-12 14:58 henry Status new => closed
2021-04-12 14:58 henry Resolution open => no change required
2021-04-12 14:58 henry Note Added: 0011991