Skip to content

Commit

Permalink
Merged feature/p-UAMG into master
Browse files Browse the repository at this point in the history 
Former-commit-id: ef4bafe7800f6c1f6d69afa181de8b897a3d01b9
  • Loading branch information
@pmatalon
pmatalon committed Nov 29, 2021
2 parents e2046b2 + 59ea011 commit 5831458
Show file tree
Hide file tree
Showing 3 changed files with 154 additions and 19 deletions.
3 changes: 2 additions & 1 deletion src/Program.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -464,8 +464,9 @@ class ProgramDim : public Program
else if (args.Solver.SolverCode.compare("uamg") == 0)
{
Diffusion_HHO<Dim>* hhoProblem = dynamic_cast<Diffusion_HHO<Dim>*>(problem);
UncondensedAMG* mg = new UncondensedAMG(hhoProblem->HHO->nCellUnknowns, hhoProblem->HHO->nFaceUnknowns, 0.25, args.Solver.MG.UAMGFaceProlong, args.Solver.MG.UAMGCoarseningProlong, args.Solver.MG.UAMGMultigridProlong, args.Solver.MG.Levels);
UncondensedAMG* mg = new UncondensedAMG(Dim, hhoProblem->HHO->FaceBasis->GetDegree(), hhoProblem->HHO->nCellUnknowns, hhoProblem->HHO->nFaceUnknowns, 0.25, args.Solver.MG.UAMGFaceProlong, args.Solver.MG.UAMGCoarseningProlong, args.Solver.MG.UAMGMultigridProlong, args.Solver.MG.Levels);
SetMultigridParameters(mg, args, blockSize);
mg->CoarsePolyDegree = 0;
solver = mg;
}
else if (args.Solver.SolverCode.compare("aggregamg") == 0)
Expand Down
29 changes: 22 additions & 7 deletions src/Solver/Multigrid/UncondensedAMG/UncondensedAMG.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,18 @@ class UncondensedAMG : public Multigrid
UAMGFaceProlongation _faceProlong = UAMGFaceProlongation::FaceAggregates;
UAMGProlongation _coarseningProlong = UAMGProlongation::FaceProlongation;
UAMGProlongation _multigridProlong = UAMGProlongation::ReconstructSmoothedTraceOrInject;
int _dim;
int _degree;
int _cellBlockSize;
int _faceBlockSize;
double _strongCouplingThreshold;
public:

UncondensedAMG(int cellBlockSize, int faceBlockSize, double strongCouplingThreshold, UAMGFaceProlongation faceProlong, UAMGProlongation coarseningProlong, UAMGProlongation mgProlong, int nLevels = 0)
UncondensedAMG(int dim, int degree, int cellBlockSize, int faceBlockSize, double strongCouplingThreshold, UAMGFaceProlongation faceProlong, UAMGProlongation coarseningProlong, UAMGProlongation mgProlong, int nLevels = 0)
: Multigrid(nLevels)
{
this->_dim = dim;
this->_degree = degree;
this->_cellBlockSize = cellBlockSize;
this->_faceBlockSize = faceBlockSize;
this->_strongCouplingThreshold = strongCouplingThreshold;
Expand Down Expand Up @@ -103,16 +107,27 @@ class UncondensedAMG : public Multigrid
protected:
Level* CreateFineLevel() const override
{
return new UncondensedLevel(0, _cellBlockSize, _faceBlockSize, _strongCouplingThreshold, _faceProlong, _coarseningProlong, _multigridProlong);
return new UncondensedLevel(0, _degree, _cellBlockSize, _faceBlockSize, _strongCouplingThreshold, _faceProlong, _coarseningProlong, _multigridProlong);
}

Level* CreateCoarseLevel(Level* fineLevel, CoarseningType coarseningType, int coarseDegree) override
{
if (coarseningType != CoarseningType::H)
Utils::FatalError("Only h-coarsening allowed for this multigrid.");
if (coarseningType == CoarseningType::HP)
Utils::FatalError("hp-coarsening not allowed for this multigrid.");

UncondensedLevel* coarseLevel = new UncondensedLevel(fineLevel->Number + 1, _cellBlockSize, _faceBlockSize, _strongCouplingThreshold, _faceProlong, _coarseningProlong, _multigridProlong);
coarseLevel->OperatorMatrix = &dynamic_cast<UncondensedLevel*>(fineLevel)->Ac;
return coarseLevel;
UncondensedLevel* fine = dynamic_cast<UncondensedLevel*>(fineLevel);

if (coarseningType == CoarseningType::P)
{
int coarseCellBlockSize = Utils::Binomial(coarseDegree + _dim , coarseDegree);
int coarseFaceBlockSize = Utils::Binomial(coarseDegree + _dim - 1, coarseDegree);
UncondensedLevel* coarseLevel = new UncondensedLevel(fine->Number + 1, coarseDegree, coarseCellBlockSize, coarseFaceBlockSize, _strongCouplingThreshold, _faceProlong, _coarseningProlong, _multigridProlong);
return coarseLevel;
}
else
{
UncondensedLevel* coarseLevel = new UncondensedLevel(fine->Number + 1, fine->PolynomialDegree(), fine->CellBlockSize(), fine->FaceBlockSize(), _strongCouplingThreshold, _faceProlong, _coarseningProlong, _multigridProlong);
return coarseLevel;
}
}
};
141 changes: 130 additions & 11 deletions src/Solver/Multigrid/UncondensedAMG/UncondensedLevel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ class UncondensedLevel : public Level
UAMGFaceProlongation _faceProlong = UAMGFaceProlongation::FaceAggregates;
UAMGProlongation _coarseningProlong = UAMGProlongation::FaceProlongation;
UAMGProlongation _multigridProlong = UAMGProlongation::ReconstructSmoothedTraceOrInject;
int _degree;
int _cellBlockSize;
int _faceBlockSize;
double _strongCouplingThreshold;
Expand All @@ -29,12 +30,13 @@ class UncondensedLevel : public Level
SparseMatrix Ac;

public:
UncondensedLevel(int number, int cellBlockSize, int faceBlockSize, double strongCouplingThreshold, UAMGFaceProlongation faceProlong, UAMGProlongation coarseningProlong, UAMGProlongation mgProlong)
UncondensedLevel(int number, int degree, int cellBlockSize, int faceBlockSize, double strongCouplingThreshold, UAMGFaceProlongation faceProlong, UAMGProlongation coarseningProlong, UAMGProlongation mgProlong)
: Level(number)
{
if (cellBlockSize > 1 || faceBlockSize > 1)
Utils::Warning("This multigrid is efficient if cellBlockSize = faceBlockSize = 1. It may converge badly.");
//if (cellBlockSize > 1 || faceBlockSize > 1)
//Utils::Warning("This multigrid is efficient if cellBlockSize = faceBlockSize = 1. It may converge badly.");

this->_degree = degree;
this->_cellBlockSize = cellBlockSize;
this->_faceBlockSize = faceBlockSize;
this->_strongCouplingThreshold = strongCouplingThreshold;
Expand All @@ -50,7 +52,22 @@ class UncondensedLevel : public Level

int PolynomialDegree() override
{
return 0;
return _degree;
}

int BlockSizeForBlockSmoothers() override
{
return _faceBlockSize;
}

int CellBlockSize()
{
return _cellBlockSize;
}

int FaceBlockSize()
{
return _faceBlockSize;
}

void ExportVector(const Vector& v, string suffix, int levelNumber) override
Expand Down Expand Up @@ -628,7 +645,43 @@ class UncondensedLevel : public Level
this->OperatorMatrix = new SparseMatrix(fine->Q_F.transpose() * *(fine->OperatorMatrix) * fine->Q_F);
}

public:
void SetupOperatorByBlockExtraction()
{
UncondensedLevel* fine = dynamic_cast<UncondensedLevel*>(FinerLevel);

auto nElems = fine->A_T_F->rows() / fine->_cellBlockSize;
auto nFaces = fine->A_T_F->cols() / fine->_faceBlockSize;

this->OperatorMatrix = ExtractCoarseMatrix(*this->FinerLevel->OperatorMatrix, nFaces, nFaces, fine->_faceBlockSize, fine->_faceBlockSize, this->_faceBlockSize, this->_faceBlockSize);
this->A_T_T = ExtractCoarseMatrix(*fine->A_T_T, nElems, nElems, fine->_cellBlockSize, fine->_cellBlockSize, this->_cellBlockSize, this->_cellBlockSize);
this->A_T_F = ExtractCoarseMatrix(*fine->A_T_F, nElems, nFaces, fine->_cellBlockSize, fine->_faceBlockSize, this->_cellBlockSize, this->_faceBlockSize);
this->A_F_F = ExtractCoarseMatrix(*fine->A_F_F, nFaces, nFaces, fine->_faceBlockSize, fine->_faceBlockSize, this->_faceBlockSize, this->_faceBlockSize);
}

private:
static SparseMatrix* ExtractCoarseMatrix(const SparseMatrix& fineMatrix, BigNumber nBlockRows, BigNumber nBlockCols, int fineBlockRows, int fineBlockCols, int coarseBlockRows, int coarseBlockCols)
{
NumberParallelLoop<CoeffsChunk> parallelLoop(nBlockRows);
parallelLoop.Execute([&fineMatrix, fineBlockRows, fineBlockCols, coarseBlockRows, coarseBlockCols](BigNumber i, ParallelChunk<CoeffsChunk>* chunk)
{
for (int k = 0; k < coarseBlockRows; k++)
{
for (RowMajorSparseMatrix::InnerIterator it(fineMatrix, i*fineBlockRows + k); it; ++it)
{
auto j = it.col() / fineBlockCols;
int l = it.col() - j * fineBlockCols;
if (l < coarseBlockCols)
chunk->Results.Coeffs.Add(i*coarseBlockRows + k, j*coarseBlockCols + l, it.value());
}
}
});
SparseMatrix* coarseMatrix = new SparseMatrix(nBlockRows * coarseBlockRows, nBlockCols * coarseBlockCols);
parallelLoop.Fill(*coarseMatrix);
return coarseMatrix;
}


// Cell prolongation Q_T with only one 1 coefficient per row
SparseMatrix BuildQ_T(const HybridAlgebraicMesh& mesh)
{
Expand Down Expand Up @@ -868,6 +921,7 @@ class UncondensedLevel : public Level
public:
void OnStartSetup() override
{
cout << "\t\tk = " << this->PolynomialDegree() << endl;
cout << "\t\tMesh : " << this->A_T_T->rows() / _cellBlockSize << " elements, " << this->A_T_F->cols() / _faceBlockSize << " faces";
if (!this->IsFinestLevel())
{
Expand All @@ -884,20 +938,85 @@ class UncondensedLevel : public Level

void SetupRestriction() override
{
double scalingFactor = 1.0;
//scalingFactor = 1.0 / 4.0;
R = (scalingFactor * P.transpose()).eval();
//R = scalingFactor * P.transpose();
if (this->CoarserLevel->ComesFrom == CoarseningType::P)
{
// nothing to do
}
else
{
double scalingFactor = 1.0;
//scalingFactor = 1.0 / 4.0;
R = (scalingFactor * P.transpose()).eval();
//R = scalingFactor * P.transpose();
}
}

Vector Prolong(Vector& vectorOnTheCoarserLevel) override
{
if (this->CoarserLevel->ComesFrom == CoarseningType::P)
{
// This is possible only if the basis is hierarchical
UncondensedLevel* coarseLevel = dynamic_cast<UncondensedLevel*>(CoarserLevel);
auto nHigherDegreeUnknowns = this->_faceBlockSize;
auto nLowerDegreeUnknowns = coarseLevel->_faceBlockSize;
auto nFaces = vectorOnTheCoarserLevel.rows() / nLowerDegreeUnknowns;
Vector vectorOnThisLevel = Vector::Zero(nFaces * nHigherDegreeUnknowns);
for (BigNumber i = 0; i < nFaces; i++)
vectorOnThisLevel.segment(i*nHigherDegreeUnknowns, nLowerDegreeUnknowns) = vectorOnTheCoarserLevel.segment(i*nLowerDegreeUnknowns, nLowerDegreeUnknowns);
return vectorOnThisLevel;
}
else
return Level::Prolong(vectorOnTheCoarserLevel);
}

Vector Restrict(Vector& vectorOnThisLevel) override
{
if (this->CoarserLevel->ComesFrom == CoarseningType::P)
{
// This makes sense only if the basis is hierarchical and orthogonal
UncondensedLevel* coarseLevel = dynamic_cast<UncondensedLevel*>(CoarserLevel);
auto nHigherDegreeUnknowns = this->_faceBlockSize;
auto nLowerDegreeUnknowns = coarseLevel->_faceBlockSize;
auto nFaces = vectorOnThisLevel.rows() / nHigherDegreeUnknowns;
Vector vectorOnTheCoarserLevel(nFaces * nLowerDegreeUnknowns);
for (BigNumber i = 0; i < nFaces; i++)
vectorOnTheCoarserLevel.segment(i*nLowerDegreeUnknowns, nLowerDegreeUnknowns) = vectorOnThisLevel.segment(i*nHigherDegreeUnknowns, nLowerDegreeUnknowns);
return vectorOnTheCoarserLevel;
}
else
return Level::Restrict(vectorOnThisLevel);
}

Flops ProlongCost() override
{
if (this->CoarserLevel->ComesFrom == CoarseningType::P)
return 0;
else
return Level::ProlongCost();
}

Flops RestrictCost() override
{
if (this->CoarserLevel->ComesFrom == CoarseningType::P)
return 0;
else
return Level::RestrictCost();
}

void OnEndSetup() override
{
if (this->CoarserLevel)
{
UncondensedLevel* coarse = dynamic_cast<UncondensedLevel*>(this->CoarserLevel);
coarse->A_T_T = &A_T_Tc;
coarse->A_T_F = &A_T_Fc;
coarse->A_F_F = &A_F_Fc;
if (this->CoarserLevel->ComesFrom == CoarseningType::P)
coarse->SetupOperatorByBlockExtraction();
else
{
coarse->OperatorMatrix = &Ac;
coarse->A_T_T = &A_T_Tc;
coarse->A_T_F = &A_T_Fc;
coarse->A_F_F = &A_F_Fc;
}
}
}

Expand Down

0 comments on commit 5831458

Please sign in to comment.