View Issue Details

IDProjectCategoryView StatusLast Update
0003338OpenFOAMBugpublic2019-08-25 19:55
Reporterwyldckat Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Summary0003338: porosityModel::adjustNegativeResistance doesn't properly check for when all values are negative
DescriptionThis is a somewhat peculiar issue, in method 'porosityModel::adjustNegativeResistance':

    scalar maxCmpt = max(0, cmptMax(resist.value()));

    if (maxCmpt < 0)
    {
        FatalErrorInFunction
            << "Cannot have all resistances set to negative, resistance = "
            << resist
            << exit(FatalError);
    }

Which means that 'maxCmpt' is always '>= 0', due to how the 'max' operation is being done. Therefore, someone could have used negative resistance values, gotten no resistance out of it and no fatal error accordingly...

This affects all versions of OpenFOAM, as far as I could search back into the past (seems to have been added somewhere between 1.4 and 1.5).

Attached at the following files as a proposal to fix this issue:

 - proposal_v1.patch - this provides the proposed changes, to properly both check if all values are not negative and then use the above 0 values only when needed.

 - porosityModel.C - modified file from OpenFOAM-dev to update 'src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C', which works on both OpenFOAM-dev, 7 and 6.
TagsNo tags attached.

Activities

wyldckat

2019-08-24 17:34

updater  

porosityModel.C (5,700 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2012-2019 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 "porosityModel.H"
#include "volFields.H"

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

namespace Foam
{
    defineTypeNameAndDebug(porosityModel, 0);
    defineRunTimeSelectionTable(porosityModel, mesh);
}


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

void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
{
    scalar maxCmpt = cmptMax(resist.value());

    if (maxCmpt < 0)
    {
        FatalErrorInFunction
            << "Cannot have all resistances set to negative, resistance = "
            << resist
            << exit(FatalError);
    }
    else
    {
        maxCmpt = max(0, maxCmpt);
        vector& val = resist.value();
        for (label cmpt = 0; cmpt < vector::nComponents; cmpt++)
        {
            if (val[cmpt] < 0)
            {
                val[cmpt] *= -maxCmpt;
            }
        }
    }
}


Foam::label Foam::porosityModel::fieldIndex(const label i) const
{
    label index = 0;
    if (!coordSys_.R().uniform())
    {
        index = i;
    }
    return index;
}


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

Foam::porosityModel::porosityModel
(
    const word& name,
    const word& modelType,
    const fvMesh& mesh,
    const dictionary& dict,
    const word& cellZoneName
)
:
    regIOobject
    (
        IOobject
        (
            name,
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        )
    ),
    name_(name),
    mesh_(mesh),
    dict_(dict),
    coeffs_(dict.optionalSubDict(modelType + "Coeffs")),
    active_(true),
    zoneName_(cellZoneName),
    cellZoneIDs_(),
    coordSys_(coordinateSystem::New(mesh, coeffs_))
{
    if (zoneName_ == word::null)
    {
        dict.readIfPresent("active", active_);
        dict_.lookup("cellZone") >> zoneName_;
    }

    cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);

    Info<< "    creating porous zone: " << zoneName_ << endl;

    bool foundZone = !cellZoneIDs_.empty();
    reduce(foundZone, orOp<bool>());

    if (!foundZone && Pstream::master())
    {
        FatalErrorInFunction
            << "cannot find porous cellZone " << zoneName_
            << exit(FatalError);
    }
}


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

Foam::porosityModel::~porosityModel()
{}


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

void Foam::porosityModel::transformModelData()
{
    if (!mesh_.upToDatePoints(*this))
    {
        calcTransformModelData();

        // set model up-to-date wrt points
        mesh_.setUpToDatePoints(*this);
    }
}


Foam::tmp<Foam::vectorField> Foam::porosityModel::porosityModel::force
(
    const volVectorField& U,
    const volScalarField& rho,
    const volScalarField& mu
)
{
    transformModelData();

    tmp<vectorField> tforce(new vectorField(U.size(), Zero));

    if (!cellZoneIDs_.empty())
    {
        this->calcForce(U, rho, mu, tforce.ref());
    }

    return tforce;
}


void Foam::porosityModel::addResistance(fvVectorMatrix& UEqn)
{
    if (cellZoneIDs_.empty())
    {
        return;
    }

    transformModelData();
    this->correct(UEqn);
}


void Foam::porosityModel::addResistance
(
    fvVectorMatrix& UEqn,
    const volScalarField& rho,
    const volScalarField& mu
)
{
    if (cellZoneIDs_.empty())
    {
        return;
    }

    transformModelData();
    this->correct(UEqn, rho, mu);
}


void Foam::porosityModel::addResistance
(
    const fvVectorMatrix& UEqn,
    volTensorField& AU,
    bool correctAUprocBC
)
{
    if (cellZoneIDs_.empty())
    {
        return;
    }

    transformModelData();
    this->correct(UEqn, AU);

    if (correctAUprocBC)
    {
        // Correct the boundary conditions of the tensorial diagonal to ensure
        // processor boundaries are correctly handled when AU^-1 is interpolated
        // for the pressure equation.
        AU.correctBoundaryConditions();
    }
}


bool Foam::porosityModel::writeData(Ostream& os) const
{
    return true;
}


bool Foam::porosityModel::read(const dictionary& dict)
{
    dict.readIfPresent("active", active_);

    coeffs_ = dict.optionalSubDict(type() + "Coeffs");

    dict.lookup("cellZone") >> zoneName_;
    cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);

    return true;
}


// ************************************************************************* //
porosityModel.C (5,700 bytes)   
proposal_v1.patch (1,486 bytes)   
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
index e491bb2..9b50a9d 100644
--- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,16 +39,18 @@ namespace Foam
 
 void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
 {
-    scalar maxCmpt = max(0, cmptMax(resist.value()));
+    scalar maxCmpt = cmptMax(resist.value());
 
     if (maxCmpt < 0)
     {
         FatalErrorInFunction
-            << "Negative resistances are invalid, resistance = " << resist
+            << "Cannot have all resistances set to negative, resistance = "
+            << resist
             << exit(FatalError);
     }
     else
     {
+        maxCmpt = max(0, maxCmpt);
         vector& val = resist.value();
         for (label cmpt = 0; cmpt < vector::nComponents; cmpt++)
         {
proposal_v1.patch (1,486 bytes)   

henry

2019-08-25 19:55

manager   ~0010698

Thanks Bruno
Resolved in OpenFOAM-7 by commit ce67ac7f21d1edab96957a01a0b9cb5fbaff76a3
Resolved in OpenFOAM-dev by commit 9847140a8f021957a9471168e52b123136ee51e6

Issue History

Date Modified Username Field Change
2019-08-24 17:34 wyldckat New Issue
2019-08-24 17:34 wyldckat Status new => assigned
2019-08-24 17:34 wyldckat Assigned To => henry
2019-08-24 17:34 wyldckat File Added: porosityModel.C
2019-08-24 17:34 wyldckat File Added: proposal_v1.patch
2019-08-25 19:55 henry Status assigned => resolved
2019-08-25 19:55 henry Resolution open => fixed
2019-08-25 19:55 henry Fixed in Version => 7
2019-08-25 19:55 henry Note Added: 0010698