Work on relaxation factors

This commit is contained in:
Henrik Rusche 2016-05-27 15:10:00 +02:00
parent fd95374293
commit f89267604a
14 changed files with 272 additions and 111 deletions

View file

@ -10,7 +10,7 @@
// Get under-relaxation factor // Get under-relaxation factor
scalar UUrf = scalar UUrf =
mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter())); mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()));
if (pimple.momentumPredictor()) if (pimple.momentumPredictor())
{ {

View file

@ -10,7 +10,7 @@ fvVectorMatrix ddtUEqn(fvm::ddt(U));
// Get under-relaxation factor // Get under-relaxation factor
scalar UUrf = scalar UUrf =
mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter())); mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()));
if (pimple.momentumPredictor()) if (pimple.momentumPredictor())
{ {

View file

@ -7,7 +7,7 @@
); );
// Get under-relaxation factor // Get under-relaxation factor
const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name()); const scalar UUrf = mesh.solutionDict().equationRelaxationFactor(U.name());
// Momentum solution // Momentum solution
solve solve

View file

@ -8,7 +8,7 @@
); );
// Get under-relaxation factor // Get under-relaxation factor
const scalar UUrf = mesh.solutionDict().relaxationFactor(Urel.name()); const scalar UUrf = mesh.solutionDict().equationRelaxationFactor(Urel.name());
// Momentum solution // Momentum solution
solve(relax(HUrelEqn(), UUrf) == -fvc::grad(p)); solve(relax(HUrelEqn(), UUrf) == -fvc::grad(p));

View file

@ -548,16 +548,18 @@ void faMatrix<Type>::relax(const scalar alpha)
template<class Type> template<class Type>
void faMatrix<Type>::relax() void faMatrix<Type>::relax()
{ {
scalar alpha = 0;
if (psi_.mesh().solutionDict().relax(psi_.name())) if (psi_.mesh().solutionDict().relax(psi_.name()))
{ {
alpha = psi_.mesh().solutionDict().relaxationFactor(psi_.name()); relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()));
} }
else
if (alpha > 0)
{ {
relax(alpha); if (debug)
{
InfoIn("void faMatrix<Type>::relax()")
<< "Relaxation factor for field " << psi_.name()
<< " not found. Relaxation will not be used." << endl;
}
} }
} }

View file

@ -45,7 +45,11 @@ void Foam::solutionControl::storePrevIter() const
size_t prevIterField = fName.find("PrevIter"); size_t prevIterField = fName.find("PrevIter");
if ((prevIterField == word::npos) && mesh_.solutionDict().relax(fName)) if
(
(prevIterField == word::npos)
&& mesh_.solutionDict().relaxField(fName)
)
{ {
if (debug) if (debug)
{ {

View file

@ -696,9 +696,9 @@ void Foam::fvMatrix<Type>::relax(const scalar alpha)
template<class Type> template<class Type>
void Foam::fvMatrix<Type>::relax() void Foam::fvMatrix<Type>::relax()
{ {
if (psi_.mesh().solutionDict().relax(psi_.name())) if (psi_.mesh().solutionDict().relaxEquation(psi_.name()))
{ {
relax(psi_.mesh().solutionDict().relaxationFactor(psi_.name())); relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()));
} }
else else
{ {

View file

@ -666,7 +666,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
if (debug) if (debug)
{ {
Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : " Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
"constructing as copy resetting IO params" "constructing as copy resetting IO params and patch types"
<< endl << this->info() << endl; << endl << this->info() << endl;
} }
@ -911,6 +911,20 @@ bool Foam::GeometricField<Type, PatchField, GeoMesh>::needReference() const
template<class Type, template<class> class PatchField, class GeoMesh> template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha) void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha)
{ {
if(alpha <= 0)
{
return;
}
//if (debug)
{
InfoIn
(
"GeometricField<Type, PatchField, GeoMesh>::relax"
"(const scalar alpha)"
) << "Relaxing" << endl << this->info() << " by " << alpha << endl;
}
operator==(prevIter() + alpha*(*this - prevIter())); operator==(prevIter() + alpha*(*this - prevIter()));
} }
@ -918,16 +932,11 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha)
template<class Type, template<class> class PatchField, class GeoMesh> template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::relax() void Foam::GeometricField<Type, PatchField, GeoMesh>::relax()
{ {
scalar alpha = 0; word name = this->name();
if (this->mesh().solutionDict().relax(this->name())) if (this->mesh().solutionDict().relaxField(name))
{ {
alpha = this->mesh().solutionDict().relaxationFactor(this->name()); relax(this->mesh().solutionDict().fieldRelaxationFactor(name));
}
if (alpha > 0)
{
relax(alpha);
} }
} }
@ -949,6 +958,19 @@ Foam::word Foam::GeometricField<Type, PatchField, GeoMesh>::select
} }
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::writeMinMax
(
Ostream& os
) const
{
os << "min/max(" << this->name() << ") = "
<< Foam::min(*this).value() << ", "
<< Foam::max(*this).value()
<< endl;
}
// writeData member function required by regIOobject // writeData member function required by regIOobject
template<class Type, template<class> class PatchField, class GeoMesh> template<class Type, template<class> class PatchField, class GeoMesh>
bool Foam::GeometricField<Type, PatchField, GeoMesh>:: bool Foam::GeometricField<Type, PatchField, GeoMesh>::

View file

@ -491,14 +491,17 @@ public:
// alpha = 0 : do nothing // alpha = 0 : do nothing
void relax(const scalar alpha); void relax(const scalar alpha);
//- Relax field (for steady-state solution).
// alpha is read from controlDict
void relax();
//- Select the final iteration parameters if 'final' is true //- Select the final iteration parameters if 'final' is true
// by returning the field name + "Final" // by returning the field name + "Final"
// otherwise the standard parameters by returning the field name // otherwise the standard parameters by returning the field name
word select(bool final) const; word select(bool final) const;
//- Relax field (for steady-state solution). //- Helper function to write the min and max to an Ostream
// alpha is read from controlDict void writeMinMax(Ostream& os) const;
void relax();
// Member function *this operators // Member function *this operators

View file

@ -49,6 +49,82 @@ static const Foam::List<Foam::word> subDictNames
//! @endcond localScope //! @endcond localScope
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::solution::read(const dictionary& dict)
{
if (dict.found("cache"))
{
cache_ = dict.subDict("cache");
caching_ = cache_.lookupOrDefault("active", true);
}
if (dict.found("relaxationFactors"))
{
const dictionary& relaxDict(dict.subDict("relaxationFactors"));
if (relaxDict.found("fields") || relaxDict.found("equations"))
{
if (relaxDict.found("fields"))
{
fieldRelaxDict_ = relaxDict.subDict("fields");
}
if (relaxDict.found("equations"))
{
eqnRelaxDict_ = relaxDict.subDict("equations");
}
}
else
{
// backwards compatibility
fieldRelaxDict_.clear();
const wordList entryNames(relaxDict.toc());
forAll(entryNames, i)
{
const word& e = entryNames[i];
scalar value = readScalar(relaxDict.lookup(e));
if (e(0, 1) == "p")
{
fieldRelaxDict_.add(e, value);
}
else if (e.length() >= 3)
{
if (e(0, 3) == "rho")
{
fieldRelaxDict_.add(e, value);
}
}
}
eqnRelaxDict_ = relaxDict;
}
fieldRelaxDefault_ =
fieldRelaxDict_.lookupOrDefault<scalar>("default", 0.0);
eqnRelaxDefault_ =
eqnRelaxDict_.lookupOrDefault<scalar>("default", 0.0);
if (debug)
{
Info<< "relaxation factors:" << nl
<< "fields: " << fieldRelaxDict_ << nl
<< "equations: " << eqnRelaxDict_ << endl;
}
}
if (dict.found("solvers"))
{
solvers_ = dict.subDict("solvers");
upgradeSolverDict(solvers_);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solution::solution(const objectRegistry& obr, const fileName& dictName) Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
@ -66,12 +142,11 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
), ),
cache_(dictionary::null), cache_(dictionary::null),
caching_(false), caching_(false),
relaxationFactors_ fieldRelaxDict_(dictionary::null),
( eqnRelaxDict_(dictionary::null),
ITstream("relaxationFactors", tokenList())() fieldRelaxDefault_(0),
), eqnRelaxDefault_(0),
defaultRelaxationFactor_(0), solvers_(dictionary::null)
solvers_(ITstream("solvers", tokenList())())
{ {
if (!headerOk()) if (!headerOk())
{ {
@ -86,7 +161,7 @@ Foam::solution::solution(const objectRegistry& obr, const fileName& dictName)
} }
} }
read(); read(solutionDict());
} }
@ -170,43 +245,6 @@ Foam::label Foam::solution::upgradeSolverDict
} }
bool Foam::solution::read()
{
bool readOk = false;
if (headerOk())
{
readOk = regIOobject::read();
}
if (readOk)
{
if (found("cache"))
{
cache_ = subDict("cache");
caching_ = cache_.lookupOrDefault("active", true);
}
if (found("relaxationFactors"))
{
relaxationFactors_ = subDict("relaxationFactors");
}
relaxationFactors_.readIfPresent("default", defaultRelaxationFactor_);
if (found("solvers"))
{
solvers_ = subDict("solvers");
upgradeSolverDict(solvers_);
}
return true;
}
return readOk;
}
bool Foam::solution::cache(const word& name) const bool Foam::solution::cache(const word& name) const
{ {
if (caching_) if (caching_)
@ -226,41 +264,52 @@ bool Foam::solution::cache(const word& name) const
} }
bool Foam::solution::relax(const word& name) const bool Foam::solution::relaxField(const word& name) const
{ {
if (debug) if (debug)
{ {
Info<< "Find relax for " << name << endl; Info<< "Field relaxation factor for " << name
<< " is " << (fieldRelaxDict_.found(name) ? "set" : "unset")
<< endl;
} }
return return fieldRelaxDict_.found(name) || fieldRelaxDict_.found("default");
relaxationFactors_.found(name)
|| relaxationFactors_.found("default");
} }
Foam::scalar Foam::solution::relaxationFactor(const word& name) const bool Foam::solution::relaxEquation(const word& name) const
{ {
if (debug) if (debug)
{ {
Info<< "Lookup relaxationFactor for " << name << endl; Info<< "Find equation relaxation factor for " << name << endl;
} }
if (relaxationFactors_.found(name)) return eqnRelaxDict_.found(name) || eqnRelaxDict_.found("default");
{
return readScalar(relaxationFactors_.lookup(name));
} }
else if (defaultRelaxationFactor_ > SMALL)
Foam::scalar Foam::solution::fieldRelaxationFactor(const word& name) const
{ {
return defaultRelaxationFactor_; if (debug)
{
Info<< "Lookup variable relaxation factor for " << name << endl;
}
if (fieldRelaxDict_.found(name))
{
return readScalar(fieldRelaxDict_.lookup(name));
}
else if (fieldRelaxDefault_ > SMALL)
{
return fieldRelaxDefault_;
} }
else else
{ {
FatalIOErrorIn FatalIOErrorIn
( (
"Foam::solution::relaxationFactor(const word&)", "Foam::solution::fieldRelaxationFactor(const word&)",
relaxationFactors_ fieldRelaxDict_
) << "Cannot find relaxationFactor for '" << name ) << "Cannot find variable relaxation factor for '" << name
<< "' or a suitable default value." << "' or a suitable default value."
<< exit(FatalIOError); << exit(FatalIOError);
@ -269,6 +318,49 @@ Foam::scalar Foam::solution::relaxationFactor(const word& name) const
} }
Foam::scalar Foam::solution::equationRelaxationFactor(const word& name) const
{
if (debug)
{
Info<< "Lookup equation relaxation factor for " << name << endl;
}
if (eqnRelaxDict_.found(name))
{
return readScalar(eqnRelaxDict_.lookup(name));
}
else if (eqnRelaxDefault_ > SMALL)
{
return eqnRelaxDefault_;
}
else
{
FatalIOErrorIn
(
"Foam::solution::eqnRelaxationFactor(const word&)",
eqnRelaxDict_
) << "Cannot find equation relaxation factor for '" << name
<< "' or a suitable default value."
<< exit(FatalIOError);
return 0;
}
}
const Foam::dictionary& Foam::solution::solutionDict() const
{
if (found("select"))
{
return subDict(word(lookup("select")));
}
else
{
return *this;
}
}
const Foam::dictionary& Foam::solution::solverDict(const word& name) const const Foam::dictionary& Foam::solution::solverDict(const word& name) const
{ {
if (debug) if (debug)
@ -293,15 +385,22 @@ const Foam::dictionary& Foam::solution::solver(const word& name) const
} }
bool Foam::solution::read()
{
if (regIOobject::read())
{
read(solutionDict());
return true;
}
else
{
return false;
}
}
bool Foam::solution::writeData(Ostream& os) const bool Foam::solution::writeData(Ostream& os) const
{ {
// Write dictionaries
os << nl << "solvers";
solvers_.write(os, true);
os << nl << "relaxationFactors";
relaxationFactors_.write(os, true);
// Write direct entries of the solution dictionary // Write direct entries of the solution dictionary
// HJ, 16/Feb/2010 // HJ, 16/Feb/2010
dictionary::write(os, false); dictionary::write(os, false);

View file

@ -60,10 +60,16 @@ class solution
bool caching_; bool caching_;
//- Dictionary of relaxation factors for all the fields //- Dictionary of relaxation factors for all the fields
dictionary relaxationFactors_; dictionary fieldRelaxDict_;
//- Dictionary of relaxation factors for all the equations
dictionary eqnRelaxDict_;
//- Optional default relaxation factor for all the fields //- Optional default relaxation factor for all the fields
scalar defaultRelaxationFactor_; scalar fieldRelaxDefault_;
//- Optional default relaxation factor for all the equations
scalar eqnRelaxDefault_;
//- Dictionary of solver parameters for all the fields //- Dictionary of solver parameters for all the fields
dictionary solvers_; dictionary solvers_;
@ -117,10 +123,20 @@ public:
); );
//- Return true if the relaxation factor is given for the field //- Return true if the relaxation factor is given for the field
bool relax(const word& name) const; bool relaxField(const word& name) const;
//- Return true if the relaxation factor is given for the equation
bool relaxEquation(const word& name) const;
//- Return the relaxation factor for the given field //- Return the relaxation factor for the given field
scalar relaxationFactor(const word& name) const; scalar fieldRelaxationFactor(const word& name) const;
//- Return the relaxation factor for the given eqation
scalar equationRelaxationFactor(const word& name) const;
//- Return the selected sub-dictionary of solvers if the "select"
// keyword is given, otherwise return the complete dictionary
const dictionary& solutionDict() const;
//- Return the solver controls dictionary for the given field //- Return the solver controls dictionary for the given field
const dictionary& solverDict(const word& name) const; const dictionary& solverDict(const word& name) const;
@ -131,10 +147,16 @@ public:
// Edit // Edit
//- Return access to relaxation factors dictionary //- Return access to field relaxation factors dictionary
dictionary& relaxationFactors() dictionary& fieldRelaxationFactors()
{ {
return relaxationFactors_; return fieldRelaxDict_;
}
//- Return access to equation relaxation factors dictionary
dictionary& equationRelaxationFactors()
{
return eqnRelaxDict_;
} }
//- Return access to solvers dictionary //- Return access to solvers dictionary

View file

@ -62,7 +62,7 @@ bool tolerances::read()
{ {
if (regIOobject::read()) if (regIOobject::read())
{ {
word toleranceSetName(lookup("toleranceSet")); const word toleranceSetName(lookup("toleranceSet"));
const dictionary& toleranceSet(subDict(toleranceSetName)); const dictionary& toleranceSet(subDict(toleranceSetName));
if (toleranceSet.found("relaxationFactors")) if (toleranceSet.found("relaxationFactors"))

View file

@ -215,16 +215,18 @@ void tetFemMatrix<Type>::relax(const scalar alpha)
template<class Type> template<class Type>
void tetFemMatrix<Type>::relax() void tetFemMatrix<Type>::relax()
{ {
scalar alpha = 0;
if (psi_.mesh().solutionDict().relax(psi_.name())) if (psi_.mesh().solutionDict().relax(psi_.name()))
{ {
alpha = psi_.mesh().solutionDict().relaxationFactor(psi_.name()); relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()));
} }
else
if (alpha > 0)
{ {
relax(alpha); if (debug)
{
InfoIn("void tetFemMatrix<Type>::relax()")
<< "Relaxation factor for field " << psi_.name()
<< " not found. Relaxation will not be used." << endl;
}
} }
} }

View file

@ -73,13 +73,20 @@ SIMPLE
} }
relaxationFactors relaxationFactors
{
fields
{ {
p 0.3; p 0.3;
}
equations
{
U 0.7; U 0.7;
k 0.7; k 0.7;
epsilon 0.7; epsilon 0.7;
R 0.7; R 0.7;
nuTilda 0.7; nuTilda 0.7;
} }
}
// ************************************************************************* // // ************************************************************************* //