ENH: GAMGSolver: initial working version

This commit is contained in:
mattijs
2013-04-03 16:46:19 +01:00
parent f966d8a584
commit 8acd81694e
5 changed files with 943 additions and 532 deletions

View File

@ -142,8 +142,8 @@ void Foam::GAMGSolver::Vcycle
const label coarsestLevel = matrixLevels_.size() - 1;
// Restrict finest grid residual for the next level up
agglomeration_.restrictField(coarseSources[0], finestResidual, 0);
// Restrict finest grid residual for the next level up.
agglomeration_.restrictField(coarseSources[0], finestResidual, 0, true);
if (debug >= 2 && nPreSweeps_)
{
@ -154,66 +154,81 @@ void Foam::GAMGSolver::Vcycle
// Residual restriction (going to coarser levels)
for (label leveli = 0; leveli < coarsestLevel; leveli++)
{
// If the optional pre-smoothing sweeps are selected
// smooth the coarse-grid field for the restriced source
if (nPreSweeps_)
Pout<< "Restriction for level:" << leveli << endl;
if (coarseSources.set(leveli + 1))
{
coarseCorrFields[leveli] = 0.0;
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
coarseSources[leveli],
cmpt,
min
(
nPreSweeps_ + preSweepsLevelMultiplier_*leveli,
maxPreSweeps_
)
);
scalarField::subField ACf
(
Apsi,
coarseCorrFields[leveli].size()
);
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
if (scaleCorrection_ && leveli < coarsestLevel - 1)
// If the optional pre-smoothing sweeps are selected
// smooth the coarse-grid field for the restriced source
if (nPreSweeps_)
{
scale
coarseCorrFields[leveli] = 0.0;
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
const_cast<scalarField&>(ACf.operator const scalarField&()),
matrixLevels_[leveli],
coarseSources[leveli],
cmpt,
min
(
nPreSweeps_ + preSweepsLevelMultiplier_*leveli,
maxPreSweeps_
)
);
scalarField::subField ACf
(
Apsi,
coarseCorrFields[leveli].size()
);
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
if (scaleCorrection_ && leveli < coarsestLevel - 1)
{
int comm = matrixLevels_[leveli].mesh().comm();
scale
(
coarseCorrFields[leveli],
const_cast<scalarField&>
(
ACf.operator const scalarField&()
),
matrixLevels_[leveli],
comm,
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
}
// Correct the residual with the new solution
matrixLevels_[leveli].Amul
(
const_cast<scalarField&>
(
ACf.operator const scalarField&()
),
coarseCorrFields[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
coarseSources[leveli] -= ACf;
}
// Correct the residual with the new solution
matrixLevels_[leveli].Amul
// Residual is equal to source
agglomeration_.restrictField
(
const_cast<scalarField&>(ACf.operator const scalarField&()),
coarseCorrFields[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
cmpt
coarseSources[leveli + 1],
coarseSources[leveli],
leveli + 1,
true
);
coarseSources[leveli] -= ACf;
}
// Residual is equal to source
agglomeration_.restrictField
(
coarseSources[leveli + 1],
coarseSources[leveli],
leveli + 1
);
}
if (debug >= 2 && nPreSweeps_)
@ -223,12 +238,15 @@ void Foam::GAMGSolver::Vcycle
// Solve Coarsest level with either an iterative or direct solver
solveCoarsestLevel
(
coarseCorrFields[coarsestLevel],
coarseSources[coarsestLevel]
);
if (coarseCorrFields.set(coarsestLevel))
{
Pout<< "Coarsest solve for level:" << coarsestLevel << endl;
solveCoarsestLevel
(
coarseCorrFields[coarsestLevel],
coarseSources[coarsestLevel]
);
}
if (debug >= 2)
{
@ -237,99 +255,141 @@ void Foam::GAMGSolver::Vcycle
// Smoothing and prolongation of the coarse correction fields
// (going to finer levels)
scalarField dummyField(0);
for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
{
// Create a field for the pre-smoothed correction field
// as a sub-field of the finestCorrection which is not
// currently being used
scalarField::subField preSmoothedCoarseCorrField
(
finestCorrection,
coarseCorrFields[leveli].size()
);
Pout<< "Smoothing and prolongation for level:" << leveli << endl;
// Only store the preSmoothedCoarseCorrField if pre-smoothing is used
if (nPreSweeps_)
if (coarseCorrFields.set(leveli))
{
preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
}
Pout<< "Prolonging from " << leveli + 1 << " up to "
<< leveli << endl;
//// Create a field for the pre-smoothed correction field
//// as a sub-field of the finestCorrection which is not
//// currently being used
//scalarField::subField preSmoothedCoarseCorrField
//(
// finestCorrection,
// coarseCorrFields[leveli].size()
//);
scalarField preSmoothedCoarseCorrField;
agglomeration_.prolongField
(
coarseCorrFields[leveli],
coarseCorrFields[leveli + 1],
leveli + 1
);
// Only store the preSmoothedCoarseCorrField if pre-smoothing is
// used
if (nPreSweeps_)
{
//preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
preSmoothedCoarseCorrField = coarseCorrFields[leveli];
}
// Create A.psi for this coarse level as a sub-field of Apsi
scalarField::subField ACf
(
Apsi,
coarseCorrFields[leveli].size()
);
scalarField& ACfRef =
const_cast<scalarField&>(ACf.operator const scalarField&());
if (interpolateCorrection_)
{
interpolate
agglomeration_.prolongField
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
(
coarseCorrFields.set(leveli + 1)
? coarseCorrFields[leveli + 1]
: dummyField // dummy value
),
leveli + 1,
true
);
}
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
if (scaleCorrection_ && leveli < coarsestLevel - 1)
{
scale
Pout<< "Doing stuff at level " << leveli << endl;
//// Create A.psi for this coarse level as a sub-field of Apsi
//scalarField::subField ACf
//(
// Apsi,
// coarseCorrFields[leveli].size()
//);
//scalarField& ACfRef =
// const_cast<scalarField&>(ACf.operator const scalarField&());
scalarField ACfRef(coarseCorrFields[leveli].size());
if (interpolateCorrection_)
{
Pout<< "doing interpolate." << endl;
interpolate
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
Pout<< "done interpolate." << endl;
}
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
if (scaleCorrection_ && leveli < coarsestLevel - 1)
{
//int comm =
//(
// matrixLevels_.set(leveli+1)
// ? matrixLevels_[leveli+1].mesh().comm()
// : matrixLevels_[leveli].mesh().comm()
//);
int comm = matrixLevels_[leveli].mesh().comm();
Pout<< "doing scale with comm:" << comm << endl;
scale
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
comm,
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
Pout<< "done scale with comm:" << comm << endl;
}
// Only add the preSmoothedCoarseCorrField if pre-smoothing is
// used
if (nPreSweeps_)
{
coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
}
Pout<< "doing smooth." << endl;
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
cmpt,
min
(
nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
maxPostSweeps_
)
);
Pout<< "done smooth." << endl;
}
// Only add the preSmoothedCoarseCorrField if pre-smoothing is used
if (nPreSweeps_)
{
coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
}
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
coarseSources[leveli],
cmpt,
min
(
nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
maxPostSweeps_
)
);
}
// Prolong the finest level correction
Pout<< "Doing Prolong to finest level" << endl;
agglomeration_.prolongField
(
finestCorrection,
coarseCorrFields[0],
0
0,
true //false // no proc agglomeration for now
);
Pout<< "Done Prolong to finest level" << endl;
if (interpolateCorrection_)
{
Pout<< "doing interpolate on finest level" << endl;
interpolate
(
finestCorrection,
@ -340,21 +400,26 @@ void Foam::GAMGSolver::Vcycle
finestResidual,
cmpt
);
Pout<< "done interpolate on finest level" << endl;
}
if (scaleCorrection_)
{
// Scale the finest level correction
int comm = matrix_.mesh().comm();
Pout<< "doing scale on finest level with comm:" << comm << endl;
scale
(
finestCorrection,
Apsi,
matrix_,
comm,
interfaceBouCoeffs_,
interfaces_,
finestResidual,
cmpt
);
Pout<< "done scale on finest level with comm:" << comm << endl;
}
forAll(psi, i)
@ -362,6 +427,7 @@ void Foam::GAMGSolver::Vcycle
psi[i] += finestCorrection[i];
}
Pout<< "Doing smooth on finest level" << endl;
smoothers[0].smooth
(
psi,
@ -369,6 +435,7 @@ void Foam::GAMGSolver::Vcycle
cmpt,
nFinestSweeps_
);
Pout<< "Done smooth on finest level" << endl;
}
@ -400,37 +467,47 @@ void Foam::GAMGSolver::initVcycle
forAll(matrixLevels_, leveli)
{
coarseCorrFields.set
(
leveli,
new scalarField
(
agglomeration_.meshLevel(leveli + 1).lduAddr().size()
)
);
if (agglomeration_.nCells(leveli) >= 0)
{
label nCoarseCells = agglomeration_.nCells(leveli);
coarseSources.set
(
leveli,
new scalarField
(
agglomeration_.meshLevel(leveli + 1).lduAddr().size()
)
);
Pout<< "initVCucle level:" << leveli << " nCoarseCells:"
<< nCoarseCells << endl;
smoothers.set
(
leveli + 1,
lduMatrix::smoother::New
coarseSources.set(leveli, new scalarField(nCoarseCells));
//if (!matrixLevels_.set(leveli))
//{
// coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
//}
}
if (matrixLevels_.set(leveli))
{
const lduMatrix& mat = matrixLevels_[leveli];
label nCoarseCells = mat.diag().size();
Pout<< "initVCucle level:" << leveli << " matrix size:"
<< nCoarseCells << endl;
coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
//coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
//coarseSources.set(leveli, new scalarField(nCoarseCells));
smoothers.set
(
fieldName_,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevelsIntCoeffs_[leveli],
interfaceLevels_[leveli],
controlDict_
)
);
leveli + 1,
lduMatrix::smoother::New
(
fieldName_,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevelsIntCoeffs_[leveli],
interfaceLevels_[leveli],
controlDict_
)
);
}
}
}
@ -456,106 +533,111 @@ void Foam::GAMGSolver::solveCoarsestLevel
coarsestCorrField = coarsestSource;
coarsestLUMatrixPtr_->solve(coarsestCorrField);
}
else if (processorAgglomerate_)
{
const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
(
coarsestLevel
);
scalarField allSource;
globalIndex cellOffsets;
if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
{
cellOffsets.offsets() = agglomeration_.cellOffsets(coarsestLevel);
}
cellOffsets.gather
(
coarseComm,
agglomProcIDs,
coarsestSource,
allSource
);
scalarField allCorrField;
solverPerformance coarseSolverPerf;
label solveComm = agglomeration_.procCommunicator(coarsestLevel);
label oldWarn = UPstream::warnComm;
UPstream::warnComm = solveComm;
if (Pstream::myProcNo(solveComm) != -1)
{
const lduMatrix& allMatrix = allMatrixPtr_();
allCorrField.setSize(allSource.size(), 0);
{
Pout<< "** Master:Solving on comm:" << solveComm
<< " with procs:" << UPstream::procID(solveComm) << endl;
if (allMatrix.asymmetric())
{
coarseSolverPerf = BICCG
(
"coarsestLevelCorr",
allMatrix,
allInterfaceBouCoeffs_,
allInterfaceIntCoeffs_,
allInterfaces_,
tolerance_,
relTol_
).solve
(
allCorrField,
allSource
);
}
else
{
coarseSolverPerf = ICCG
(
"coarsestLevelCorr",
allMatrix,
allInterfaceBouCoeffs_,
allInterfaceIntCoeffs_,
allInterfaces_,
tolerance_,
relTol_
).solve
(
allCorrField,
allSource
);
}
}
}
UPstream::warnComm = oldWarn;
Pout<< "done master solve." << endl;
// Scatter to all processors
coarsestCorrField.setSize(coarsestSource.size());
cellOffsets.scatter
(
coarseComm,
agglomProcIDs,
allCorrField,
coarsestCorrField
);
if (debug >= 2)
{
coarseSolverPerf.print(Info(coarseComm));
}
Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl;
Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
}
//else if
//(
// agglomeration_.processorAgglomerate()
// && procMatrixLevels_.set(coarsestLevel)
//)
//{
// //const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
// //(
// // coarsestLevel
// //);
// //
// //scalarField allSource;
// //
// //globalIndex cellOffsets;
// //if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
// //{
// // cellOffsets.offsets() =
// // agglomeration_.cellOffsets(coarsestLevel);
// //}
// //
// //cellOffsets.gather
// //(
// // coarseComm,
// // agglomProcIDs,
// // coarsestSource,
// // allSource
// //);
// //
// //scalarField allCorrField;
// //solverPerformance coarseSolverPerf;
//
// label solveComm = agglomeration_.procCommunicator(coarsestLevel);
// label oldWarn = UPstream::warnComm;
// UPstream::warnComm = solveComm;
//
//
// coarsestCorrField = 0;
// solverPerformance coarseSolverPerf;
//
// if (Pstream::myProcNo(solveComm) != -1)
// {
// const lduMatrix& allMatrix = procMatrixLevels_[coarsestLevel];
//
// {
// Pout<< "** Master:Solving on comm:" << solveComm
// << " with procs:" << UPstream::procID(solveComm) << endl;
//
// if (allMatrix.asymmetric())
// {
// coarseSolverPerf = BICCG
// (
// "coarsestLevelCorr",
// allMatrix,
// procInterfaceLevelsBouCoeffs_[coarsestLevel],
// procInterfaceLevelsIntCoeffs_[coarsestLevel],
// procInterfaceLevels_[coarsestLevel],
// tolerance_,
// relTol_
// ).solve
// (
// coarsestCorrField,
// coarsestSource
// );
// }
// else
// {
// coarseSolverPerf = ICCG
// (
// "coarsestLevelCorr",
// allMatrix,
// procInterfaceLevelsBouCoeffs_[coarsestLevel],
// procInterfaceLevelsIntCoeffs_[coarsestLevel],
// procInterfaceLevels_[coarsestLevel],
// tolerance_,
// relTol_
// ).solve
// (
// coarsestCorrField,
// coarsestSource
// );
// }
// }
// }
//
// UPstream::warnComm = oldWarn;
// Pout<< "done master solve." << endl;
//
// //// Scatter to all processors
// //coarsestCorrField.setSize(coarsestSource.size());
// //cellOffsets.scatter
// //(
// // coarseComm,
// // agglomProcIDs,
// // allCorrField,
// // coarsestCorrField
// //);
//
// if (debug >= 2)
// {
// coarseSolverPerf.print(Info(coarseComm));
// }
//
// Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl;
// Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
//}
else
{
coarsestCorrField = 0;