38 lduMatrix::solver::addsymMatrixConstructorToTable<PBiCGStab>
41 lduMatrix::solver::addasymMatrixConstructorToTable<PBiCGStab>
48 Foam::PBiCGStab::PBiCGStab
50 const word& fieldName,
86 const label nCells =
psi.
size();
88 solveScalar* __restrict__ psiPtr =
psi.
begin();
91 solveScalar* __restrict__ pAPtr = pA.begin();
94 solveScalar* __restrict__ yAPtr = yA.begin();
97 matrix_.Amul(yA,
psi, interfaceBouCoeffs_, interfaces_, cmpt);
101 solveScalar* __restrict__ rAPtr = rA.begin();
111 const solveScalar normFactor = this->normFactor(
psi, source, yA, pA);
115 Info<<
" Normalisation factor = " << normFactor <<
endl;
119 solverPerf.initialResidual() =
122 solverPerf.finalResidual() = solverPerf.initialResidual();
128 || !solverPerf.checkConvergence(tolerance_, relTol_, log_)
132 solveScalar* __restrict__ AyAPtr = AyA.begin();
135 solveScalar* __restrict__ sAPtr = sA.begin();
138 solveScalar* __restrict__ zAPtr = zA.begin();
141 solveScalar* __restrict__ tAPtr = tA.begin();
147 solveScalar rA0rA = 0;
148 solveScalar
alpha = 0;
149 solveScalar omega = 0;
165 const solveScalar rA0rAold = rA0rA;
170 if (solverPerf.checkSingularity(
mag(rA0rA)))
176 if (solverPerf.nIterations() == 0)
178 for (label cell=0; cell<nCells; cell++)
180 pAPtr[cell] = rAPtr[cell];
186 if (solverPerf.checkSingularity(
mag(omega)))
191 const solveScalar
beta = (rA0rA/rA0rAold)*(
alpha/omega);
193 for (label cell=0; cell<nCells; cell++)
196 rAPtr[cell] +
beta*(pAPtr[cell] - omega*AyAPtr[cell]);
201 preconPtr_->precondition(yA, pA, cmpt);
204 matrix_.Amul(AyA, yA, interfaceBouCoeffs_, interfaces_, cmpt);
206 const solveScalar rA0AyA =
209 alpha = rA0rA/rA0AyA;
212 for (label cell=0; cell<nCells; cell++)
214 sAPtr[cell] = rAPtr[cell] -
alpha*AyAPtr[cell];
218 solverPerf.finalResidual() =
223 solverPerf.nIterations() >= minIter_
224 && solverPerf.checkConvergence(tolerance_, relTol_, log_)
227 for (label cell=0; cell<nCells; cell++)
229 psiPtr[cell] +=
alpha*yAPtr[cell];
232 solverPerf.nIterations()++;
238 preconPtr_->precondition(zA, sA, cmpt);
241 matrix_.Amul(tA, zA, interfaceBouCoeffs_, interfaces_, cmpt);
243 const solveScalar tAtA =
gSumSqr(tA, matrix().
mesh().comm());
250 for (label cell=0; cell<nCells; cell++)
252 psiPtr[cell] +=
alpha*yAPtr[cell] + omega*zAPtr[cell];
253 rAPtr[cell] = sAPtr[cell] - omega*tAPtr[cell];
256 solverPerf.finalResidual() =
262 ++solverPerf.nIterations() < maxIter_
263 && !solverPerf.checkConvergence(tolerance_, relTol_, log_)
265 || solverPerf.nIterations() < minIter_
271 preconPtr_->setFinished(solverPerf);
276 ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
void size(const label n)
Older name for setAddressableSize.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Field< solveScalar > solveScalarField
A list of keyword definitions, which are a keyword followed by a number of values (eg...
void setResidualField(const scalarField &residual, const word &fieldName, const bool initial) const
Set the residual field using an IOField on the object registry if it exists.
scalarProduct< Type, Type >::type gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
Ostream & endl(Ostream &os)
Add newline and flush stream.
A non-const Field/List wrapper with possible data conversion.
lduMatrix::solver::addsymMatrixConstructorToTable< PBiCGStab > addPBiCGStabSymMatrixConstructorToTable_
A field of fields is a PtrList of fields with reference counting.
A class for handling words, derived from Foam::string.
static autoPtr< preconditioner > New(const solver &sol, const dictionary &solverControls)
Return a new preconditioner.
A const Field/List wrapper with possible data conversion.
lduMatrix::solver::addasymMatrixConstructorToTable< PBiCGStab > addPBiCGStabAsymMatrixConstructorToTable_
iterator begin() noexcept
Return an iterator to begin traversing the UList.
virtual solverPerformance solve(scalarField &psi, const scalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.
int debug
Static debugging option.
Container< Type > & ref() const
Return non-const reference to the contents of a non-null managed pointer.
defineTypeNameAndDebug(combustionModel, 0)
outerProduct1< Type >::type gSumSqr(const UList< Type > &f, const label comm)
typeOfMag< Type >::type gSumMag(const FieldField< Field, Type > &f)
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
static word getName(const dictionary &)
Find the preconditioner name (directly or from a sub-dictionary)
const volScalarField & psi
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Preconditioned bi-conjugate gradient stabilized solver for asymmetric lduMatrices using a run-time se...
virtual solverPerformance scalarSolve(solveScalarField &psi, const solveScalarField &source, const direction cmpt=0) const
Solve the matrix with this solver.