View Issue Details

IDProjectCategoryView StatusLast Update
0002561OpenFOAMBugpublic2022-05-18 11:14
Reporterniklas.wikstrom Assigned Towill  
PrioritynormalSeveritycrashReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSOtherOS Version(please specify)
Product Versiondev 
Fixed in Versiondev 
Summary0002561: cyclicACMI makes rhoCentral(DyM)Foam crash due to surface interpolation issue
DescriptionrhoCentralFoam crashes (sigFpe) immediately in time loop on a cases containing cyclicACMI patches. The cause of the sigFpe seems to be that "interpolate(rho,neg)" produces a zero rho_neg field, which is then used in the denominator to calculate U_neg (around line 88, rhoCentralFoam.C).

Thank you.
Steps To ReproduceUnpack attached small case.
run included makeMesh.sh script
run rhoCentralDyMFoam (or rhoCentralFoam)
Additional InformationCase runs as expected in OpenFOAM 3.x.
Case runs in OpenFOAM 4.x/dev using "linear" interpolation.
Case runs in OpenFOAM 4.x/dev using e.g. sonicFoam

(For completeness: Tested on Fedora-25;gcc-6.3.1 and Centos-7.3;gcc-4.8.5)
TagsACMI, cyclicACMI

Activities

niklas.wikstrom

2017-05-24 09:52

reporter  

niklas.wikstrom

2017-05-24 10:04

reporter   ~0008185

I should perhaps mention that the same issue appears in version 4.x

MattijsJ

2017-06-23 16:30

reporter   ~0008238

The issue is in cyclicACMIFvPatchField<Type>::patchNeighbourField().

Uncoupled faces assume a patchNeighbour value of 0. Attached version assumes a patchNeighbour value of the blockage patch (for uncoupled faces). Needs a bit of testing.

MattijsJ

2017-06-23 16:30

reporter  

cyclicACMIFvPatchField.C (8,252 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "cyclicACMIFvPatchField.H"
#include "transformField.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF
)
:
    cyclicACMILduInterfaceField(),
    coupledFvPatchField<Type>(p, iF),
    cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{}


template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const dictionary& dict
)
:
    cyclicACMILduInterfaceField(),
    coupledFvPatchField<Type>(p, iF, dict, dict.found("value")),
    cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{
    if (!isA<cyclicACMIFvPatch>(p))
    {
        FatalIOErrorInFunction
        (
            dict
        )   << "    patch type '" << p.type()
            << "' not constraint type '" << typeName << "'"
            << "\n    for patch " << p.name()
            << " of field " << this->internalField().name()
            << " in file " << this->internalField().objectPath()
            << exit(FatalIOError);
    }

    if (!dict.found("value") && this->coupled())
    {
        this->evaluate(Pstream::commsTypes::blocking);
    }
}


template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
    const cyclicACMIFvPatchField<Type>& ptf,
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    cyclicACMILduInterfaceField(),
    coupledFvPatchField<Type>(ptf, p, iF, mapper),
    cyclicACMIPatch_(refCast<const cyclicACMIFvPatch>(p))
{
    if (!isA<cyclicACMIFvPatch>(this->patch()))
    {
        FatalErrorInFunction
            << "' not constraint type '" << typeName << "'"
            << "\n    for patch " << p.name()
            << " of field " << this->internalField().name()
            << " in file " << this->internalField().objectPath()
            << exit(FatalIOError);
    }
}



template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
    const cyclicACMIFvPatchField<Type>& ptf
)
:
    cyclicACMILduInterfaceField(),
    coupledFvPatchField<Type>(ptf),
    cyclicACMIPatch_(ptf.cyclicACMIPatch_)
{}


template<class Type>
Foam::cyclicACMIFvPatchField<Type>::cyclicACMIFvPatchField
(
    const cyclicACMIFvPatchField<Type>& ptf,
    const DimensionedField<Type, volMesh>& iF
)
:
    cyclicACMILduInterfaceField(),
    coupledFvPatchField<Type>(ptf, iF),
    cyclicACMIPatch_(ptf.cyclicACMIPatch_)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Type>
bool Foam::cyclicACMIFvPatchField<Type>::coupled() const
{
    return cyclicACMIPatch_.coupled();
}


template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::cyclicACMIFvPatchField<Type>::patchNeighbourField() const
{
    const Field<Type>& iField = this->primitiveField();
    const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();
    tmp<Field<Type>> tpnf
    (
        cyclicACMIPatch_.interpolate
        (
            Field<Type>
            (
                iField,
                cpp.neighbPatch().faceCells()
            )
        )
    );

    if (doTransform())
    {
        tpnf.ref() = transform(forwardT(), tpnf());
    }

    return cpp.mask()*tpnf + (1.0-cpp.mask())*nonOverlapPatchField();
    //return tpnf;
}


template<class Type>
const Foam::cyclicACMIFvPatchField<Type>&
Foam::cyclicACMIFvPatchField<Type>::neighbourPatchField() const
{
    const GeometricField<Type, fvPatchField, volMesh>& fld =
        static_cast<const GeometricField<Type, fvPatchField, volMesh>&>
        (
            this->primitiveField()
        );

    return refCast<const cyclicACMIFvPatchField<Type>>
    (
        fld.boundaryField()[cyclicACMIPatch_.neighbPatchID()]
    );
}


template<class Type>
const Foam::fvPatchField<Type>&
Foam::cyclicACMIFvPatchField<Type>::nonOverlapPatchField() const
{
    const GeometricField<Type, fvPatchField, volMesh>& fld =
        static_cast<const GeometricField<Type, fvPatchField, volMesh>&>
        (
            this->primitiveField()
        );

    return fld.boundaryField()[cyclicACMIPatch_.nonOverlapPatchID()];
}


template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
(
    scalarField& result,
    const scalarField& psiInternal,
    const scalarField& coeffs,
    const direction cmpt,
    const Pstream::commsTypes
) const
{
    const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();

    // note: only applying coupled contribution

    const labelUList& nbrFaceCellsCoupled =
        cpp.neighbPatch().faceCells();

    scalarField pnf(psiInternal, nbrFaceCellsCoupled);

    // Transform according to the transformation tensors
    transformCoupleField(pnf, cmpt);

    const labelUList& faceCells = cyclicACMIPatch_.faceCells();

    pnf = cyclicACMIPatch_.interpolate(pnf);

    forAll(faceCells, elemI)
    {
        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
    }
}


template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
(
    Field<Type>& result,
    const Field<Type>& psiInternal,
    const scalarField& coeffs,
    const Pstream::commsTypes
) const
{
    const cyclicACMIPolyPatch& cpp = cyclicACMIPatch_.cyclicACMIPatch();

    // note: only applying coupled contribution

    const labelUList& nbrFaceCellsCoupled = cpp.neighbPatch().faceCells();

    Field<Type> pnf(psiInternal, nbrFaceCellsCoupled);

    // Transform according to the transformation tensors
    transformCoupleField(pnf);

    const labelUList& faceCells = cyclicACMIPatch_.faceCells();

    pnf = cyclicACMIPatch_.interpolate(pnf);

    forAll(faceCells, elemI)
    {
        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
    }
}


template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
(
    fvMatrix<Type>& matrix
)
{
    const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();

    // nothing to be done by the AMI, but re-direct to non-overlap patch
    // with non-overlap patch weights
    const fvPatchField<Type>& npf = nonOverlapPatchField();

    const_cast<fvPatchField<Type>&>(npf).manipulateMatrix(matrix, 1.0 - mask);
}


template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateCoeffs()
{
    // Update non-overlap patch - some will implement updateCoeffs, and
    // others will implement evaluate

    // Pass in (1 - mask) to give non-overlap patch the chance to do
    // manipulation of non-face based data

    const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
    const fvPatchField<Type>& npf = nonOverlapPatchField();
    const_cast<fvPatchField<Type>&>(npf).updateWeightedCoeffs(1.0 - mask);
}


template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::write(Ostream& os) const
{
    fvPatchField<Type>::write(os);
    this->writeEntry("value", os);
}


// ************************************************************************* //
cyclicACMIFvPatchField.C (8,252 bytes)   

MattijsJ

2017-06-30 10:42

reporter   ~0008294

Above fix adjusts only the patchNeighbour value so will unfortunately make e.g. Gauss gradient incorrect on the partially overlapping bits since the coupled part should only include coupled contributions and the blockage part only the 'wall' contributions. We probably have to have the value field calculated separately, instead of being calculated from internal and neighbour field.

will

2022-05-18 11:14

manager   ~0012593

Numerical issues with ACMI-based patch couplings have been resolved by the new Non-Conformal Coupled (NCC) development. See the following commit for explanation, and instructions for how to use NCC and convert cases from AMI.

https://github.com/OpenFOAM/OpenFOAM-dev/commit/569fa31d09f98e29d1aaf84d40bb16043f104ec6

Issue History

Date Modified Username Field Change
2017-05-24 09:52 niklas.wikstrom New Issue
2017-05-24 09:52 niklas.wikstrom File Added: slidingBlockTestCase.tgz
2017-05-24 09:52 niklas.wikstrom Tag Attached: ACMI
2017-05-24 09:52 niklas.wikstrom Tag Attached: cyclicACMI
2017-05-24 10:04 niklas.wikstrom Note Added: 0008185
2017-06-23 16:30 MattijsJ Note Added: 0008238
2017-06-23 16:30 MattijsJ File Added: cyclicACMIFvPatchField.C
2017-06-30 10:42 MattijsJ Note Added: 0008294
2022-05-18 11:14 will Assigned To => will
2022-05-18 11:14 will Status new => resolved
2022-05-18 11:14 will Resolution open => fixed
2022-05-18 11:14 will Fixed in Version => dev
2022-05-18 11:14 will Note Added: 0012593