OpenFOAM – errorEstimate

In this article, we propose a slightly modified version of an error estimator proposed by the fork foam-extend.

It is a residual error estimator that estimates the discretization error due to the finite volume method for 2D/3D meshes
For more information on error estimators, please refer to the work of H. Jsak or to CODE_ASTER documentation.

This version of the error estimator recoded to be compatible with OpenFoam v5.0 and available here is for information purposes only and must still be validated.

 

  • files
LIB = $(FOAM_LIBBIN)/liberrorEstimation

 

  • option
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

LIB_LIBS = \
    -lfiniteVolume \
    -lmeshTools

 

  • errorEstimate.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::errorEstimate

Description
    Residual error estimation

SourceFiles
    errorEstimate.C

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

#ifndef errorEstimate_H
#define errorEstimate_H

#include "volFields.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                        Class errorEstimate Declaration
\*---------------------------------------------------------------------------*/

template<class Type>
class errorEstimate
:
    public refCount
{
    // Private data

        // Reference to GeometricField<Type, fvPatchField, volMesh>
        const GeometricField<Type, fvPatchField, volMesh>& psi_;

        //- Dimension set
        dimensionSet dimensions_;

        //- Cell residual pointer
        Field<Type> residual_;

        //- Normalisation factor
        scalarField normFactor_;


    // Private Member Functions

        //- Return boundary condition types for the error field
        wordList errorBCTypes() const;

public:

    // Static data members

    ClassName("errorEstimate");


    // Constructors

        //- Construct from components
        errorEstimate
        (
            const GeometricField<Type, fvPatchField, volMesh>& psi,
            const dimensionSet& ds,
            const Field<Type>& res,
            const scalarField& norm
        );

        //- Construct as copy
        errorEstimate(const errorEstimate<Type>&);


    // Destructor

        ~errorEstimate();


    // Member Functions

        // Access

            //- Return field
            const GeometricField<Type, fvPatchField, volMesh>& psi() const
            {
                return psi_;
            }

            //- Return residual dimensions
            const dimensionSet& dimensions() const
            {
                return dimensions_;
            }

            //- Construct and return a clone
       	    tmp<errorEstimate<Type> > clone() const;
        // Raw residual (for calculus)

            Field<Type>& res()
            {
                return residual_;
            }

            const Field<Type>& res() const
            {
                return residual_;
            }


        // Error Estimate

            //- Cell residual (volume intensive)
            tmp<GeometricField<Type, fvPatchField, volMesh> > residual() const;

            //- Normalisation factor
            tmp<volScalarField> normFactor() const;

            //- Error estimate
            tmp<GeometricField<Type, fvPatchField, volMesh> > error() const;


    // Member Operators

        void operator=(const errorEstimate<Type>&);
        void operator=(const tmp<errorEstimate<Type> >&);

        void negate();

        void operator+=(const errorEstimate<Type>&);
        void operator+=(const tmp<errorEstimate<Type> >&);

        void operator-=(const errorEstimate<Type>&);
        void operator-=(const tmp<errorEstimate<Type> >&);

        void operator+=(const GeometricField<Type,fvPatchField,volMesh>&);
        void operator+=(const tmp<GeometricField<Type,fvPatchField,volMesh> >&);

        void operator-=(const GeometricField<Type,fvPatchField,volMesh>&);
        void operator-=(const tmp<GeometricField<Type,fvPatchField,volMesh> >&);

        void operator+=(const dimensioned<Type>&);
        void operator-=(const dimensioned<Type>&);

        void operator*=(const volScalarField&);
        void operator*=(const tmp<volScalarField>&);

        void operator*=(const dimensioned<scalar>&);


    // Friend Functions

    // Friend Operators
};


// * * * * * * * * * * * * * * * Global functions  * * * * * * * * * * * * * //

template<class Type>
void checkMethod
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&,
    const char*
);

template<class Type>
void checkMethod
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&,
    const char*
);

template<class Type>
void checkMethod
(
    const errorEstimate<Type>&,
    const dimensioned<Type>&,
    const char*
);


// * * * * * * * * * * * * * * * Global operators  * * * * * * * * * * * * * //

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >&,
    const dimensioned<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const dimensioned<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >&,
    const dimensioned<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const dimensioned<Type>&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const GeometricField<Type, fvPatchField, volMesh>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>&,
    const dimensioned<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >&,
    const dimensioned<Type>&
);


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField&,
    const tmp<errorEstimate<Type> >&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>&,
    const tmp<errorEstimate<Type> >&
);


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>&,
    const errorEstimate<Type>&
);

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>&,
    const tmp<errorEstimate<Type> >&
);


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "errorEstimate.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

 

  • errorEstimate.C
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

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

#include "errorEstimate.H"
#include "zeroGradientFvPatchField.H"
#include "fixedValueFvPatchField.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //


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

template<class Type>
Foam::wordList Foam::errorEstimate<Type>::errorBCTypes() const
{
    // Make the boundary condition type list
    // Default types get over-ridden anyway
    wordList ebct
    (
        psi_.boundaryField().size(),
        zeroGradientFvPatchField<Type>::typeName
    );

    forAll (psi_.boundaryField(), patchI)
    {
        if (psi_.boundaryField()[patchI].fixesValue())
        {
            ebct[patchI] = fixedValueFvPatchField<Type>::typeName;
        }
    }

    return ebct;
}


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

// Construct from components
template<class Type>
Foam::errorEstimate<Type>::errorEstimate
(
    const GeometricField<Type, fvPatchField, volMesh>& psi,
    const dimensionSet& ds,
    const Field<Type>& res,
    const scalarField& norm
)
:
    psi_(psi),
    dimensions_(ds),
    residual_(res),
    normFactor_(norm)
{}


// Construct as copy
template<class Type>
Foam::errorEstimate<Type>::errorEstimate(const Foam::errorEstimate<Type>& ee)
:
    refCount(),
    psi_(ee.psi_),
    dimensions_(ee.dimensions_),
    residual_(ee.residual_),
    normFactor_(ee.normFactor_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class Type>
Foam::errorEstimate<Type>::~errorEstimate()
{}


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

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::errorEstimate<Type>::residual() const
{
    tmp<GeometricField<Type, fvPatchField, volMesh> > tres
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "residual" + psi_.name(),
                psi_.mesh().time().timeName(),
                psi_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            psi_.dimensions()/dimTime,
            errorBCTypes()
        )
    );

    GeometricField<Type, fvPatchField, volMesh>& res = tres();

    res.internalField() = residual_;
    res.boundaryField() == pTraits<Type>::zero;

    res.correctBoundaryConditions();

    return tres;
}


template<class Type>
Foam::tmp<Foam::volScalarField> Foam::errorEstimate<Type>::normFactor() const
{
    tmp<volScalarField> tnormFactor
    (
        new volScalarField
        (
            IOobject
            (
                "normFactor" + psi_.name(),
                psi_.mesh().time().timeName(),
                psi_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            dimless/dimTime,
            errorBCTypes()
        )
    );

    volScalarField& normFactor = tnormFactor();

    normFactor.internalField() = normFactor_;
    normFactor.boundaryField() == pTraits<Type>::zero;

    normFactor.correctBoundaryConditions();

    return tnormFactor;
}

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::errorEstimate<Type>::error() const
{
    tmp<GeometricField<Type, fvPatchField, volMesh> > tresError
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "resError" + psi_.name(),
                psi_.mesh().time().timeName(),
                psi_.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            psi_.dimensions(),
            errorBCTypes()
        )
    );

    GeometricField<Type, fvPatchField, volMesh> resError = tresError();

    resError.internalField() == residual_/normFactor_;
    resError.boundaryField() = pTraits<Type>::zero;	// [ 0 0 0 ] 

    resError.correctBoundaryConditions();

    return tresError;
}
// * * * * * * * * * * * * * * *     clone         * * * * * * * * * * * * * //
//- Construct and return a clone
template<class Type>
Foam::tmp<Foam::errorEstimate<Type> > Foam::errorEstimate<Type>::clone() const
        {

	tmp<errorEstimate<Type> > terror
	(
		new errorEstimate<Type>(*this)
	);

	return terror;

        } 
// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //

template<class Type>
void Foam::errorEstimate<Type>::operator=(const Foam::errorEstimate<Type>& rhs)
{
    // Check for assignment to self
    if (this == &rhs)
    {
        FatalErrorIn
        (
            "errorEstimate<Type>::operator=(const Foam::errorEstimate<Type>&)"
        )   << "Attempted assignment to self"
            << abort(FatalError);
    }

    if (&psi_ != &(rhs.psi_))
    {
        FatalErrorIn
        (
            "errorEstimate<Type>::operator=(const errorEstimate<Type>&)"
        )   << "different fields"
            << abort(FatalError);
    }

    residual_ = rhs.residual_;
    normFactor_ = rhs.normFactor_;
}


template<class Type>
void Foam::errorEstimate<Type>::operator=(const tmp<errorEstimate<Type> >& teev)
{
    operator=(teev());
    teev.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::negate()
{
    residual_.negate();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=(const errorEstimate<Type>& eev)
{
    checkMethod(*this, eev, "+=");

    dimensions_ += eev.dimensions_;

    residual_ += eev.residual_;
    normFactor_ += eev.normFactor_;
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const tmp<errorEstimate<Type> >& teev
)
{
    operator+=(teev());
    teev.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=(const errorEstimate<Type>& eev)
{
    checkMethod(*this, eev, "+=");

    dimensions_ -= eev.dimensions_;
    residual_ -= eev.residual_;
    normFactor_ += eev.normFactor_;
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=(const tmp<errorEstimate<Type> >& teev)
{
    operator-=(teev());
    teev.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(*this, su, "+=");
    residual_ -= su.internalField();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    operator+=(tsu());
    tsu.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=
(
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(*this, su, "-=");
    residual_ += su.internalField();
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    operator-=(tsu());
    tsu.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator+=
(
    const dimensioned<Type>& su
)
{
    residual_ -= su;
}


template<class Type>
void Foam::errorEstimate<Type>::operator-=
(
    const dimensioned<Type>& su
)
{
    residual_ += su;
}


template<class Type>
void Foam::errorEstimate<Type>::operator*=
(
    const volScalarField& vsf
)
{
    dimensions_ *= vsf.dimensions();
    residual_ *= vsf.internalField();
    normFactor_ *= vsf.internalField();
}


template<class Type>
void Foam::errorEstimate<Type>::operator*=
(
    const tmp<volScalarField>& tvsf
)
{
    operator*=(tvsf());
    tvsf.clear();
}


template<class Type>
void Foam::errorEstimate<Type>::operator*=
(
    const dimensioned<scalar>& ds
)
{
    dimensions_ *= ds.dimensions();
    residual_ *= ds.value();
    normFactor_ *= ds.value();
}


// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //

template<class Type>
void Foam::checkMethod
(
    const errorEstimate<Type>& ee1,
    const errorEstimate<Type>& ee2,
    const char* op
)
{
    if (&ee1.psi() != &ee2.psi())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, "
            "const errorEstimate<Type>&)"
        )   << "incompatible fields for operation "
            << endl << "    "
            << "[" << ee1.psi().name() << "] "
            << op
            << " [" << ee2.psi().name() << "]"
            << abort(FatalError);
    }

    if (dimensionSet::debug && ee1.dimensions() != ee2.dimensions())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, "
            "const errorEstimate<Type>&)"
        )   << "incompatible dimensions for operation "
            << endl << "    "
            << "[" << ee1.psi().name() << ee1.dimensions()/dimVolume << " ] "
            << op
            << " [" << ee2.psi().name() << ee2.dimensions()/dimVolume << " ]"
            << abort(FatalError);
    }
}


template<class Type>
void Foam::checkMethod
(
    const errorEstimate<Type>& ee,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const char* op
)
{
    if (dimensionSet::debug && ee.dimensions()/dimVolume != vf.dimensions())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, "
            "const GeometricField<Type, fvPatchField, volMesh>&)"
        )   <<  "incompatible dimensions for operation "
            << endl << "    "
            << "[" << ee.psi().name() << ee.dimensions()/dimVolume << " ] "
            << op
            << " [" << vf.name() << vf.dimensions() << " ]"
            << abort(FatalError);
    }
}


template<class Type>
void Foam::checkMethod
(
    const errorEstimate<Type>& ee,
    const dimensioned<Type>& dt,
    const char* op
)
{
    if (dimensionSet::debug && ee.dimensions()/dimVolume != dt.dimensions())
    {
        FatalErrorIn
        (
            "checkMethod(const errorEstimate<Type>&, const dimensioned<Type>&)"
        )   << "incompatible dimensions for operation "
            << endl << "    "
            << "[" << ee.psi().name() << ee.dimensions()/dimVolume << " ] "
            << op
            << " [" << dt.name() << dt.dimensions() << " ]"
            << abort(FatalError);
    }
}


// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //

namespace Foam
{

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const errorEstimate<Type>& B
)
{
    checkMethod(A, B, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() += B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const errorEstimate<Type>& B
)
{
    checkMethod(tA(), B, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() += B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(A, tB(), "+");
    tmp<errorEstimate<Type> > tC(tB.ptr());
    tC() += A;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(tA(), tB(), "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() += tB();
    tB.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const errorEstimate<Type>& B
)
{
    checkMethod(A, B, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() -= B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const errorEstimate<Type>& B
)
{
    checkMethod(tA(), B, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() -= B;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(A, tB(), "-");
    tmp<errorEstimate<Type> > tC(tB.ptr());
    tC() -= A;
    tC().negate();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(tA(), tB(), "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() -= tB();
    tB.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const errorEstimate<Type>& B
)
{
    checkMethod(A, B, "==");
    return (A - B);
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const errorEstimate<Type>& B
)
{
    checkMethod(tA(), B, "==");
    return (tA - B);
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(A, tB(), "==");
    return (A - tB);
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<errorEstimate<Type> >& tB
)
{
    checkMethod(tA(), tB(), "==");
    return (tA - tB);
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(A, tsu(), "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(tA(), tsu(), "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, tsu(), "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), tsu(), "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(A, tsu(), "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(tA(), tsu(), "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    tC().res() -= su.internalField();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const GeometricField<Type, fvPatchField, volMesh>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    tC().res() -= su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, tsu(), "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), tsu(), "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    tC().res() -= tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const errorEstimate<Type>& A,
    const dimensioned<Type>& su
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const tmp<errorEstimate<Type> >& tA,
    const dimensioned<Type>& su
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const dimensioned<Type>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "+");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator+
(
    const dimensioned<Type>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "+");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const errorEstimate<Type>& A,
    const dimensioned<Type>& su
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const tmp<errorEstimate<Type> >& tA,
    const dimensioned<Type>& su
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const dimensioned<Type>& su,
    const errorEstimate<Type>& A
)
{
    checkMethod(A, su, "-");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().negate();
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator-
(
    const dimensioned<Type>& su,
    const tmp<errorEstimate<Type> >& tA
)
{
    checkMethod(tA(), su, "-");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().negate();
    tC().res() -= su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(A, su, "==");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const GeometricField<Type, fvPatchField, volMesh>& su
)
{
    checkMethod(tA(), su, "==");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.internalField();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(A, tsu(), "==");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
)
{
    checkMethod(tA(), tsu(), "==");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += tsu().internalField();
    tsu.clear();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const errorEstimate<Type>& A,
    const dimensioned<Type>& su
)
{
    checkMethod(A, su, "==");
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator==
(
    const tmp<errorEstimate<Type> >& tA,
    const dimensioned<Type>& su
)
{
    checkMethod(tA(), su, "==");
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC().res() += su.value();
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField& vsf,
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() *= vsf;
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>& tvsf,
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() *= tvsf;
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const volScalarField& vsf,
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() *= vsf;
    return tC;
}

template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const tmp<volScalarField>& tvsf,
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() *= tvsf;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>& ds,
    const errorEstimate<Type>& A
)
{
    tmp<errorEstimate<Type> > tC(new errorEstimate<Type>(A));
    tC() *= ds;
    return tC;
}


template<class Type>
tmp<errorEstimate<Type> > operator*
(
    const dimensioned<scalar>& ds,
    const tmp<errorEstimate<Type> >& tA
)
{
    tmp<errorEstimate<Type> > tC(tA.ptr());
    tC() *= ds;
    return tC;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //

 

  • resError.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Namespace
    Foam::resError

Description
    Namespace for residual error estimate operators.

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

#ifndef resError_H
#define resError_H

#include "resErrorDiv.H"
#include "resErrorLaplacian.H"
#include "resErrorSup.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

 

  • resErrorDiv.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

InNamespace
    Foam::resError

Description
    Residual error estimate for the fv convection operators.

SourceFiles
    resErrorDiv.C

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

#ifndef resErrorDiv_H
#define resErrorDiv_H

#include "errorEstimate.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

namespace resError
{
    // Divergence terms

        template<class Type>
        tmp<errorEstimate<Type> > div
        (
            const surfaceScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > div
        (
            const tmp<surfaceScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

}  // End namespace resError


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "resErrorDiv.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

 

  • resErrorDiv.C
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description
    Residual error estimate for the fv convection operators.

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

#include "resErrorDiv.H"
#include "fvc.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace resError
{

template<class Type>
tmp<errorEstimate<Type> >
div
(
    const surfaceScalarField& flux,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    const scalarField& vols = mesh.V();
    const surfaceVectorField& faceCentres = mesh.Cf();
    const volVectorField& cellCentres = mesh.C();
    const fvPatchList& patches = mesh.boundary();
    const unallocLabelList& owner = mesh.owner();
    const unallocLabelList& neighbour = mesh.neighbour();

    Field<Type> res(vols.size(), pTraits<Type>::zero);
    scalarField aNorm(vols.size(), 0.0);

    // Get sign of flux
    const surfaceScalarField signF = pos(flux);

    // Calculate gradient of the solution
    // Change of return type due to gradient cacheing.  HJ, 22/Apr/2016
    const tmp
    <
        GeometricField
        <
            typename outerProduct<vector, Type>::type,
            fvPatchField,
            volMesh
        >
    > tgradVf = fvc::grad(vf);

    const GeometricField
    <
        typename outerProduct<vector, Type>::type,
        fvPatchField,
        volMesh
    >& gradVf = tgradVf();

    // Internal faces
    forAll (owner, faceI)
    {
        // Calculate the centre of the face
        const vector& curFaceCentre = faceCentres[faceI];

        // Owner
        vector ownD = curFaceCentre - cellCentres[owner[faceI]];

        // Subtract convection
        res[owner[faceI]] -=
        (
            vf[owner[faceI]]
          + (ownD & gradVf[owner[faceI]])
        )*flux[faceI];

        aNorm[owner[faceI]] += signF[faceI]*flux[faceI];

        // Neighbour
        vector neiD = curFaceCentre - cellCentres[neighbour[faceI]];

        // Subtract convection
        res[neighbour[faceI]] +=
        (
            vf[neighbour[faceI]]
          + (neiD & gradVf[neighbour[faceI]])
        )*flux[faceI];

        aNorm[neighbour[faceI]] -= (1.0 - signF[faceI])*flux[faceI];
    }

    forAll (patches, patchI)
    {
        const vectorField& patchFaceCentres =
            faceCentres.boundaryField()[patchI];

        const scalarField& patchFlux = flux.boundaryField()[patchI];
        const scalarField& patchSignFlux = signF.boundaryField()[patchI];

        const labelList& fCells = patches[patchI].faceCells();

        forAll (fCells, faceI)
        {
            vector d =
                patchFaceCentres[faceI] - cellCentres[fCells[faceI]];

            // Subtract convection
            res[fCells[faceI]] -=
            (
                vf[fCells[faceI]]
              + (d & gradVf[fCells[faceI]])
            )*patchFlux[faceI];

            aNorm[fCells[faceI]] += patchSignFlux[faceI]*patchFlux[faceI];
        }
    }

    res /= vols;
    aNorm /= vols;

    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            flux.dimensions()*vf.dimensions(),
            res,
            aNorm
        )
    );
}


template<class Type>
tmp<errorEstimate<Type> >
div
(
    const tmp<surfaceScalarField>& tflux,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Div(resError::div(tflux(), vf));
    tflux.clear();
    return Div;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace resError

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 

  • resErrorLaplacian.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

InNamespace
    Foam::resError

Description
    Residual error estimate for the fv laplacian operators

SourceFiles
    resErrorLaplacian.C

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

#ifndef resErrorLaplacian_H
#define resErrorLaplacian_H

#include "errorEstimate.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

namespace resError
{
    // Laplacian terms

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const dimensionedScalar&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const volScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<volScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const surfaceScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<surfaceScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const volTensorField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<volTensorField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const surfaceTensorField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > laplacian
        (
            const tmp<surfaceTensorField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "resErrorLaplacian.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

 

  • resErrorLaplacian.C
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description
    Residual error estimate for the fv laplacian operators.

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

#include "resErrorLaplacian.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace resError
{

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    surfaceScalarField Gamma
    (
        IOobject
        (
            "gamma",
            vf.time().constant(),
            vf.db(),
            IOobject::NO_READ
        ),
        vf.mesh(),
        dimensionedScalar("1", dimless, 1.0)
    );

    return resError::laplacian(Gamma, vf);
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const dimensionedScalar& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    surfaceScalarField Gamma
    (
        IOobject
        (
            gamma.name(),
            vf.time().timeName(),
            vf.db(),
            IOobject::NO_READ
        ),
        vf.mesh(),
        gamma
    );

    return resError::laplacian(Gamma, vf);
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const volScalarField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return resError::laplacian(fvc::interpolate(gamma), vf);
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<volScalarField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Laplacian(resError::laplacian(tgamma(), vf));
    tgamma.clear();
    return Laplacian;
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const surfaceScalarField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    const scalarField& vols = mesh.V();
    const surfaceVectorField& Sf = mesh.Sf();
    const surfaceScalarField magSf = mesh.magSf();
    const fvPatchList& patches = mesh.boundary();
    const unallocLabelList& owner = mesh.owner();
    const unallocLabelList& neighbour = mesh.neighbour();

    const surfaceScalarField& delta =
        mesh.surfaceInterpolation::deltaCoeffs();

    Field<Type> res(vols.size(), pTraits<Type>::zero);
    scalarField aNorm(vols.size(), 0.0);

    // Calculate gradient of the solution.
    // Change of return type due to gradient cacheing.  HJ, 22/Apr/2016
    const tmp
    <
        GeometricField
        <
            typename outerProduct<vector, Type>::type,
            fvPatchField,
            volMesh
        >
    > tgradVf = fvc::grad(vf);

    const GeometricField
    <
        typename outerProduct<vector, Type>::type,
        fvPatchField,
        volMesh
    >& gradVf = tgradVf();

    // Internal faces
    forAll (owner, faceI)
    {
        // Owner

        // Subtract diffusion
        res[owner[faceI]] -=
            gamma[faceI]*(Sf[faceI] & gradVf[owner[faceI]]);

        aNorm[owner[faceI]] += delta[faceI]*gamma[faceI]*magSf[faceI];

        // Neighbour

        // Subtract diffusion
        res[neighbour[faceI]] +=
            gamma[faceI]*(Sf[faceI] & gradVf[neighbour[faceI]]);

        aNorm[neighbour[faceI]] += delta[faceI]*gamma[faceI]*magSf[faceI];

    }

    forAll (patches, patchI)
    {
        const vectorField& patchSf = Sf.boundaryField()[patchI];
        const scalarField& patchMagSf = magSf.boundaryField()[patchI];
        const scalarField& patchGamma = gamma.boundaryField()[patchI];
        const scalarField& patchDelta = delta.boundaryField()[patchI];

        const labelList& fCells = patches[patchI].faceCells();

        forAll (fCells, faceI)
        {
            // Subtract diffusion
            res[fCells[faceI]] -=
                patchGamma[faceI]*
                (
                    patchSf[faceI] & gradVf[fCells[faceI]]
                );

            aNorm[fCells[faceI]] +=
                patchDelta[faceI]*patchGamma[faceI]*patchMagSf[faceI];
        }
    }

    res /= vols;
    aNorm /= vols;

    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            delta.dimensions()*gamma.dimensions()*magSf.dimensions()
            *vf.dimensions(),
            res,
            aNorm
        )
    );
}

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<surfaceScalarField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > tresError(resError::laplacian(tgamma(), vf));
    tgamma.clear();
    return tresError;
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const volTensorField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    return resError::laplacian
    (
        (mesh.Sf() & fvc::interpolate(gamma) & mesh.Sf())
        /sqr(mesh.magSf()),
        vf
    );
}

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<volTensorField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Laplacian = resError::laplacian(tgamma(), vf);
    tgamma.clear();
    return Laplacian;
}


template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const surfaceTensorField& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    return resError::laplacian
    (
        (mesh.Sf() & gamma & mesh.Sf())/sqr(mesh.magSf()),
        vf
    );
}

template<class Type>
tmp<errorEstimate<Type> >
laplacian
(
    const tmp<surfaceTensorField>& tgamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > Laplacian = resError::laplacian(tgamma(), vf);
    tgamma.clear();
    return Laplacian;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace resError

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

 

  • resErrorSup.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

InNamespace
    Foam::resError

Description
    Residual error estimate for the fv source operators

SourceFiles
    resErrorSup.C

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

#ifndef resErrorSup_H
#define resErrorSup_H

#include "errorEstimate.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{


namespace resError
{
    // Implicit source

        template<class Type>
        tmp<errorEstimate<Type> > Sp
        (
            const volScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > Sp
        (
            const tmp<volScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );


        template<class Type>
        tmp<errorEstimate<Type> > Sp
        (
            const dimensionedScalar&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );


    // Implicit/Explicit source depending on sign of coefficient

        template<class Type>
        tmp<errorEstimate<Type> > SuSp
        (
            const volScalarField&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        template<class Type>
        tmp<errorEstimate<Type> > SuSp
        (
            const tmp<volScalarField>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#ifdef NoRepository
#   include "resErrorSup.C"
#endif

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

 

  • resErrorSup.C
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description
    Residual error estimate for the fv source operators.

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

#include "resErrorSup.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace resError
{

template<class Type>
tmp<errorEstimate<Type> >
Sp
(
    const volScalarField& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            sp.dimensions()*vf.dimensions()*dimVolume,
            sp.internalField()*vf.internalField(),
            scalarField(vf.internalField().size(), 0)
        )
    );
}

template<class Type>
tmp<errorEstimate<Type> >
Sp
(
    const tmp<volScalarField>& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > tee = resError::Sp(tsp(), vf);
    tsp.clear();
    return tee;
}


template<class Type>
tmp<errorEstimate<Type> >
Sp
(
    const dimensionedScalar& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return tmp<errorEstimate<Type> >
    (
        new errorEstimate<Type>
        (
            vf,
            sp.dimensions()*vf.dimensions()*dimVolume,
            sp.value()*vf.internalField(),
            scalarField(vf.internalField().size(), 0)
        )
    );
}


template<class Type>
tmp<errorEstimate<Type> >
SuSp
(
    const volScalarField& sp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return Sp(sp, vf);
}

template<class Type>
tmp<errorEstimate<Type> >
SuSp
(
    const tmp<volScalarField>& tsp,
    GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<errorEstimate<Type> > tee = resError::SuSp(tsp(), vf);
    tsp.clear();
    return tee;
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace resError

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //