View Issue Details

IDProjectCategoryView StatusLast Update
0001824OpenFOAMBugpublic2015-11-10 09:10
Reportermichael2015 Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSOtherOS Version2.3.0
Summary0001824: solve ( ...) for a tmp<fvVectorMatrix> in case of a coupled matrix solver gives back an empty object of type solverPerformance
DescriptionDoes anybody know wy when I use a segregated approuch for the velocity solve ( ...) gives back a non empty opject of the type solverPerformance but when I solve it coupled solve (...) gives me back the standard constructor containing zeros for nItierations, initalResidual, finialResidual and so one.
TagsNo tags attached.

Activities

wyldckat

2015-10-18 18:30

updater  

fvMatrixSolve.C (9,806 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2014 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 "LduMatrix.H"
#include "diagTensorField.H"

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

template<class Type>
void Foam::fvMatrix<Type>::setComponentReference
(
    const label patchi,
    const label facei,
    const direction cmpt,
    const scalar value
)
{
    if (psi_.needReference())
    {
        if (Pstream::master())
        {
            internalCoeffs_[patchi][facei].component(cmpt) +=
                diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]];

            boundaryCoeffs_[patchi][facei].component(cmpt) +=
                diag()[psi_.mesh().boundary()[patchi].faceCells()[facei]]
               *value;
        }
    }
}


template<class Type>
Foam::solverPerformance Foam::fvMatrix<Type>::solve
(
    const dictionary& solverControls
)
{
    if (debug)
    {
        Info.masterStream(this->mesh().comm())
            << "fvMatrix<Type>::solve(const dictionary& solverControls) : "
               "solving fvMatrix<Type>"
            << endl;
    }

    label maxIter = -1;
    if (solverControls.readIfPresent("maxIter", maxIter))
    {
        if (maxIter == 0)
        {
            return solverPerformance();
        }
    }

    word type(solverControls.lookupOrDefault<word>("type", "segregated"));

    if (type == "segregated")
    {
        return solveSegregated(solverControls);
    }
    else if (type == "coupled")
    {
        return solveCoupled(solverControls);
    }
    else
    {
        FatalIOErrorIn
        (
            "fvMatrix<Type>::solve(const dictionary& solverControls)",
            solverControls
        )   << "Unknown type " << type
            << "; currently supported solver types are segregated and coupled"
            << exit(FatalIOError);

        return solverPerformance();
    }
}


template<class Type>
Foam::solverPerformance Foam::fvMatrix<Type>::solveSegregated
(
    const dictionary& solverControls
)
{
    if (debug)
    {
        Info.masterStream(this->mesh().comm())
            << "fvMatrix<Type>::solveSegregated"
               "(const dictionary& solverControls) : "
               "solving fvMatrix<Type>"
            << endl;
    }

    GeometricField<Type, fvPatchField, volMesh>& psi =
       const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);

    solverPerformance solverPerfVec
    (
        "fvMatrix<Type>::solveSegregated",
        psi.name()
    );

    scalarField saveDiag(diag());

    Field<Type> source(source_);

    // At this point include the boundary source from the coupled boundaries.
    // This is corrected for the implict part by updateMatrixInterfaces within
    // the component loop.
    addBoundarySource(source);

    typename Type::labelType validComponents
    (
        pow
        (
            psi.mesh().solutionD(),
            pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
        )
    );

    for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
    {
        if (validComponents[cmpt] == -1) continue;

        // copy field and source

        scalarField psiCmpt(psi.internalField().component(cmpt));
        addBoundaryDiag(diag(), cmpt);

        scalarField sourceCmpt(source.component(cmpt));

        FieldField<Field, scalar> bouCoeffsCmpt
        (
            boundaryCoeffs_.component(cmpt)
        );

        FieldField<Field, scalar> intCoeffsCmpt
        (
            internalCoeffs_.component(cmpt)
        );

        lduInterfaceFieldPtrsList interfaces =
            psi.boundaryField().scalarInterfaces();

        // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
        // bouCoeffsCmpt for the explicit part of the coupled boundary
        // conditions
        initMatrixInterfaces
        (
            bouCoeffsCmpt,
            interfaces,
            psiCmpt,
            sourceCmpt,
            cmpt
        );

        updateMatrixInterfaces
        (
            bouCoeffsCmpt,
            interfaces,
            psiCmpt,
            sourceCmpt,
            cmpt
        );

        solverPerformance solverPerf;

        // Solver call
        solverPerf = lduMatrix::solver::New
        (
            psi.name() + pTraits<Type>::componentNames[cmpt],
            *this,
            bouCoeffsCmpt,
            intCoeffsCmpt,
            interfaces,
            solverControls
        )->solve(psiCmpt, sourceCmpt, cmpt);

        if (solverPerformance::debug)
        {
            solverPerf.print(Info.masterStream(this->mesh().comm()));
        }

        solverPerfVec = max(solverPerfVec, solverPerf);
        solverPerfVec.solverName() = solverPerf.solverName();

        psi.internalField().replace(cmpt, psiCmpt);
        diag() = saveDiag;
    }

    psi.correctBoundaryConditions();

    psi.mesh().setSolverPerformance(psi.name(), solverPerfVec);

    return solverPerfVec;
}


template<class Type>
Foam::solverPerformance Foam::fvMatrix<Type>::solveCoupled
(
    const dictionary& solverControls
)
{
    if (debug)
    {
        Info.masterStream(this->mesh().comm())
            << "fvMatrix<Type>::solveCoupled"
               "(const dictionary& solverControls) : "
               "solving fvMatrix<Type>"
            << endl;
    }

    GeometricField<Type, fvPatchField, volMesh>& psi =
       const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);

    LduMatrix<Type, scalar, scalar> coupledMatrix(psi.mesh());
    coupledMatrix.diag() = diag();
    coupledMatrix.upper() = upper();
    coupledMatrix.lower() = lower();
    coupledMatrix.source() = source();

    addBoundaryDiag(coupledMatrix.diag(), 0);
    addBoundarySource(coupledMatrix.source(), false);

    coupledMatrix.interfaces() = psi.boundaryField().interfaces();
    coupledMatrix.interfacesUpper() = boundaryCoeffs().component(0);
    coupledMatrix.interfacesLower() = internalCoeffs().component(0);

    autoPtr<typename LduMatrix<Type, scalar, scalar>::solver>
    coupledMatrixSolver
    (
        LduMatrix<Type, scalar, scalar>::solver::New
        (
            psi.name(),
            coupledMatrix,
            solverControls
        )
    );

    SolverPerformance<Type> solverPerf
    (
        coupledMatrixSolver->solve(psi)
    );

    if (SolverPerformance<Type>::debug)
    {
        solverPerf.print(Info.masterStream(this->mesh().comm()));
    }

    psi.correctBoundaryConditions();

    //create a summary of the maximum residuals and iterations
    solverPerformance solverPerfSummary = summary(solverPerf);

    psi.mesh().setSolverPerformance(psi.name(), solverPerfSummary);

    return solverPerfSummary;
}


template<class Type>
Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver>
Foam::fvMatrix<Type>::solver()
{
    return solver
    (
        psi_.mesh().solverDict
        (
            psi_.select
            (
                psi_.mesh().data::template lookupOrDefault<bool>
                ("finalIteration", false)
            )
        )
    );
}


template<class Type>
Foam::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve()
{
    return solve
    (
        fvMat_.psi_.mesh().solverDict
        (
            fvMat_.psi_.select
            (
                fvMat_.psi_.mesh().data::template lookupOrDefault<bool>
                ("finalIteration", false)
            )
        )
    );
}


template<class Type>
Foam::solverPerformance Foam::fvMatrix<Type>::solve()
{
    return solve
    (
        psi_.mesh().solverDict
        (
            psi_.select
            (
                psi_.mesh().data::template lookupOrDefault<bool>
                ("finalIteration", false)
            )
        )
    );
}


template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::residual() const
{
    tmp<Field<Type> > tres(new Field<Type>(source_));
    Field<Type>& res = tres();

    addBoundarySource(res);

    // Loop over field components
    for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
    {
        scalarField psiCmpt(psi_.internalField().component(cmpt));

        scalarField boundaryDiagCmpt(psi_.size(), 0.0);
        addBoundaryDiag(boundaryDiagCmpt, cmpt);

        FieldField<Field, scalar> bouCoeffsCmpt
        (
            boundaryCoeffs_.component(cmpt)
        );

        res.replace
        (
            cmpt,
            lduMatrix::residual
            (
                psiCmpt,
                res.component(cmpt) - boundaryDiagCmpt*psiCmpt,
                bouCoeffsCmpt,
                psi_.boundaryField().scalarInterfaces(),
                cmpt
            )
        );
    }

    return tres;
}


// ************************************************************************* //
fvMatrixSolve.C (9,806 bytes)   

wyldckat

2015-10-18 18:30

updater  

SolverPerformance.C (6,577 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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 "SolverPerformance.H"
#include "IOstreams.H"

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

template<class Type>
typename Foam::SolverPerformance<Type>::componentType
Foam::SolverPerformance<Type>::maxInitialResidual() const
{
    componentType maxValue = pTraits<componentType>::min;

    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
    {
        if(maxValue < component(initialResidual_, cmpt))
        {
            maxValue = component(initialResidual_, cmpt);
        }
    }

    return maxValue;
}


template<class Type>
typename Foam::SolverPerformance<Type>::componentType
Foam::SolverPerformance<Type>::maxFinalResidual() const
{
    componentType maxValue = pTraits<componentType>::min;

    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
    {
        if(maxValue < component(finalResidual_, cmpt))
        {
            maxValue = component(finalResidual_, cmpt);
        }
    }

    return maxValue;
}


template<class Type>
bool Foam::SolverPerformance<Type>::checkSingularity
(
    const Type& wApA
)
{
    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
    {
        singular_[cmpt] =
            component(wApA, cmpt) < vsmall_;
    }

    return singular();
}


template<class Type>
bool Foam::SolverPerformance<Type>::singular() const
{
    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
    {
        if (!singular_[cmpt]) return false;
    }

    return true;
}


template<class Type>
bool Foam::SolverPerformance<Type>::checkConvergence
(
    const Type& Tolerance,
    const Type& RelTolerance
)
{
    if (debug >= 2)
    {
        Info<< solverName_
            << ":  Iteration " << noIterations_
            << " residual = " << finalResidual_
            << endl;
    }

    if
    (
        finalResidual_ < Tolerance
     || (
            RelTolerance
          > small_*pTraits<Type>::one
         && finalResidual_ < cmptMultiply(RelTolerance, initialResidual_)
        )
    )
    {
        converged_ = true;
    }
    else
    {
        converged_ = false;
    }

    return converged_;
}


template<class Type>
void Foam::SolverPerformance<Type>::print
(
    Ostream& os
) const
{
    for(direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
    {
        if (pTraits<Type>::nComponents == 1)
        {
            os  << solverName_ << ":  Solving for " << fieldName_;

        }
        else
        {
            os  << solverName_ << ":  Solving for "
                << word(fieldName_ + pTraits<Type>::componentNames[cmpt]);
        }

        if (singular_[cmpt])
        {
            os  << ":  solution singularity" << endl;
        }
        else
        {
            os  << ", Initial residual = " << component(initialResidual_, cmpt)
                << ", Final residual = " << component(finalResidual_, cmpt)
                << ", No Iterations " << noIterations_
                << endl;
        }
    }
}


template<class Type>
bool Foam::SolverPerformance<Type>::operator!=
(
    const SolverPerformance<Type>& sp
) const
{
    return
    (
        solverName()      != sp.solverName()
     || fieldName()       != sp.fieldName()
     || initialResidual() != sp.initialResidual()
     || finalResidual()   != sp.finalResidual()
     || nIterations()     != sp.nIterations()
     || converged()       != sp.converged()
     || singular()        != sp.singular()
    );
}


template<class Type>
Foam::SolverPerformance<typename Foam::pTraits<Type>::cmptType>
Foam::summary(const SolverPerformance<Type>& sp)
{
  return 
      SolverPerformance<typename pTraits<Type>::cmptType>
      (
          sp.solverName(),
          sp.fieldName(),
          sp.maxInitialResidual(),
          sp.maxFinalResidual(),
          sp.nIterations(),
          sp.converged(),
          sp.singular()
      );
}


template<class Type>
typename Foam::SolverPerformance<Type> Foam::max
(
    const typename Foam::SolverPerformance<Type>& sp1,
    const typename Foam::SolverPerformance<Type>& sp2
)
{
    return SolverPerformance<Type>
    (
        sp1.solverName(),
        sp1.fieldName_,
        max(sp1.initialResidual(), sp2.initialResidual()),
        max(sp1.finalResidual(), sp2.finalResidual()),
        max(sp1.nIterations(), sp2.nIterations()),
        sp1.converged() && sp2.converged(),
        sp1.singular() || sp2.singular()
    );
}


template<class Type>
Foam::Istream& Foam::operator>>
(
    Istream& is,
    typename Foam::SolverPerformance<Type>& sp
)
{
    is.readBeginList("SolverPerformance<Type>");
    is  >> sp.solverName_
        >> sp.fieldName_
        >> sp.initialResidual_
        >> sp.finalResidual_
        >> sp.noIterations_
        >> sp.converged_
        >> sp.singular_;
    is.readEndList("SolverPerformance<Type>");

    return is;
}


template<class Type>
Foam::Ostream& Foam::operator<<
(
    Ostream& os,
    const typename Foam::SolverPerformance<Type>& sp
)
{
    os  << token::BEGIN_LIST
        << sp.solverName_ << token::SPACE
        << sp.fieldName_ << token::SPACE
        << sp.initialResidual_ << token::SPACE
        << sp.finalResidual_ << token::SPACE
        << sp.noIterations_ << token::SPACE
        << sp.converged_ << token::SPACE
        << sp.singular_ << token::SPACE
        << token::END_LIST;

    return os;
}


// ************************************************************************* //
SolverPerformance.C (6,577 bytes)   

wyldckat

2015-10-18 18:30

updater  

SolverPerformance.H (8,252 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2012-2013 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/>.

Class
    Foam::SolverPerformance

Description
    SolverPerformance is the class returned by the LduMatrix solver
    containing performance statistics.

SourceFiles
    SolverPerformance.C

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

#ifndef SolverPerformance_H
#define SolverPerformance_H

#include "word.H"
#include "FixedList.H"

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

namespace Foam
{

// Forward declaration of friend functions and operators

template<class Type>
class SolverPerformance;

//- Return a summary version of the solver performance with the
//  maximum values
template<class Type>
SolverPerformance<typename Foam::pTraits<Type>::cmptType>
summary(const SolverPerformance<Type>& sp);

template<class Type>
SolverPerformance<Type> max
(
    const SolverPerformance<Type>&,
    const SolverPerformance<Type>&
);

template<class Type>
Istream& operator>>
(
    Istream&,
    SolverPerformance<Type>&
);

template<class Type>
Ostream& operator<<
(
    Ostream&,
    const SolverPerformance<Type>&
);


/*---------------------------------------------------------------------------*\
                       Class SolverPerformance Declaration
\*---------------------------------------------------------------------------*/

template<class Type>
class SolverPerformance
{
    word   solverName_;
    word   fieldName_;
    Type   initialResidual_;
    Type   finalResidual_;
    label  noIterations_;
    bool   converged_;
    FixedList<bool, pTraits<Type>::nComponents> singular_;

public:

    //Typedefs

        typedef typename pTraits<Type>::cmptType componentType;

    // Static data

        // Declare name of the class and its debug switch
        ClassName("SolverPerformance");

        //- Large Type for the use in solvers
        static const scalar great_;

        //- Small Type for the use in solvers
        static const scalar small_;

        //- Very small Type for the use in solvers
        static const scalar vsmall_;


    // Constructors

        SolverPerformance()
        :
            initialResidual_(pTraits<Type>::zero),
            finalResidual_(pTraits<Type>::zero),
            noIterations_(0),
            converged_(false),
            singular_(false)
        {}


        SolverPerformance
        (
            const word&  solverName,
            const word&  fieldName,
            const Type&  iRes = pTraits<Type>::zero,
            const Type&  fRes = pTraits<Type>::zero,
            const label  nIter = 0,
            const bool   converged = false,
            const bool   singular = false
        )
        :
            solverName_(solverName),
            fieldName_(fieldName),
            initialResidual_(iRes),
            finalResidual_(fRes),
            noIterations_(nIter),
            converged_(converged),
            singular_(singular)
        {}


    // Member functions

        //- Return solver name
        const word& solverName() const
        {
            return solverName_;
        }

        //- Return solver name
        word& solverName()
        {
            return solverName_;
        }


        //- Return field name
        const word& fieldName() const
        {
            return fieldName_;
        }


        //- Return initial residual
        const Type& initialResidual() const
        {
            return initialResidual_;
        }

        //- Return initial residual
        Type& initialResidual()
        {
            return initialResidual_;
        }

        //- Return the maximum initial residual from the list
        componentType maxInitialResidual() const;


        //- Return final residual
        const Type& finalResidual() const
        {
            return finalResidual_;
        }

        //- Return final residual
        Type& finalResidual()
        {
            return finalResidual_;
        }

        //- Return the maximum final residual from the list
        componentType maxFinalResidual() const;


        //- Return number of iterations
        label nIterations() const
        {
            return noIterations_;
        }

        //- Return number of iterations
        label& nIterations()
        {
            return noIterations_;
        }


        //- Has the solver converged?
        bool converged() const
        {
            return converged_;
        }

        //- Is the matrix singular?
        bool singular() const;

        //- Check, store and return convergence
        bool checkConvergence
        (
            const Type& tolerance,
            const Type& relTolerance
        );

        //- Singularity test
        bool checkSingularity(const Type& residual);

        //- Print summary of solver performance to the given stream
        void print(Ostream& os) const;

    // Member Operators

        bool operator!=(const SolverPerformance<Type>&) const;


    // Friend functions

        //- Return the element-wise maximum of two SolverPerformance<Type>s
        friend SolverPerformance<Type> max <Type>
        (
            const SolverPerformance<Type>&,
            const SolverPerformance<Type>&
        );


    // Ostream Operator

        friend Istream& operator>> <Type>
        (
            Istream&,
            SolverPerformance<Type>&
        );

        friend Ostream& operator<< <Type>
        (
            Ostream&,
            const SolverPerformance<Type>&
        );
};


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

} // End namespace Foam

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

#define makeSolverPerformance(Type)                                           \
                                                                              \
typedef Foam::SolverPerformance<Type>                                         \
    solverPerformance##Type;                                                  \
                                                                              \
defineNamedTemplateTypeNameAndDebug(solverPerformance##Type, 0);              \
                                                                              \
template<>                                                                    \
const scalar solverPerformance##Type::great_(1e20);                           \
                                                                              \
template<>                                                                    \
const scalar solverPerformance##Type::small_(1e-20);                          \
                                                                              \
template<>                                                                    \
const scalar solverPerformance##Type::vsmall_(VSMALL);                        \


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

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

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

#endif

// ************************************************************************* //
SolverPerformance.H (8,252 bytes)   

wyldckat

2015-10-18 18:45

updater   ~0005426

I was curious about this issue and I went looking into it. So the story seems to be something along these lines:

  The coupled solver was implemented just fine, but there was a limitation in the current framework that OpenFOAM uses for the matrix solvers: the solver performance can only report 1 initial residual and 1 final residual, but the coupled solver is able to report as many residuals as the components the equation has got.
  In the tutorial "incompressible/pisoFoam/ras/cavityCoupledU" test case, the solver performance is reported by the internal coupled solver as having 3 values for each residual value.

  The solution probably wasn't implemented back then, because a general implementation wasn't straight forward. Which is why the empty result was returned.


Either way, I was curious about this and I also needed to hone my skills in C++ templates, so attached is a proposition for this.


@Henry: Attached are the following 3 files:

 - fvMatrixSolve.C - for the folder "src/finiteVolume/fvMatrices/fvMatrix/" - it uses the new "summary" function that is defined in the other two files.

 - SolverPerformance.H and SolverPerformance.C - for the folder "src/OpenFOAM/matrices/LduMatrix/LduMatrix/" - the changes are as follows:

  - 2 new public methods named "maxInitialResidual()" and "maxFinalResidual()", which rely on the type defined in "pTraits<Type>::cmptType" for returning the maximum value from the respective arrays.

  - 1 new function named "summary" which is a convenience function for returning a summary instance of "SolverPerformance<scalar>", which technically is defined as "SolverPerformance<pTraits<Type>::cmptType>".

This isn't the most elegantly possible solution for this problem, but it's the best I could do for now.

henry

2015-10-18 19:02

manager   ~0005427

This is all pending the completion of the rewrite of the solvers in the new templated framework after which the residual information would be handled in a consistent manner. There are still some choices to be made, for example should the segregated approach be maintained? If so even with non-aligned symmetry planes, cyclics etc. for which the coupled approach would be better. There are many details still to be worked on for which significant funding will be required.

wyldckat

2015-10-18 19:19

updater   ~0005428

I expected as much, namely that this feature was only a precursor to something a lot greater.
But in the meantime, it's a bit sad that the returned value is an empty instance of "solverPerformance".

I didn't propose just one single function that did the summary, since it would look like a much bigger hack than the one I proposed... hence the two auxiliary methods. And it's consistent with how the results for the segregated solver are given for the U equation.

henry

2015-10-20 10:35

manager   ~0005443

I have reviewed your changes and I think this takes the code in the wrong direction. What is required is that fvMatrix<Type>::solve returns a SolverPerformance<Type> rather than a solverPerformance. This is more work but at least consistent with the next set of developments in this area of the code.

wyldckat

2015-10-20 10:44

updater   ~0005444

Last edited: 2015-10-20 10:47

Understood! The addition of the "max*Residual" methods are adding methods to the class that would only be temporary, and any other hack would just be a clunky way to postpone the problem. Therefore, if anything is done, might as well be in the right direction; and OpenFOAM-dev is for the next major version, which doesn't bother as much with retro-compatibility.

OK, then the previously attached code can at least help out those who are eager to use this feature.

I can look into using "SolverPerformance<Type>" in the upcoming weekend. It's a great opportunity to learn a few more details on some of the code's innermost workings.

henry

2015-10-20 11:52

manager   ~0005445

Agreed. For backward compatibility "SolverPerformance<Type>" can have a "max" method which returns a "SolverPerformance<scalar>" the equivalent of the current "solverPerformance" when returned from the segregated vector (or higher rank) solver.

wyldckat

2015-10-31 23:19

updater  

proposition_v2.tar.gz (16,237 bytes)

wyldckat

2015-10-31 23:33

updater   ~0005536

Attached is the tarball "proposition_v2.tar.gz" that has the following files (caution, it's meant to be unpacked within "OpenFOAM-dev"):

  applications/solvers/stressAnalysis/solidDisplacementFoam/solidDisplacementFoam.C
  applications/test/volField/Test-volField.C
  src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C
  src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H
  src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.C
  src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.H
  src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
  src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
  src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C


Summary of proposed changes:

  1. solidDisplacementFoam.C: Adapted to using the new "max" function:

    - initialResidual = DEqn.solve().initialResidual();
    + initialResidual = max(DEqn.solve()).initialResidual();

  2. Test-volField.C: The solver performance is now saved to an object and outputted to the screen. In addition, an unbalanced "symmTensor" source was added to the equation, so that the residuals and solution was not identical for all of the solved components.

  3. SolverPerformance.[CH]: New friendly function "max" that returns the maximum instance as "SolverPerformance<typename Foam::pTraits<Type>::cmptType>"

  4. solverPerformance.[CH]: Specialization added for "SolverPerformance<scalar> max<scalar>(...)".

  5. fvMatrix.[CH]: All "solverPerformance" occurrences changed to "SolverPerformance<Type>".

  6. fvMatrixSolve.C:

     - Most "solverPerformance" occurrences changed to "SolverPerformance<Type>".

     - The method "solveSegregated" now uses a temporary and somewhat convoluted way of an instance of "SolverPerformance<Type>". It's the cleanest implementation I could figure out. You'll see what I mean after unpacking and seeing the changes.

     - The method "solveCoupled" now returns the true "SolverPerformance<Type>" result, instead of an empty instance.

     - The call to "setSolverPerformance" is now done by using the new "max" function, i.e. currently only storing the maximum solver performance indicators.


I only did a few tests, namely with "Test-volField" with the respective test case and the tutorial case "incompressible/icoFoam/cavity". Everything built successfully.


The next step forward would be to update "src/OpenFOAM/meshes/data/data.*" to the generic "SolverPerformance<Type>", but I haven't figured it out yet.

After that, if I saw correctly, then only "lduMatrix" is pending a revision.

wyldckat

2015-11-01 00:26

updater  

proposition_v3.tar.gz (3,756 bytes)

wyldckat

2015-11-01 00:36

updater   ~0005537

Attached is another tarball "proposition_v3.tar.gz", which will work on top of "proposition_v2.tar.gz".

The additional files provided in this tarball are as follows:

  src/OpenFOAM/meshes/data/data.C
  src/OpenFOAM/meshes/data/data.H
  src/OpenFOAM/meshes/data/dataTemplates.C
  src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C

Summary of changes (continuing from the previous list):

  6.d. Removed the need to call "setSolverPerformance" with the "max" function.

  7. data.[CH]: The methods "setSolverPerformance" are now template based. The header now depends on the file "dataTemplates.C".

  8. dataTemplates.C: New file, which has the two template methods "setSolverPerformance".


If I'm not mistaken, this solves most of the issues for this report and pushes OpenFOAM towards the right evolutionary direction :)


The only thing that is bugging me is that the "singular_" variable in "SolverPerformance" isn't transferring properly between instances... should it be upgraded as well?

wyldckat

2015-11-02 21:49

updater   ~0005544

Yesterday I didn't find the time to do more tests, but today I've ran several tutorials using the script "Alltest" in the "tutorials" folder and it looks like the "proposition_v3" is still buggy :(
The error message:

        From function solution::solverDict(const word&)
        in file matrices/solution/solution.C at line 369
        Lookup solver for k
    smoothSolver: Solving for k, Initial residual = 1, Final residual = 0.00458206, No Iterations 1


    --> FOAM FATAL IO ERROR:
    wrong token type - expected Scalar, found on line 0 the punctuation token '('

    file: /home/ofuser/OpenFOAM/ofuser-dev/run/tutorialsTest/compressible/rhoPimpleFoam/ras/angledDuct/system/data.solverPerformance.U at line 0.

        From function operator>>(Istream&, Scalar&)
        in file lnInclude/Scalar.C at line 93.

    FOAM exiting


Looks like each major type will have to be stored in its own subdict block, to avoid these kinds of conflicts.

wyldckat

2015-11-08 20:44

updater  

proposition_v4.tar.gz (17,288 bytes)

wyldckat

2015-11-08 21:01

updater   ~0005586

OK, I think I got it down right this time. Attached is the file "proposition_v4.tar.gz" that provides all of the modified files and does not depend on the previous v2 and v3 attachments.

Again, the package is meant to be unpacked within "OpenFOAM-dev". The provided files in the package:

    applications/solvers/stressAnalysis/solidDisplacementFoam/solidDisplacementFoam.C
    applications/test/volField/Test-volField.C
    src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.C
    src/OpenFOAM/matrices/LduMatrix/LduMatrix/SolverPerformance.H
    src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.C
    src/OpenFOAM/matrices/LduMatrix/LduMatrix/solverPerformance.H
    src/OpenFOAM/meshes/data/data.C
    src/OpenFOAM/meshes/data/data.H
    src/OpenFOAM/meshes/data/dataTemplates.C
    src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
    src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
    src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C

Summary of proposed changes:

  1. solidDisplacementFoam.C: Adapted to using the new "max" function:

    - initialResidual = DEqn.solve().initialResidual();
    + initialResidual = DEqn.solve().max().initialResidual();

  2. Test-volField.C: The solver performance is now saved to an object and outputted to the screen. In addition, an unbalanced "symmTensor" source was added to the equation, so that the residuals and solution was not identical for all of the solved components.

  3. SolverPerformance.[CH]: 2 new methods:

    - "max()" that returns the maximum instance as "SolverPerformance<typename Foam::pTraits<Type>::cmptType>"

    - "replace(cmpt, sp)" is essentially the opposite of "max", which allows to replace a particular component by using a "SolverPerformance<typename Foam::pTraits<Type>::cmptType>".

  4. solverPerformance.[CH]: Specialization added for "SolverPerformance<scalar>::max()".

  5. fvMatrix.[CH]: All "solverPerformance" occurrences changed to "SolverPerformance<Type>".

  6. fvMatrixSolve.C:

     - Most "solverPerformance" occurrences changed to "SolverPerformance<Type>".

     - The critical part of the method "solveSegregated" was modified like this:

          - solverPerfVec = max(solverPerfVec, solverPerf);
          - solverPerfVec.solverName() = solverPerf.solverName();
          + solverPerfVec.replace(cmpt, solverPerf);

     - The method "solveCoupled" now returns the true "SolverPerformance<Type>" result, instead of an empty instance, and also stores the solver performance data into the dictionary.

  7. data.[CH]: A few modifications:

     - The methods "setSolverPerformance" are now template based. The header now depends on the file "dataTemplates.C".

     - And also added the method "solverPerformanceDict(const SolverPerformance<Type>&)" which returns a subdictionary into the main subdict, so that we store each type of solver performance inside a block for its own type; respective method specialization for "scalar" also added.

  8. dataTemplates.C: New file, which has the two template methods "setSolverPerformance" and the new template method "solverPerformanceDict(const SolverPerformance<Type>&)".


I've ran many of the tutorial cases by using "Alltest" in both OpenFOAM-dev (commit 2ac29110299, Oct 31) and OpenFOAM-2.4.x. There were no crashes like the ones I mentioned in the previous comment. And I didn't spot any major residual result discrepancies.

A few aesthetic details are probably missing, but I did the best I could do based on the instructions given at http://www.openfoam.org/contrib/code-style.php

henry

2015-11-09 09:26

manager   ~0005589

Excellent work Bruno. I am applying and testing these changes at the moment, if all goes well I should be able to push by the end of the day.

A couple of formatting tips:

    Set your editor to remove trailing whitespace automatically.
    This is trivial to setup in emacs and vim but I guess it can be done in all decent editors.

    Start comments with a space and capital letter i.e.
    // This....
    //not this...

One other point, could you update the copyright date range for file you change.

I have made the above changes to the files you provided so far so this is for future reference.

Thanks

henry

2015-11-10 09:10

manager   ~0005595

@Bruno: Thanks for the patches. In the form provided the residualControl and residual functionObject functionality no longer worked so I back-tracked on some of the changes and updated residualControl (backward-compatible) and the residual functionObject (upgraded to output the component residuals) and now all
the standard tests pass.

The changes are quite substantial and required a few core library changes to improve pTraits support so that all the types could be handled consistently including scalar without function specialization. You might find the way pTraits is used interesting.

Resolved in OpenFOAM-dev by commit 1944b09bb530a8b9a2546ac3b196b137668d329f

Issue History

Date Modified Username Field Change
2015-08-12 18:54 michael2015 New Issue
2015-10-18 18:30 wyldckat File Added: fvMatrixSolve.C
2015-10-18 18:30 wyldckat File Added: SolverPerformance.C
2015-10-18 18:30 wyldckat File Added: SolverPerformance.H
2015-10-18 18:45 wyldckat Note Added: 0005426
2015-10-18 18:46 wyldckat Assigned To => henry
2015-10-18 18:46 wyldckat Status new => assigned
2015-10-18 19:02 henry Note Added: 0005427
2015-10-18 19:19 wyldckat Note Added: 0005428
2015-10-20 10:35 henry Note Added: 0005443
2015-10-20 10:44 wyldckat Note Added: 0005444
2015-10-20 10:45 wyldckat Assigned To henry =>
2015-10-20 10:46 wyldckat Status assigned => new
2015-10-20 10:47 wyldckat Note Edited: 0005444
2015-10-20 11:52 henry Note Added: 0005445
2015-10-31 23:19 wyldckat File Added: proposition_v2.tar.gz
2015-10-31 23:33 wyldckat Note Added: 0005536
2015-11-01 00:26 wyldckat File Added: proposition_v3.tar.gz
2015-11-01 00:36 wyldckat Note Added: 0005537
2015-11-02 21:49 wyldckat Note Added: 0005544
2015-11-08 20:44 wyldckat File Added: proposition_v4.tar.gz
2015-11-08 21:01 wyldckat Note Added: 0005586
2015-11-08 21:01 wyldckat Assigned To => henry
2015-11-08 21:01 wyldckat Status new => assigned
2015-11-09 09:26 henry Note Added: 0005589
2015-11-10 09:10 henry Note Added: 0005595
2015-11-10 09:10 henry Status assigned => resolved
2015-11-10 09:10 henry Resolution open => fixed