zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2013-11-08 12:07:35 +00:00
parent 6840fa315e 3dec6f0101
commit 6359588079
31 changed files with 1345 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
View File

@ -54,6 +54,8 @@
{
phi = phiHbyA - p_rghEqn.flux();
p_rgh.relax();
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);

3
src/ODE/Make/files
View File

@ -9,10 +9,13 @@ ODESolvers/RKF45/RKF45.C
ODESolvers/RKCK45/RKCK45.C
ODESolvers/RKDP45/RKDP45.C
ODESolvers/Rosenbrock21/Rosenbrock21.C
ODESolvers/Rosenbrock32/Rosenbrock32.C
ODESolvers/Rosenbrock43/Rosenbrock43.C
ODESolvers/rodas32/rodas32.C
ODESolvers/rodas43/rodas43.C
ODESolvers/SIBS/SIBS.C
ODESolvers/SIBS/SIMPR.C
ODESolvers/SIBS/polyExtrapolate.C
ODESolvers/seulex/seulex.C
LIB = $(FOAM_LIBBIN)/libODE

4
src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C
View File

@ -42,8 +42,8 @@ const scalar
Rosenbrock21::b2 = (1.0/2.0)/gamma,
Rosenbrock21::e1 = b1 - 1.0/gamma,
Rosenbrock21::e2 = b2,
Rosenbrock21::d1 = 1,
Rosenbrock21::d2 = -1;
Rosenbrock21::d1 = gamma,
Rosenbrock21::d2 = -gamma;
}

161
src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C Normal file
View File

@ -0,0 +1,161 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 "Rosenbrock32.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(Rosenbrock32, 0);
addToRunTimeSelectionTable(ODESolver, Rosenbrock32, dictionary);
const scalar
Rosenbrock32::a21 = 1,
Rosenbrock32::a31 = 1,
Rosenbrock32::a32 = 0,
Rosenbrock32::c21 = -1.0156171083877702091975600115545,
Rosenbrock32::c31 = 4.0759956452537699824805835358067,
Rosenbrock32::c32 = 9.2076794298330791242156818474003,
Rosenbrock32::b1 = 1,
Rosenbrock32::b2 = 6.1697947043828245592553615689730,
Rosenbrock32::b3 = -0.4277225654321857332623837380651,
Rosenbrock32::e1 = 0.5,
Rosenbrock32::e2 = -2.9079558716805469821718236208017,
Rosenbrock32::e3 = 0.2235406989781156962736090927619,
Rosenbrock32::gamma = 0.43586652150845899941601945119356,
Rosenbrock32::c2 = 0.43586652150845899941601945119356,
Rosenbrock32::d1 = 0.43586652150845899941601945119356,
Rosenbrock32::d2 = 0.24291996454816804366592249683314,
Rosenbrock32::d3 = 2.1851380027664058511513169485832;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
:
ODESolver(ode, dict),
adaptiveSolver(ode, dict),
k1_(n_),
k2_(n_),
k3_(n_),
err_(n_),
dydx_(n_),
dfdx_(n_),
dfdy_(n_, n_),
a_(n_, n_),
pivotIndices_(n_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::Rosenbrock32::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y
) const
{
ode.jacobian(x0, y0, dfdx_, dfdy_);
for (register label i=0; i<n_; i++)
{
for (register label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
}
a_[i][i] += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
// Calculate k1:
forAll(k1_, i)
{
k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, k1_);
// Calculate k2:
forAll(y, i)
{
y[i] = y0[i] + a21*k1_[i];
}
ode.derivatives(x0 + c2*dx, y, dydx_);
forAll(k2_, i)
{
k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
}
LUBacksubstitute(a_, pivotIndices_, k2_);
// Calculate k3:
forAll(k3_, i)
{
k3_[i] = dydx_[i] + dx*d3*dfdx_[i]
+ (c31*k1_[i] + c32*k2_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, k3_);
// Calculate error and update state:
forAll(y, i)
{
y[i] = y0[i] + b1*k1_[i] + b2*k2_[i] + b3*k3_[i];
err_[i] = e1*k1_[i] + e2*k2_[i] + e3*k3_[i];
}
return normalizeError(y0, y, err_);
}
void Foam::Rosenbrock32::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
}
// ************************************************************************* //

135
src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H Normal file
View File

@ -0,0 +1,135 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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::Rosenbrock32
Description
L-stable embedded Rosenbrock ODE solver of order (3)4.
References:
\verbatim
Sandu et al,
"Benchmarking stiff ODE solvers for atmospheric chemistry problems II
Rosenbrock solvers",
A. Sandu,
J.G. Verwer,
J.G. Blom,
E.J. Spee,
G.R. Carmichael,
F.A. Potra,
Atmospheric Environment, Volume 31, 1997, Issue 20, Pages 3459-3472
\endverbatim
SourceFiles
Rosenbrock32.C
\*---------------------------------------------------------------------------*/
#ifndef Rosenbrock32_H
#define Rosenbrock32_H
#include "ODESolver.H"
#include "adaptiveSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class Rosenbrock32 Declaration
\*---------------------------------------------------------------------------*/
class Rosenbrock32
:
public ODESolver,
public adaptiveSolver
{
// Private data
mutable scalarField k1_;
mutable scalarField k2_;
mutable scalarField k3_;
mutable scalarField err_;
mutable scalarField dydx_;
mutable scalarField dfdx_;
mutable scalarSquareMatrix dfdy_;
mutable scalarSquareMatrix a_;
mutable labelList pivotIndices_;
static const scalar
a21, a31, a32,
c21, c31, c32,
b1, b2, b3,
e1, e2, e3,
gamma,
c2, c3,
d1, d2, d3;
public:
//- Runtime type information
TypeName("Rosenbrock32");
// Constructors
//- Construct from ODE
Rosenbrock32(const ODESystem& ode, const dictionary& dict);
// Member Functions
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y
) const;
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C
View File

@ -91,7 +91,7 @@ const scalar
Rosenbrock43::e3 = 0,
Rosenbrock43::e4 = 125.0/108.0,
Rosenbrock43::gamma = 0.5,
Rosenbrock43::gamma = 1.0/2.0,
Rosenbrock43::c2 = 1,
Rosenbrock43::c3 = 3.0/5.0,

166
src/ODE/ODESolvers/rodas32/rodas32.C Normal file
View File

@ -0,0 +1,166 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 "rodas32.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(rodas32, 0);
addToRunTimeSelectionTable(ODESolver, rodas32, dictionary);
const scalar
rodas32::c3 = 1,
rodas32::d1 = 1.0/2.0,
rodas32::d2 = 3.0/2.0,
rodas32::a31 = 2,
rodas32::a41 = 2,
rodas32::c21 = 4,
rodas32::c31 = 1,
rodas32::c32 = -1,
rodas32::c41 = 1,
rodas32::c42 = -1,
rodas32::c43 = -8.0/3.0,
rodas32::gamma = 1.0/2.0;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict)
:
ODESolver(ode, dict),
adaptiveSolver(ode, dict),
k1_(n_),
k2_(n_),
k3_(n_),
dy_(n_),
err_(n_),
dydx_(n_),
dfdx_(n_),
dfdy_(n_, n_),
a_(n_, n_),
pivotIndices_(n_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::rodas32::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y
) const
{
ode.jacobian(x0, y0, dfdx_, dfdy_);
for (register label i=0; i<n_; i++)
{
for (register label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
}
a_[i][i] += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
// Calculate k1:
forAll(k1_, i)
{
k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, k1_);
// Calculate k2:
forAll(k2_, i)
{
k2_[i] = dydx0[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
}
LUBacksubstitute(a_, pivotIndices_, k2_);
// Calculate k3:
forAll(y, i)
{
dy_[i] = a31*k1_[i];
y[i] = y0[i] + dy_[i];
}
ode.derivatives(x0 + dx, y, dydx_);
forAll(k3_, i)
{
k3_[i] = dydx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, k3_);
// Calculate new state and error
forAll(y, i)
{
dy_[i] += k3_[i];
y[i] = y0[i] + dy_[i];
}
ode.derivatives(x0 + dx, y, dydx_);
forAll(err_, i)
{
err_[i] = dydx_[i] + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, err_);
forAll(y, i)
{
y[i] = y0[i] + dy_[i] + err_[i];
}
return normalizeError(y0, y, err_);
}
void Foam::rodas32::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
}
// ************************************************************************* //

135
src/ODE/ODESolvers/rodas32/rodas32.H Normal file
View File

@ -0,0 +1,135 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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::rodas32
Description
L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (2)3.
References:
\verbatim
Sandu et al,
"Benchmarking stiff ODE solvers for atmospheric chemistry problems II
Rosenbrock solvers",
A. Sandu,
J.G. Verwer,
J.G. Blom,
E.J. Spee,
G.R. Carmichael,
F.A. Potra,
Atmospheric Environment, Volume 31, 1997, Issue 20, Pages 3459-3472
\endverbatim
SourceFiles
rodas32.C
\*---------------------------------------------------------------------------*/
#ifndef rodas32_H
#define rodas32_H
#include "ODESolver.H"
#include "adaptiveSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class rodas32 Declaration
\*---------------------------------------------------------------------------*/
class rodas32
:
public ODESolver,
public adaptiveSolver
{
// Private data
mutable scalarField k1_;
mutable scalarField k2_;
mutable scalarField k3_;
mutable scalarField dy_;
mutable scalarField err_;
mutable scalarField dydx_;
mutable scalarField dfdx_;
mutable scalarSquareMatrix dfdy_;
mutable scalarSquareMatrix a_;
mutable labelList pivotIndices_;
static const scalar
c3,
d1, d2,
a31,
a41,
c21, c31, c32,
c41, c42, c43,
gamma;
public:
//- Runtime type information
TypeName("rodas32");
// Constructors
//- Construct from ODE
rodas32(const ODESystem& ode, const dictionary& dict);
// Member Functions
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y
) const;
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

3
src/ODE/ODESolvers/rodas43/rodas43.C
View File

@ -37,9 +37,6 @@ const scalar
rodas43::c2 = 0.386,
rodas43::c3 = 0.21,
rodas43::c4 = 0.63,
rodas43::bet2p = 0.0317,
rodas43::bet3p = 0.0635,
rodas43::bet4p = 0.3438,
rodas43::d1 = 0.25,
rodas43::d2 = -0.1043,
rodas43::d3 = 0.1035,

3
src/ODE/ODESolvers/rodas43/rodas43.H
View File

@ -25,7 +25,7 @@ Class
Foam::rodas43
Description
Stiffly-stable embedded Rosenbrock ODE solver of order (3)4.
L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (3)4.
References:
\verbatim
@ -79,7 +79,6 @@ class rodas43
static const scalar
c2, c3, c4,
bet2p, bet3p, bet4p,
d1, d2, d3, d4,
a21, a31, a32,
a41, a42, a43,

501
src/ODE/ODESolvers/seulex/seulex.C Normal file
View File

@ -0,0 +1,501 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 "seulex.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(seulex, 0);
addToRunTimeSelectionTable(ODESolver, seulex, dictionary);
const scalar
seulex::EPS = VSMALL,
seulex::STEPFAC1 = 0.6,
seulex::STEPFAC2 = 0.93,
seulex::STEPFAC3 = 0.1,
seulex::STEPFAC4 = 4.0,
seulex::STEPFAC5 = 0.5,
seulex::KFAC1 = 0.7,
seulex::KFAC2 = 0.9;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
:
ODESolver(ode, dict),
ode_(ode),
nSeq_(IMAXX),
cost_(IMAXX),
dfdx_(n_),
factrl_(IMAXX),
table_(KMAXX,n_),
fSave_((IMAXX-1)*(IMAXX+1)/2 + 2,n_),
dfdy_(n_),
calcJac_(false),
dense_(false),
a_(n_),
coeff_(IMAXX,IMAXX),
pivotIndices_(n_, 0.0),
hOpt_(IMAXX),
work_(IMAXX),
ySaved_(n_),
ySequence_(n_),
scale_(n_),
del_(n_),
yTemp_(n_),
dyTemp_(n_),
delTemp_(n_)
{
const scalar costfunc = 1.0, costjac = 5.0, costlu = 1.0, costsolve = 1.0;
jacRedo_ = min(1.0e-4, min(relTol_));
theta_ = 2.0*jacRedo_;
nSeq_[0] = 2;
nSeq_[1] = 3;
for (label i = 2; i < IMAXX; i++)
{
nSeq_[i] = 2*nSeq_[i-2];
}
cost_[0] = costjac + costlu + nSeq_[0]*(costfunc + costsolve);
for (label k=0;k<KMAXX;k++)
{
cost_[k+1] = cost_[k] + (nSeq_[k+1]-1)*(costfunc + costsolve) + costlu;
}
hnext_ =- VGREAT;
//NOTE: the first element of relTol_ and absTol_ are used here.
scalar logfact =- log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
kTarg_ = min(1, min(KMAXX - 1, label(logfact)));
for (label k = 0; k < IMAXX; k++)
{
for (label l = 0; l < k; l++)
{
scalar ratio = scalar(nSeq_[k])/nSeq_[l];
coeff_[k][l] = 1.0/(ratio - 1.0);
}
}
factrl_[0] = 1.0;
for (label k = 0; k < IMAXX - 1; k++)
{
factrl_[k+1] = (k+1)*factrl_[k];
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::seulex::solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
step(dxTry, x, y);
dxTry = hnext_;
}
void Foam::seulex::step(const scalar& htry, scalar& x, scalarField& y) const
{
static bool first_step = true, last_step = false;
static bool forward, reject = false, prev_reject = false;
static scalar errold;
label i, k;
scalar fac, h, hnew, err;
bool firstk;
work_[0] = 1.e30;
h = htry;
forward = h > 0 ? true : false;
ySaved_ = y;
if (h != hnext_ && !first_step)
{
last_step = true;
}
if (reject)
{
prev_reject = true;
last_step = false;
theta_=2.0*jacRedo_;
}
for (i = 0; i < n_; i++)
{
scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
}
reject = false;
firstk = true;
hnew = mag(h);
if (theta_ > jacRedo_ && !calcJac_)
{
ode_.jacobian(x, y, dfdx_, dfdy_);
calcJac_ = true;
}
while (firstk || reject)
{
h = forward ? hnew : -hnew;
firstk = false;
reject = false;
if (mag(h) <= mag(x)*EPS)
{
WarningIn("seulex::step(const scalar")
<< "step size underflow in step :" << h << endl;
}
label ipt=-1;
for (k = 0; k <= kTarg_ + 1; k++)
{
bool success = dy(x, ySaved_, h, k, ySequence_, ipt, scale_);
if (!success)
{
reject = true;
hnew = mag(h)*STEPFAC5;
break;
}
if (k == 0)
{
y = ySequence_;
}
else
{
for (i = 0; i < n_; i++)
table_[k-1][i] = ySequence_[i];
}
if (k != 0)
{
polyextr(k, table_, y);
err = 0.0;
for (i = 0; i < n_; i++)
{
scale_[i] = absTol_[i] + relTol_[i]*mag(ySaved_[i]);
err += sqr((y[i] - table_[0][i])/scale_[i]);
}
err = sqrt(err/n_);
if (err > 1.0/EPS || (k > 1 && err >= errold))
{
reject = true;
hnew = mag(h)*STEPFAC5;
break;
}
errold = min(4.0*err, 1.0);
scalar expo = 1.0/(k + 1);
scalar facmin = pow(STEPFAC3, expo);
if (err == 0.0)
{
fac = 1.0/facmin;
}
else
{
fac = STEPFAC2/pow(err/STEPFAC1, expo);
fac = max(facmin/STEPFAC4, min(1.0/facmin, fac));
}
hOpt_[k] = mag(h*fac);
work_[k] = cost_[k]/hOpt_[k];
if ((first_step || last_step) && err <= 1.0)
{
break;
}
if (k == kTarg_-1 && !prev_reject && !first_step && !last_step)
{
if (err <= 1.0)
{
break;
}
else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0)
{
reject = true;
kTarg_ = k;
if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
{
kTarg_--;
}
hnew = hOpt_[kTarg_];
break;
}
}
if (k == kTarg_)
{
if (err <= 1.0)
{
break;
}
else if (err > nSeq_[k + 1]*2.0)
{
reject = true;
if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
{
kTarg_--;
}
hnew = hOpt_[kTarg_];
break;
}
}
if (k == kTarg_+1)
{
if (err > 1.0)
{
reject = true;
if (kTarg_ > 1 && work_[kTarg_-1] < KFAC1*work_[kTarg_])
{
kTarg_--;
}
hnew = hOpt_[kTarg_];
}
break;
}
}
}
if (reject)
{
prev_reject = true;
if (!calcJac_)
{
theta_ = 2.0*jacRedo_;
if (theta_ > jacRedo_ && !calcJac_)
{
ode_.jacobian(x, y, dfdx_, dfdy_);
calcJac_ = true;
}
}
}
}
calcJac_ = false;
x += h;
hdid_ = h;
first_step = false;
label kopt;
if (k == 1)
{
kopt = 2;
}
else if (k <= kTarg_)
{
kopt=k;
if (work_[k-1] < KFAC1*work_[k])
{
kopt = k - 1;
}
else if (work_[k] < KFAC2*work_[k - 1])
{
kopt = min(k + 1, KMAXX - 1);
}
}
else
{
kopt = k - 1;
if (k > 2 && work_[k-2] < KFAC1*work_[k - 1])
{
kopt = k - 2;
}
if (work_[k] < KFAC2*work_[kopt])
{
kopt = min(k, KMAXX - 1);
}
}
if (prev_reject)
{
kTarg_ = min(kopt, k);
hnew = min(mag(h), hOpt_[kTarg_]);
prev_reject = false;
}
else
{
if (kopt <= k)
{
hnew = hOpt_[kopt];
}
else
{
if (k < kTarg_ && work_[k] < KFAC2*work_[k - 1])
{
hnew = hOpt_[k]*cost_[kopt + 1]/cost_[k];
}
else
{
hnew = hOpt_[k]*cost_[kopt]/cost_[k];
}
}
kTarg_ = kopt;
}
if (forward)
{
hnext_ = hnew;
}
else
{
hnext_ =- hnew;
}
}
bool Foam::seulex::dy
(
const scalar& x,
scalarField& y,
const scalar htot,
const label k,
scalarField& yend,
label& ipt,
scalarField& scale
) const
{
label nstep = nSeq_[k];
scalar h = htot/nstep;
for (label i = 0; i < n_; i++)
{
for (label j = 0; j < n_; j++) a_[i][j] = -dfdy_[i][j];
a_[i][i] += 1.0/h;
}
LUDecompose(a_, pivotIndices_);
scalar xnew = x + h;
ode_.derivatives(xnew, y, del_);
yTemp_ = y;
LUBacksubstitute(a_, pivotIndices_, del_);
if (dense_ && nstep == k + 1)
{
ipt++;
for (label i = 0; i < n_; i++)
{
fSave_[ipt][i] = del_[i];
}
}
for (label nn = 1; nn < nstep; nn++)
{
for (label i=0;i<n_;i++)
{
yTemp_[i] += del_[i];
}
xnew += h;
ode_.derivatives(xnew, yTemp_, yend);
if (nn == 1 && k<=1)
{
scalar del1=0.0;
for (label i = 0; i < n_; i++)
{
del1 += sqr(del_[i]/scale[i]);
}
del1 = sqrt(del1);
ode_.derivatives(x+h, yTemp_, dyTemp_);
for (label i=0;i<n_;i++)
{
del_[i] = dyTemp_[i] - del_[i]/h;
}
LUBacksubstitute(a_, pivotIndices_, del_);
scalar del2 = 0.0;
for (label i = 0; i <n_ ; i++)
{
del2 += sqr(del_[i]/scale[i]);
}
del2 = sqrt(del2);
theta_ = del2 / min(1.0, del1 + SMALL);
if (theta_ > 1.0)
{
return false;
}
}
delTemp_ = yend;
LUBacksubstitute(a_, pivotIndices_, delTemp_);
del_ = delTemp_;
if (dense_ && nn >= nstep-k-1)
{
ipt++;
for (label i=0;i<n_;i++)
{
fSave_[ipt][i]=del_[i];
}
}
}
yend = yTemp_ + del_;
return true;
}
void Foam::seulex::polyextr
(
const label k,
scalarRectangularMatrix& table,
scalarField& last
) const
{
label l=last.size();
for (label j = k - 1; j > 0; j--)
{
for (label i=0; i<l; i++)
{
table[j-1][i] =
table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
}
}
for (label i = 0; i < l; i++)
{
last[i] = table[0][i] + coeff_[k][0]*(table[0][i] - last[i]);
}
}
// ************************************************************************* //

160
src/ODE/ODESolvers/seulex/seulex.H Normal file
View File

@ -0,0 +1,160 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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::seulex
Description
An extrapolation-algorithm, based on the linearly implicit Euler method
with step size control and order selection.
The implementation is based on the SEULEX code in
\verbatim
"Solving Ordinary Differential Equations II: Stiff
and Differential-Algebraic Problems, second edition",
Hairer, E.,
Nørsett, S.,
Wanner, G.,
Springer-Verlag, Berlin. 1996.
\endverbatim
SourceFiles
seulex.C
\*---------------------------------------------------------------------------*/
#ifndef seulex_H
#define seulex_H
#include "ODESolver.H"
#include "scalarFieldField.H"
#include "scalarMatrices.H"
#include "labelField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class seulex Declaration
\*---------------------------------------------------------------------------*/
class seulex
:
public ODESolver
{
// Private data
//- seulex constants
static const scalar EPS;
static const label KMAXX=12, IMAXX = KMAXX + 1;
static const scalar
STEPFAC1, STEPFAC2, STEPFAC3, STEPFAC4, STEPFAC5, KFAC1, KFAC2;
//- Reference to ODESystem
const ODESystem& ode_;
// Temporary fields
mutable label kTarg_;
mutable labelField nSeq_;
mutable scalarField cost_, dfdx_, factrl_;
mutable scalarRectangularMatrix table_, fSave_;
mutable scalarSquareMatrix dfdy_;
mutable scalar jacRedo_, theta_;
mutable bool calcJac_, dense_;
mutable scalarSquareMatrix a_, coeff_;
mutable scalar hnext_, hdid_;
mutable labelList pivotIndices_;
// Fields space for "step" function
mutable scalarField hOpt_ ,work_, ySaved_, ySequence_, scale_;
// Fields used in "dy" function
mutable scalarField del_, yTemp_, dyTemp_, delTemp_;
// Private Member Functions
void step(const scalar& htry, scalar& x, scalarField& y) const;
bool dy
(
const scalar& x,
scalarField& y,
const scalar htot,
const label k,
scalarField& yend,
label& ipt,
scalarField& scale
) const;
void polyextr
(
const label k,
scalarRectangularMatrix& table,
scalarField& last
) const;
public:
//- Runtime type information
TypeName("seulex");
// Constructors
//- Construct from ODE
seulex(const ODESystem& ode, const dictionary& dict);
// Member Functions
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

4
src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
View File

@ -59,10 +59,10 @@ Foam::basicSymmetryPointPatchField<Type>::basicSymmetryPointPatchField
const basicSymmetryPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(p, iF)
pointPatchField<Type>(ptf, p, iF, mapper)
{}

4
src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
View File

@ -68,10 +68,10 @@ calculatedPointPatchField<Type>::calculatedPointPatchField
const calculatedPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(p, iF)
pointPatchField<Type>(ptf, p, iF, mapper)
{}

4
src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
View File

@ -61,10 +61,10 @@ coupledPointPatchField<Type>::coupledPointPatchField
const coupledPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(p, iF)
pointPatchField<Type>(ptf, p, iF, mapper)
{}

2
src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
View File

@ -108,7 +108,7 @@ Foam::valuePointPatchField<Type>::valuePointPatchField
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(p, iF),
pointPatchField<Type>(ptf, p, iF, mapper),
Field<Type>(ptf, mapper)
{}

4
src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
View File

@ -61,10 +61,10 @@ zeroGradientPointPatchField<Type>::zeroGradientPointPatchField
const zeroGradientPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(p, iF)
pointPatchField<Type>(ptf, p, iF, mapper)
{}

4
src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
View File

@ -77,10 +77,10 @@ emptyPointPatchField<Type>::emptyPointPatchField
const emptyPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(p, iF)
pointPatchField<Type>(ptf, p, iF, mapper)
{
if (!isType<emptyPointPatch>(this->patch()))
{

4
src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
View File

@ -74,10 +74,10 @@ Foam::wedgePointPatchField<Type>::wedgePointPatchField
const wedgePointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
const pointPatchFieldMapper& mapper
)
:
pointPatchField<Type>(p, iF)
pointPatchField<Type>(ptf, p, iF, mapper)
{
if (!isType<wedgePointPatch>(this->patch()))
{

16
src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C
View File

@ -63,6 +63,22 @@ pointPatchField<Type>::pointPatchField
{}
template<class Type>
Foam::pointPatchField<Type>::pointPatchField
(
const pointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
const pointPatchFieldMapper&
)
:
patch_(p),
internalField_(iF),
updated_(false),
patchType_(ptf.patchType_)
{}
template<class Type>
pointPatchField<Type>::pointPatchField
(

9
src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
View File

@ -167,6 +167,15 @@ public:
const dictionary&
);
//- Construct by mapping given patchField<Type> onto a new patch
pointPatchField
(
const pointPatchField<Type>&,
const pointPatch&,
const DimensionedField<Type, pointMesh>&,
const pointPatchFieldMapper&
);
//- Construct as copy
pointPatchField(const pointPatchField<Type>&);

4
src/OpenFOAM/matrices/solution/solution.C
View File

@ -262,7 +262,9 @@ bool Foam::solution::relaxField(const word& name) const
{
if (debug)
{
Info<< "Find variable relaxation factor for " << name << endl;
Info<< "Field relaxation factor for " << name
<< " is " << (fieldRelaxDict_.found(name) ? "set" : "unset")
<< endl;
}
return fieldRelaxDict_.found(name) || fieldRelaxDict_.found("default");

2
src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
View File

@ -79,7 +79,7 @@ inline bool Foam::pimpleControl::firstIter() const
inline bool Foam::pimpleControl::finalIter() const
{
return converged_ || (corr_ == corrPISO_);
return converged_ || (corr_ == nCorrPIMPLE_);
}

2
src/postProcessing/functionObjects/forces/forces/forces.C
View File

@ -616,7 +616,7 @@ void Foam::forces::read(const dictionary& dict)
{
initialised_ = false;
log_ = dict.lookupOrDefault<Switch>("log", true);
log_ = dict.lookupOrDefault<Switch>("log", false);
if (log_)
{

2
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
View File

@ -740,7 +740,7 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
RR_[i][celli] = (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
}
}
Info << "deltaTMin " << deltaTMin << endl;
return deltaTMin;
}

6
tutorials/combustion/chemFoam/README
View File

@ -1,13 +1,11 @@
Stiff chemistry solver validation test cases
gri : GRI-Mech 3.0. CH4 combustion, 53 species, 325 reactions
h2 : H2 combustion, 10 species, 27 reactions
nc7h16 : n-Heptane combustion, 544 species, 2446 reactions
ic8h18 : iso-Octane combustion, 874 species, 3796 reactions
Results interpreted in 'validation' sub-folder, where OpenFOAM results
are compared against those predicted by Senkin (CHEMKIN II)
are compared against those predicted by CHEMKIN II.
Overall the best performing ODE solver is seulex followed closely by rodas32.

6
tutorials/combustion/chemFoam/gri/constant/chemistryProperties
View File

@ -23,7 +23,7 @@ chemistryType
chemistry on;
initialChemicalTimeStep 1e-6;
initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
@ -33,9 +33,9 @@ EulerImplicitCoeffs
odeCoeffs
{
solver Rosenbrock43;
solver seulex;
absTol 1e-12;
relTol 1e-2;
relTol 1e-1;
}

12
tutorials/combustion/chemFoam/h2/constant/chemistryProperties
View File

@ -23,13 +23,19 @@ chemistryType
chemistry on;
initialChemicalTimeStep 1e-10;
initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs
{
solver SIBS;
solver seulex; //rodas32;
absTol 1e-12;
relTol 1e-2;
relTol 1e-1;
}

12
tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties
View File

@ -23,13 +23,19 @@ chemistryType
chemistry on;
initialChemicalTimeStep 1e-10;
initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
cTauChem 1;
equilibriumRateLimiter off;
}
odeCoeffs
{
solver SIBS;
solver seulex;
absTol 1e-12;
relTol 1e-3;
relTol 1e-1;
}

12
tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties
View File

@ -23,13 +23,19 @@ chemistryType
chemistry on;
initialChemicalTimeStep 1e-10;
initialChemicalTimeStep 1e-7;
EulerImplicitCoeffs
{
cTauChem 10;
equilibriumRateLimiter no;
}
odeCoeffs
{
solver SIBS;
solver seulex;
absTol 1e-14;
relTol 1e-4;
relTol 1e-1;
}

4
tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement
View File

@ -10,7 +10,6 @@ FoamFile
version 2.0;
format ascii;
class pointVectorField;
location "0.01";
object pointDisplacement;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -37,8 +36,9 @@ boundaryField
centreOfMass (0.5 0.5 0.5);
momentOfInertia (0.08622222 0.08622222 0.144);
mass 9.6;
rhoInf 1; // needed only for solvers solving for kinematic pressure
rhoInf 1;
report on;
accelerationRelaxation 0.3;
value uniform (0 0 0);
}
}