Skip to content

Commit

Permalink
Complex-valued multigrid smoothing, enabled by PR #47
Browse files Browse the repository at this point in the history 
Initial testing shows the expected 2x improvement in number of iterations, where each iteration now doubles in cost at the fine levels due to complex-valued MVP instead of real-valued. However, for expensive coarse solves, this is a substantial improvement and speedups look to be 25% on cavity example. Also adds back Chebyshev smoothing based on the standard 1st-kind polynomials, with runtime options for configuring which variant and the associated tolerances."
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Aug 23, 2023
1 parent 8ce9283 commit 13129b3
Show file tree
Hide file tree
Showing 22 changed files with 583 additions and 328 deletions.
14 changes: 7 additions & 7 deletions palace/drivers/drivensolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,11 +115,11 @@ void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, in
// Because the Dirichlet BC is always homogenous, no special elimination is required on
// the RHS. Assemble the linear system for the initial frequency (so we can call
// KspSolver::SetOperators). Compute everything at the first frequency step.
auto K = spaceop.GetComplexStiffnessMatrix(Operator::DIAG_ONE);
auto C = spaceop.GetComplexDampingMatrix(Operator::DIAG_ZERO);
auto M = spaceop.GetComplexMassMatrix(Operator::DIAG_ZERO);
auto A2 = spaceop.GetComplexExtraSystemMatrix(omega0, Operator::DIAG_ZERO);
auto Curl = spaceop.GetComplexCurlMatrix();
auto K = spaceop.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
auto C = spaceop.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
auto M = spaceop.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
auto A2 = spaceop.GetExtraSystemMatrix<ComplexOperator>(omega0, Operator::DIAG_ZERO);
auto Curl = spaceop.GetCurlMatrix<ComplexOperator>();

// Set up the linear solver and set operators for the first frequency step. The
// preconditioner for the complex linear system is constructed from a real approximation
Expand Down Expand Up @@ -154,7 +154,7 @@ void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, in
if (step > step0)
{
// Update frequency-dependent excitation and operators.
A2 = spaceop.GetComplexExtraSystemMatrix(omega, Operator::DIAG_ZERO);
A2 = spaceop.GetExtraSystemMatrix<ComplexOperator>(omega, Operator::DIAG_ZERO);
A = spaceop.GetSystemMatrix(std::complex<double>(1.0, 0.0), 1i * omega,
std::complex<double>(-omega * omega, 0.0), K.get(),
C.get(), M.get(), A2.get());
Expand Down Expand Up @@ -236,7 +236,7 @@ void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, i

// Allocate negative curl matrix for postprocessing the B-field and vectors for the
// high-dimensional field solution.
auto Curl = spaceop.GetComplexCurlMatrix();
auto Curl = spaceop.GetCurlMatrix<ComplexOperator>();
ComplexVector E(Curl->Width()), B(Curl->Height());
E = 0.0;
B = 0.0;
Expand Down
10 changes: 5 additions & 5 deletions palace/drivers/eigensolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,10 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
// computational range. The damping matrix may be nullptr.
timer.Lap();
SpaceOperator spaceop(iodata, mesh);
auto K = spaceop.GetComplexStiffnessMatrix(Operator::DIAG_ONE);
auto C = spaceop.GetComplexDampingMatrix(Operator::DIAG_ZERO);
auto M = spaceop.GetComplexMassMatrix(Operator::DIAG_ZERO);
auto Curl = spaceop.GetComplexCurlMatrix();
auto K = spaceop.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
auto C = spaceop.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
auto M = spaceop.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
auto Curl = spaceop.GetCurlMatrix<ComplexOperator>();
SaveMetadata(spaceop.GetNDSpace());

// Configure objects for postprocessing.
Expand Down Expand Up @@ -190,7 +190,7 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
eigen->SetInitialSpace(v0); // Copies the vector

// Debug
// auto Grad = spaceop.GetComplexGradMatrix();
// auto Grad = spaceop.GetGradMatrix<ComplexOperator>();
// ComplexVector r0(Grad->Width());
// Grad->MultTranspose(v0, r0);
// r0.Print();
Expand Down
225 changes: 141 additions & 84 deletions palace/linalg/chebyshev.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@

#include "chebyshev.hpp"

#include <vector>
#include <mfem/general/forall.hpp>
#include "linalg/rap.hpp"

Expand All @@ -20,27 +19,22 @@ void GetInverseDiagonal(const ParOperator &A, Vector &dinv)
dinv.Reciprocal();
}

void GetInverseDiagonal(const ComplexParOperator &A, Vector &dinv)
void GetInverseDiagonal(const ComplexParOperator &A, ComplexVector &dinv)
{
MFEM_VERIFY(A.HasReal() && !A.HasImag(),
"ComplexOperator for ChebyshevSmoother must be real-valued for now!");
MFEM_VERIFY(A.HasReal() || A.HasImag(),
"Invalid zero ComplexOperator for ChebyshevSmoother!");
dinv.SetSize(A.Height());
A.Real()->AssembleDiagonal(dinv);
dinv.SetSize(A.Height());
dinv = 0.0;
if (A.HasReal())
{
A.Real()->AssembleDiagonal(dinv.Real());
}
if (A.HasImag())
{
A.Imag()->AssembleDiagonal(dinv.Imag());
}
dinv.Reciprocal();
// MFEM_VERIFY(A.HasReal() || A.HasImag(),
// "Invalid zero ComplexOperator for ChebyshevSmoother!");
// dinv.SetSize(A.Height());
// dinv.SetSize(A.Height());
// dinv = 0.0;
// if (A.HasReal())
// {
// A.Real()->AssembleDiagonal(dinv.Real());
// }
// if (A.HasImag())
// {
// A.Imag()->AssembleDiagonal(dinv.Imag());
// }
// dinv.Reciprocal();
}

double GetLambdaMax(MPI_Comm comm, const Operator &A, const Vector &dinv)
Expand All @@ -50,48 +44,13 @@ double GetLambdaMax(MPI_Comm comm, const Operator &A, const Vector &dinv)
return linalg::SpectralNorm(comm, DinvA, false);
}

double GetLambdaMax(MPI_Comm comm, const ComplexOperator &A, const Vector &dinv)
double GetLambdaMax(MPI_Comm comm, const ComplexOperator &A, const ComplexVector &dinv)
{
MFEM_VERIFY(A.HasReal() && !A.HasImag(),
"ComplexOperator for ChebyshevSmoother must be real-valued for now!");
DiagonalOperator Dinv(dinv);
ProductOperator DinvA(Dinv, *A.Real());
ComplexDiagonalOperator Dinv(dinv);
ComplexProductOperator DinvA(Dinv, A);
return linalg::SpectralNorm(comm, DinvA, false);
}

} // namespace

template <typename OperType>
ChebyshevSmoother<OperType>::ChebyshevSmoother(int smooth_it, int poly_order)
: Solver<OperType>(), pc_it(smooth_it), order(poly_order), A(nullptr)
{
}

template <typename OperType>
void ChebyshevSmoother<OperType>::SetOperator(const OperType &op)
{
using ParOperType =
typename std::conditional<std::is_same<OperType, ComplexOperator>::value,
ComplexParOperator, ParOperator>::type;

A = &op;
r.SetSize(op.Height());
d.SetSize(op.Height());

const auto *PtAP = dynamic_cast<const ParOperType *>(&op);
MFEM_VERIFY(PtAP,
"ChebyshevSmoother requires a ParOperator or ComplexParOperator operator!");
GetInverseDiagonal(*PtAP, dinv);

// Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. See
// mfem::OperatorChebyshevSmoother or Adams et al., Parallel multigrid smoothing:
// polynomial versus Gauss-Seidel, JCP (2003).
lambda_max = 1.01 * GetLambdaMax(PtAP->GetComm(), *A, dinv);
}

namespace
{

template <bool Transpose = false>
inline void ApplyOp(const Operator &A, const Vector &x, Vector &y)
{
Expand Down Expand Up @@ -137,42 +96,37 @@ inline void ApplyOrder0(double sr, const Vector &dinv, const Vector &r, Vector &
const int N = d.Size();
const auto *DI = dinv.Read();
const auto *R = r.Read();
auto *D = d.ReadWrite();
auto *D = d.Write();
mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { D[i] = sr * DI[i] * R[i]; });
}

template <bool Transpose = false>
inline void ApplyOrder0(const double sr, const Vector &dinv, const ComplexVector &r,
inline void ApplyOrder0(const double sr, const ComplexVector &dinv, const ComplexVector &r,
ComplexVector &d)
{
const int N = dinv.Size();
// const auto *DIR = dinv.Real().Read();
// const auto *DII = dinv.Imag().Read();
const auto *DIR = dinv.Read();
const auto *DIR = dinv.Real().Read();
const auto *DII = dinv.Imag().Read();
const auto *RR = r.Real().Read();
const auto *RI = r.Imag().Read();
auto *DR = d.Real().ReadWrite();
auto *DI = d.Imag().ReadWrite();
auto *DR = d.Real().Write();
auto *DI = d.Imag().Write();
if constexpr (!Transpose)
{
mfem::forall(N,
[=] MFEM_HOST_DEVICE(int i)
{
// DR[i] = sr * (DIR[i] * RR[i] - DII[i] * RI[i]);
// DI[i] = sr * (DII[i] * RR[i] + DIR[i] * RI[i]);
DR[i] = sr * DIR[i] * RR[i];
DI[i] = sr * DIR[i] * RI[i];
DR[i] = sr * (DIR[i] * RR[i] - DII[i] * RI[i]);
DI[i] = sr * (DII[i] * RR[i] + DIR[i] * RI[i]);
});
}
else
{
mfem::forall(N,
[=] MFEM_HOST_DEVICE(int i)
{
// DR[i] = sr * (DIR[i] * RR[i] + DII[i] * RI[i]);
// DI[i] = sr * (-DII[i] * RR[i] + DIR[i] * RI[i]);
DR[i] = sr * DIR[i] * RR[i];
DI[i] = sr * DIR[i] * RI[i];
DR[i] = sr * (DIR[i] * RR[i] + DII[i] * RI[i]);
DI[i] = sr * (-DII[i] * RR[i] + DIR[i] * RI[i]);
});
}
}
Expand All @@ -189,13 +143,12 @@ inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv,
}

template <bool Transpose = false>
inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv,
inline void ApplyOrderK(const double sd, const double sr, const ComplexVector &dinv,
const ComplexVector &r, ComplexVector &d)
{
const int N = dinv.Size();
// const auto *DIR = dinv.Real().Read();
// const auto *DII = dinv.Imag().Read();
const auto *DIR = dinv.Read();
const auto *DIR = dinv.Real().Read();
const auto *DII = dinv.Imag().Read();
const auto *RR = r.Real().Read();
const auto *RI = r.Imag().Read();
auto *DR = d.Real().ReadWrite();
Expand All @@ -205,27 +158,54 @@ inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv,
mfem::forall(N,
[=] MFEM_HOST_DEVICE(int i)
{
// DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] - DII[i] * RI[i]);
// DI[i] = sd * DI[i] + sr * (DII[i] * RR[i] + DIR[i] * RI[i]);
DR[i] = sd * DR[i] + sr * DIR[i] * RR[i];
DI[i] = sd * DI[i] + sr * DIR[i] * RI[i];
DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] - DII[i] * RI[i]);
DI[i] = sd * DI[i] + sr * (DII[i] * RR[i] + DIR[i] * RI[i]);
});
}
else
{
mfem::forall(N,
[=] MFEM_HOST_DEVICE(int i)
{
// DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] + DII[i] * RI[i]);
// DI[i] = sd * DI[i] + sr * (-DII[i] * RR[i] + DIR[i] * RI[i]);
DR[i] = sd * DR[i] + sr * DIR[i] * RR[i];
DI[i] = sd * DI[i] + sr * DIR[i] * RI[i];
DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] + DII[i] * RI[i]);
DI[i] = sd * DI[i] + sr * (-DII[i] * RR[i] + DIR[i] * RI[i]);
});
}
}

} // namespace

template <typename OperType>
ChebyshevSmoother<OperType>::ChebyshevSmoother(int smooth_it, int poly_order, double sf_max)
: Solver<OperType>(), pc_it(smooth_it), order(poly_order), A(nullptr), lambda_max(0.0),
sf_max(sf_max)
{
MFEM_VERIFY(order > 0, "Polynomial order for Chebyshev smoothing must be positive!");
}

template <typename OperType>
void ChebyshevSmoother<OperType>::SetOperator(const OperType &op)
{
using ParOperType =
typename std::conditional<std::is_same<OperType, ComplexOperator>::value,
ComplexParOperator, ParOperator>::type;

A = &op;
r.SetSize(op.Height());
d.SetSize(op.Height());
this->height = op.Height();
this->width = op.Width();

const auto *PtAP = dynamic_cast<const ParOperType *>(&op);
MFEM_VERIFY(PtAP,
"ChebyshevSmoother requires a ParOperator or ComplexParOperator operator!");
GetInverseDiagonal(*PtAP, dinv);

// Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. See
// mfem::OperatorChebyshevSmoother or Adams et al. (2003).
lambda_max = sf_max * GetLambdaMax(PtAP->GetComm(), *A, dinv);
}

template <typename OperType>
void ChebyshevSmoother<OperType>::Mult(const VecType &x, VecType &y) const
{
Expand Down Expand Up @@ -258,7 +238,84 @@ void ChebyshevSmoother<OperType>::Mult(const VecType &x, VecType &y) const
}
}

template <typename OperType>
ChebyshevSmoother1stKind<OperType>::ChebyshevSmoother1stKind(int smooth_it, int poly_order,
double sf_max, double sf_min)
: Solver<OperType>(), pc_it(smooth_it), order(poly_order), A(nullptr), theta(0.0),
sf_max(sf_max), sf_min(sf_min)
{
MFEM_VERIFY(order > 0, "Polynomial order for Chebyshev smoothing must be positive!");
}

template <typename OperType>
void ChebyshevSmoother1stKind<OperType>::SetOperator(const OperType &op)
{
using ParOperType =
typename std::conditional<std::is_same<OperType, ComplexOperator>::value,
ComplexParOperator, ParOperator>::type;

A = &op;
r.SetSize(op.Height());
d.SetSize(op.Height());
this->height = op.Height();
this->width = op.Width();

const auto *PtAP = dynamic_cast<const ParOperType *>(&op);
MFEM_VERIFY(
PtAP,
"ChebyshevSmoother1stKind requires a ParOperator or ComplexParOperator operator!");
GetInverseDiagonal(*PtAP, dinv);

// Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. The
// optimized estimate of lambda_min comes from (2.24) of Phillips and Fischer (2022).
if (sf_min <= 0.0)
{
sf_min = 1.69 / (std::pow(order, 1.68) + 2.11 * order + 1.98);
}
const double lambda_max = sf_max * GetLambdaMax(PtAP->GetComm(), *A, dinv);
const double lambda_min = sf_min * lambda_max;
theta = 0.5 * (lambda_max + lambda_min);
delta = 0.5 * (lambda_max - lambda_min);
}

template <typename OperType>
void ChebyshevSmoother1stKind<OperType>::Mult(const VecType &x, VecType &y) const
{
// Apply smoother: y = y + p(A) (x - A y) .
for (int it = 0; it < pc_it; it++)
{
if (this->initial_guess || it > 0)
{
ApplyOp(*A, y, r);
linalg::AXPBY(1.0, x, -1.0, r);
}
else
{
r = x;
y = 0.0;
}

// 1th-kind Chebyshev smoother, from Phillips and Fischer or Adams.
ApplyOrder0(1.0 / theta, dinv, r, d);
double rhop = delta / theta;
for (int k = 1; k < order; k++)
{
y += d;
ApplyOp(*A, d, r, -1.0);
const double rho = 1.0 / (2.0 * theta / delta - rhop);
const double sd = rho * rhop;
const double sr = 2.0 * rho / delta;
ApplyOrderK(sd, sr, dinv, r, d);
rhop = rho;
}
y += d;
}
}

template class ChebyshevSmoother<Operator>;
template class ChebyshevSmoother<ComplexOperator>;

template class ChebyshevSmoother1stKind<Operator>;
template class ChebyshevSmoother1stKind<ComplexOperator>;

} // namespace palace
Loading

0 comments on commit 13129b3

Please sign in to comment.