Commit 9e38e143 authored by mattijs's avatar mattijs

ENH: exact restart. Fixes #1172.

parent 47dc4344
......@@ -53,11 +53,39 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcNut() const
const tmp<scalarField> tnuw = turbModel.nu(patchi);
const scalarField& nuw = tnuw();
return max
// Calculate new viscosity
tmp<scalarField> tnutw
(
scalar(0),
sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw
max
(
scalar(0),
sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw
)
);
if (tolerance_ != 0.01)
{
// User-specified tolerance. Find out if current nut already satisfies
// eqns.
// Run ut for one iteration
scalarField err;
tmp<scalarField> UTau(calcUTau(magGradU, 1, err));
// Preserve nutw if the initial error (err) already lower than the
// tolerance.
scalarField& nutw = tnutw.ref();
forAll(err, facei)
{
if (err[facei] < tolerance_)
{
nutw[facei] = this->operator[](facei);
}
}
}
return tnutw;
}
......@@ -65,6 +93,18 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
(
const scalarField& magGradU
) const
{
scalarField err;
return calcUTau(magGradU, maxIter_, err);
}
tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
(
const scalarField& magGradU,
const label maxIter,
scalarField& err
) const
{
const label patchi = patch().index();
......@@ -89,6 +129,9 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
tmp<scalarField> tuTau(new scalarField(patch().size(), Zero));
scalarField& uTau = tuTau.ref();
err.setSize(uTau.size());
err = 0.0;
forAll(uTau, facei)
{
scalar ut = sqrt((nutw[facei] + nuw[facei])*magGradU[facei]);
......@@ -98,7 +141,6 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
if (ut > ROOTVSMALL)
{
int iter = 0;
scalar err = GREAT;
do
{
......@@ -116,12 +158,17 @@ tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
+ 1/E_*kUu*fkUu/ut;
scalar uTauNew = ut + f/df;
err = mag((ut - uTauNew)/ut);
err[facei] = mag((ut - uTauNew)/ut);
ut = uTauNew;
//iterations_++;
} while (ut > ROOTVSMALL && err > tolerance_ && ++iter < maxIter_);
} while
(
ut > ROOTVSMALL
&& err[facei] > tolerance_
&& ++iter < maxIter
);
uTau[facei] = max(0.0, ut);
......
......@@ -71,9 +71,18 @@ Note
Suffers from non-exact restart since correctNut() (called through
turbulence->validate) returns a slightly different value every time
it is called. This is since the seed for the Newton-Raphson iteration
uses the current value of *this (= nut). Can be avoided by seeding the
NR with e.g. the laminar viscosity or tightening the convergence tolerance
to e.g. 1e-7 and the max number of iterations to 100.
uses the current value of *this (= nut).
This can be avoided by overriding the tolerance. This also switches on
a pre-detection whether the current nut already satisfies the turbulence
conditions and if so does not change it at all. This means that the nut
only changes if it 'has really changed'. This probably should be used with
a tight tolerance, e.g.
maxIter 100;
tolerance 1e-7;
to make sure to kick every iteration.
SourceFiles
nutUSpaldingWallFunctionFvPatchScalarField.C
......@@ -123,6 +132,15 @@ protected:
//- Calculate the friction velocity
virtual tmp<scalarField> calcUTau(const scalarField& magGradU) const;
//- Calculate the friction velocity and number of iterations for
// convergence
virtual tmp<scalarField> calcUTau
(
const scalarField& magGradU,
const label maxIter,
scalarField& err
) const;
//- Write local wall function variables
virtual void writeLocalEntries(Ostream&) const;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment