Skip to content

Commit

Permalink
BiCGStab solver + generate polygonal mesh from cartesian
Browse files Browse the repository at this point in the history
  • Loading branch information
@pmatalon
pmatalon committed Aug 21, 2023
1 parent 807e3e4 commit 21d06de
Show file tree
Hide file tree
Showing 9 changed files with 189 additions and 49 deletions.
10 changes: 9 additions & 1 deletion src/Mesh/PolyhedralMesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -812,7 +812,15 @@ class PolyhedralMesh : public Mesh<Dim>

if (faceCoarseningStgy == FaceCoarseningStrategy::None)
{
// do nothing
if (bdryFaceCollapsing != FaceCollapsing::Disabled)
{
ElementParallelLoop<Dim> parallelLoopCollapseFaces(coarseMesh->Elements);
parallelLoopCollapseFaces.Execute([coarseMesh, bdryFaceCollapsing](Element<Dim>* coarseElement)
{
if (!coarseElement->IsDeleted && coarseElement->IsOnBoundary())
coarseMesh->TryCollapseBoundaryFaces(coarseElement, bdryFaceCollapsing);
});
}
}
else if (faceCoarseningStgy == FaceCoarseningStrategy::InterfaceCollapsing)
{
Expand Down
10 changes: 8 additions & 2 deletions src/Mesher/MeshFactory.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,8 +198,14 @@ Mesh<2>* MeshFactory<2>::BuildMesh(ProgramArguments& args, TestCase<2>* testCase
}
else if (meshCode.compare("poly") == 0)
{
Square_GMSHUnstructTriangularMesh* triMesh = new Square_GMSHUnstructTriangularMesh(n);
fineMesh = BuildPolyhedralMesh(triMesh, args.Discretization.PolyMeshFaceCoarseningStgy, args.Discretization.PolyMeshBoundaryFaceCollapsing, args.Discretization.PolyMeshNAggregPasses);
PolyhedralMesh<2>* initMesh;
if (args.Discretization.PolyMeshInitialMesh.compare("cart") == 0)
initMesh = new Square_GMSHCartesianMesh(n);
else if (args.Discretization.PolyMeshInitialMesh.compare("stri") == 0)
initMesh = new Square_GMSHTriangularMesh(n);
else if (args.Discretization.PolyMeshInitialMesh.compare("tri") == 0)
initMesh = new Square_GMSHUnstructTriangularMesh(n);
fineMesh = BuildPolyhedralMesh(initMesh, args.Discretization.PolyMeshFaceCoarseningStgy, args.Discretization.PolyMeshBoundaryFaceCollapsing, args.Discretization.PolyMeshNAggregPasses);
}
else
Utils::FatalError("The requested mesh is not managed with this geometry.");
Expand Down
23 changes: 18 additions & 5 deletions src/Program/Program_BiHarmonic_HHO.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ class Program_BiHarmonic_HHO

Timer lapSolverSetupTimer;

cout << "Solver: " << *diffSolver << endl << endl;
cout << "Solver: " << *diffSolver << endl;
IterativeSolver* iterativeSolver = dynamic_cast<IterativeSolver*>(diffSolver);
if (iterativeSolver)
{
Expand All @@ -179,6 +179,7 @@ class Program_BiHarmonic_HHO
diffSolver->Setup(biHarPb->DiffPb().A);
lapSolverSetupTimer.Stop();
}
cout << "Setup time: CPU = " << lapSolverSetupTimer.CPU() << ", elapsed = " << lapSolverSetupTimer.Elapsed() << endl << endl;

biHarPb->SetDiffSolver(diffSolver);

Expand Down Expand Up @@ -445,18 +446,27 @@ class Program_BiHarmonic_HHO
if (args.Problem.Scheme.compare("g") == 0)
{
BiHarmonicMixedFormGlowinski_HHO<Dim>* gloScheme = static_cast<BiHarmonicMixedFormGlowinski_HHO<Dim>*>(biHarPb);

ProgramArguments precSolverArgs;
precSolverArgs.Solver.SolverCode = args.Solver.BiHarmonicPrecSolverCode;
precSolverArgs.Solver.Tolerance = args.Solver.Tolerance;
if (args.Solver.BiHarmonicPrecSolverTol > 0)
precSolverArgs.Solver.Tolerance = args.Solver.BiHarmonicPrecSolverTol;
precSolverArgs.Solver.MaxIterations = args.Solver.BiHarmonicPrecSolverMaxIter;
Solver* precSolver = SolverFactory<Dim>::CreateSolver(precSolverArgs, blockSizeForBlockSolver, out);

if (args.Solver.BiHarmonicPreconditionerCode.compare("j") == 0)
cg->Precond = new DenseBlockJacobiPreconditioner(1);
else if (args.Solver.BiHarmonicPreconditionerCode.compare("bj") == 0)
cg->Precond = new DenseBlockJacobiPreconditioner(hho->nFaceUnknowns);
else if (args.Solver.BiHarmonicPreconditionerCode.compare("s") == 0)
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::SingleFaceNeighbourhood, 0, args.Solver.NeighbourhoodDepth, false);
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::SingleFaceNeighbourhood, 0, args.Solver.NeighbourhoodDepth, precSolver, false);
else if (args.Solver.BiHarmonicPreconditionerCode.compare("ds") == 0)
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::SingleFaceNeighbourhood, 0, args.Solver.NeighbourhoodDepth, true);
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::SingleFaceNeighbourhood, 0, args.Solver.NeighbourhoodDepth, precSolver, true);
else if (args.Solver.BiHarmonicPreconditionerCode.compare("p") == 0)
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::FacePatchNeighbourhood, args.Solver.PatchSize, args.Solver.NeighbourhoodDepth, false);
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::FacePatchNeighbourhood, args.Solver.PatchSize, args.Solver.NeighbourhoodDepth, precSolver, false);
else if (args.Solver.BiHarmonicPreconditionerCode.compare("dp") == 0)
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::FacePatchNeighbourhood, args.Solver.PatchSize, args.Solver.NeighbourhoodDepth, true);
cg->Precond = new BiharPatchPreconditioner<Dim>(*gloScheme, BiharPatchPreconditioner<Dim>::Type::FacePatchNeighbourhood, args.Solver.PatchSize, args.Solver.NeighbourhoodDepth, precSolver, true);
else if (args.Solver.BiHarmonicPreconditionerCode.compare("no") == 0)
cg->Precond = new IdentityPreconditioner();
else
Expand Down Expand Up @@ -616,6 +626,9 @@ class Program_BiHarmonic_HHO
double error = biHarPb->DiffPb().BoundarySpace.RelativeL2Error(theta_f, discreteExact);
cout << "L2 Error (theta) = " << std::scientific << error << endl;
}

if (testCase->ExactSolution)
cout << "h = " << mesh->H() << endl;
}
}

Expand Down
5 changes: 5 additions & 0 deletions src/ProgramArguments.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ struct DiscretizationArguments
string Mesher = "default";
string MeshCode = "default";
double Stretch = 0.5;
string PolyMeshInitialMesh = "cart";
FaceCoarseningStrategy PolyMeshFaceCoarseningStgy = FaceCoarseningStrategy::InterfaceCollapsing;
FaceCollapsing PolyMeshBoundaryFaceCollapsing = FaceCollapsing::OnlyCollinear;
int PolyMeshNAggregPasses = 1;
Expand Down Expand Up @@ -95,8 +96,12 @@ struct SolverArguments
double RelaxationParameter = 1;
int BlockSize = -1;
int Restart = 0;

string BiHarmonicSolverCode = "fcg";
string BiHarmonicPreconditionerCode = "s";
string BiHarmonicPrecSolverCode = "bicgstab";
double BiHarmonicPrecSolverTol = 0;
double BiHarmonicPrecSolverMaxIter = 1000;
bool ComputeIterL2Error = false;
MultigridArguments MG;
bool BiHarReconstructBoundary = false;
Expand Down
63 changes: 39 additions & 24 deletions src/Solver/BiHarmonic/BiharPatchPreconditioner.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,16 +21,19 @@ class BiharPatchPreconditioner : public Preconditioner
int _facePatchSize = 3;
int _neighbourhoodDepth = 1;
bool _blockDiagPrec = false;
SparseMatrix _precondMatrix;
vector<Eigen::FullPivLU<DenseMatrix>> _invD;
EigenSparseLU _solver;
Solver* _solver;

public:
BiharPatchPreconditioner(BiHarmonicMixedFormGlowinski_HHO<Dim>& biHarPb, Type type, int facePatchSize, int neighbourhoodDepth, bool blockDiagPrec = false) :
BiharPatchPreconditioner(BiHarmonicMixedFormGlowinski_HHO<Dim>& biHarPb, Type type, int facePatchSize, int neighbourhoodDepth, Solver* solver, bool blockDiagPrec = false) :
_diffPb(biHarPb.DiffPb())
{
_type = type;
_facePatchSize = facePatchSize;
_neighbourhoodDepth = neighbourhoodDepth;
_blockDiagPrec = blockDiagPrec;
_solver = solver;
}

void Serialize(ostream& os) const override
Expand All @@ -47,18 +50,28 @@ class BiharPatchPreconditioner : public Preconditioner
os << "diagonal ";
os << "face patch neighbourhood (patch size = " << _facePatchSize << ", neighbourhood depth = " << _neighbourhoodDepth << ")";
}
if (!_blockDiagPrec)
{
os << std::endl;
os << " solver: " << (*_solver);
IterativeSolver* iterSolver = dynamic_cast<IterativeSolver*>(_solver);
if (iterSolver)
{
os << ", tol=" << iterSolver->Tolerance << ", maxIter=" << iterSolver->MaxIterations;
}
}
}

void Setup(const DenseMatrix& A) override
{
if (_type == Type::SingleFaceNeighbourhood)
Setup1(A);
Setup_SingleFaceNeighbourhood(A);
else if (_type == Type::FacePatchNeighbourhood)
Setup2(A);
Setup_FacePatchNeighbourhood(A);
}

private:
void Setup1(const DenseMatrix& A)
void Setup_SingleFaceNeighbourhood(const DenseMatrix& A)
{
ExportModule out(Utils::ProgramArgs.OutputDirectory, "", Utils::ProgramArgs.Actions.Export.ValueSeparator);

Expand Down Expand Up @@ -113,37 +126,38 @@ class BiharPatchPreconditioner : public Preconditioner
}
});

SparseMatrix mat(_diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns, _diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns);
parallelLoop.Fill(mat);
//SparseMatrix mat(_diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns, _diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns);
_precondMatrix = SparseMatrix(_diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns, _diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns);
parallelLoop.Fill(_precondMatrix);

//SparseMatrix matT = mat.triangularView<Eigen::StrictlyLower>().transpose();
//mat = mat.triangularView<Eigen::Lower>() + matT;

if (Utils::ProgramArgs.Actions.PrintDebug)
{
if (Dim <= 2)
cout << "Preconditioner matrix:" << endl << DenseMatrix(mat) << endl << endl;
cout << "Preconditioner matrix:" << endl << DenseMatrix(_precondMatrix) << endl << endl;
else
cout << "Preconditioner matrix:" << endl << DenseMatrix(mat.topLeftCorner(3 * _diffPb.HHO->nFaceUnknowns, 3 * _diffPb.HHO->nFaceUnknowns)) << endl << endl;
cout << "||A-mat|| = " << (A - DenseMatrix(mat)).norm() << endl;
cout << "Preconditioner matrix:" << endl << DenseMatrix(_precondMatrix.topLeftCorner(3 * _diffPb.HHO->nFaceUnknowns, 3 * _diffPb.HHO->nFaceUnknowns)) << endl << endl;
cout << "||A-mat|| = " << (A - DenseMatrix(_precondMatrix)).norm() << endl;
}

if (_blockDiagPrec)
{
int blockSize = _diffPb.HHO->nFaceUnknowns;
_invD = vector<Eigen::FullPivLU<DenseMatrix>>(_diffPb.HHO->nBoundaryFaces);
NumberParallelLoop<EmptyResultChunk> parallelLoop2(_diffPb.HHO->nBoundaryFaces);
parallelLoop2.Execute([this, &mat, blockSize](BigNumber i, ParallelChunk<EmptyResultChunk>* chunk)
parallelLoop2.Execute([this, blockSize](BigNumber i, ParallelChunk<EmptyResultChunk>* chunk)
{
DenseMatrix Di = mat.block(i * blockSize, i * blockSize, blockSize, blockSize);
DenseMatrix Di = _precondMatrix.block(i * blockSize, i * blockSize, blockSize, blockSize);
_invD[i].compute(Di);
});
}
else
_solver.Setup(mat);
_solver->Setup(_precondMatrix);
}

void Setup2(const DenseMatrix& A)
void Setup_FacePatchNeighbourhood(const DenseMatrix& A)
{
vector<BoundaryFacePatch<Dim>> patches;
patches.reserve(_diffPb._mesh->BoundaryFaces.size() / _facePatchSize);
Expand Down Expand Up @@ -249,17 +263,18 @@ class BiharPatchPreconditioner : public Preconditioner
}
});

SparseMatrix mat(_diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns, _diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns);
parallelLoop.Fill(mat);
//SparseMatrix mat(_diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns, _diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns);
_precondMatrix = SparseMatrix(_diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns, _diffPb.HHO->nBoundaryFaces * _diffPb.HHO->nFaceUnknowns);
parallelLoop.Fill(_precondMatrix);

if (Utils::ProgramArgs.Actions.PrintDebug)
{
if (Dim <= 2)
cout << "Preconditioner matrix:" << endl << DenseMatrix(mat) << endl << endl;
cout << "Preconditioner matrix:" << endl << DenseMatrix(_precondMatrix) << endl << endl;
else
cout << "Preconditioner matrix:" << endl << DenseMatrix(mat.topLeftCorner(3 * _diffPb.HHO->nFaceUnknowns, 3 * _diffPb.HHO->nFaceUnknowns)) << endl << endl;
cout << "Condition number = " << Utils::Cond(DenseMatrix(mat)) << endl;
cout << "||A-mat|| = " << (A - DenseMatrix(mat)).norm() << endl;
cout << "Preconditioner matrix:" << endl << DenseMatrix(_precondMatrix.topLeftCorner(3 * _diffPb.HHO->nFaceUnknowns, 3 * _diffPb.HHO->nFaceUnknowns)) << endl << endl;
cout << "Condition number = " << Utils::Cond(DenseMatrix(_precondMatrix)) << endl;
cout << "||A-mat|| = " << (A - DenseMatrix(_precondMatrix)).norm() << endl;
//cout << "A-mat: " << endl << (A - DenseMatrix(mat)) << endl << endl;
}

Expand All @@ -268,14 +283,14 @@ class BiharPatchPreconditioner : public Preconditioner
int blockSize = _diffPb.HHO->nFaceUnknowns;
_invD = vector<Eigen::FullPivLU<DenseMatrix>>(_diffPb.HHO->nBoundaryFaces);
NumberParallelLoop<EmptyResultChunk> parallelLoop2(_diffPb.HHO->nBoundaryFaces);
parallelLoop2.Execute([this, &mat, blockSize](BigNumber i, ParallelChunk<EmptyResultChunk>* chunk)
parallelLoop2.Execute([this, blockSize](BigNumber i, ParallelChunk<EmptyResultChunk>* chunk)
{
DenseMatrix Di = mat.block(i * blockSize, i * blockSize, blockSize, blockSize);
DenseMatrix Di = _precondMatrix.block(i * blockSize, i * blockSize, blockSize, blockSize);
_invD[i].compute(Di);
});
}
else
_solver.Setup(mat);
_solver->Setup(_precondMatrix);
}

Vector BiharOperator(const Vector& dirichlet, NeighbourhoodDiffusion_HHO<Dim>& nbhDiff, const SparseMatrix& S_iF_bF_transpose, const SparseMatrix& Theta_T_bF_transpose,
Expand Down Expand Up @@ -354,7 +369,7 @@ class BiharPatchPreconditioner : public Preconditioner
return res;
}
else
return _solver.Solve(r);
return _solver->Solve(r);
}

MFlops SolvingComputationalWork() override
Expand Down
2 changes: 1 addition & 1 deletion src/Solver/Direct/EigenSparseLU.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ class EigenSparseLU : public Solver
void Setup(const SparseMatrix& A) override
{
Solver::Setup(A);
_solver.isSymmetric(true);
//_solver.isSymmetric(true);
_solver.compute(A);
this->SetupComputationalWork = Cost::LUFactorization(A)*1e-6;
Eigen::ComputationInfo info = _solver.info();
Expand Down
48 changes: 48 additions & 0 deletions src/Solver/Krylov/EigenBiCGSTAB.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
#pragma once
#include "../IterativeSolver.h"
#include <Eigen/IterativeLinearSolvers>
using namespace std;

class EigenBiCGSTAB : public IterativeSolver
{
private:
Eigen::BiCGSTAB<SparseMatrix, Eigen::DiagonalPreconditioner<double>> _solver;

public:

EigenBiCGSTAB() : IterativeSolver() {}

void Serialize(ostream& os) const override
{
os << "BiCGSTAB (Eigen library)";
}

void Setup(const SparseMatrix& A) override
{
Solver::Setup(A);
_solver.compute(A);
Eigen::ComputationInfo info = _solver.info();
if (info != Eigen::ComputationInfo::Success)
{
cout << "----------------- A -------------------" << A << endl;
cout << "Error: ConjugateGradient failed to execute with the code " << info << "." << endl;
exit(EXIT_FAILURE);
}
_solver.setTolerance(this->Tolerance);
_solver.setMaxIterations(this->MaxIterations);
}

void Solve(const Vector& b, Vector& x, bool xEquals0, bool computeResidual, bool computeAx) override
{
if (xEquals0)
x = _solver.solve(b);
else
x = _solver.solveWithGuess(b, x);
this->IterationCount = _solver.iterations();
}

Vector Solve(const Vector& b) override
{
return IterativeSolver::Solve(b);
}
};
3 changes: 3 additions & 0 deletions src/Solver/SolverFactory.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "Multigrid/AggregAMG/AggregAMG.h"
#include "FixedPoint/BlockJacobi.h"
#include "Krylov/EigenCG.h"
#include "Krylov/EigenBiCGSTAB.h"
#include "Multigrid/AggregAMG/HighOrderAggregAMG.h"
#include "Multigrid/AGMG.h"

Expand Down Expand Up @@ -85,6 +86,8 @@ class SolverFactory
solver = new ConjugateGradient();
else if (args.Solver.SolverCode.compare("eigencg") == 0)
solver = new EigenCG();
else if (args.Solver.SolverCode.compare("bicgstab") == 0)
solver = new EigenBiCGSTAB();
else if (args.Solver.SolverCode.compare("agmg") == 0)
solver = new AGMG();
else if (args.Solver.SolverCode.compare("aggregamg") == 0)
Expand Down

0 comments on commit 21d06de

Please sign in to comment.