Skip to content

Commit

Permalink
Debugging: Complex-valued simulation types, driven with lumped ports …
Browse files Browse the repository at this point in the history 
…and no adaptive frequency sweeps and eigenmode solves

Bug fixes include: ComplexVector AXPBY and AXPBYPCX with uninitialized vector, ParOperator real-valued mass matrix orthogonalization for eigenvalue solves, unique_ptr<T[]> array construction, remove unused FEAST code (for now), remove unused PETSc code, SLEPc function pointer fixes, and others.
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Aug 8, 2023
1 parent 467edb5 commit c435ace
Show file tree
Hide file tree
Showing 22 changed files with 264 additions and 5,044 deletions.
7 changes: 5 additions & 2 deletions palace/drivers/drivensolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include "linalg/complex.hpp"
#include "linalg/ksp.hpp"
#include "linalg/operator.hpp"
#include "linalg/vector.hpp"
#include "models/lumpedportoperator.hpp"
#include "models/postoperator.hpp"
#include "models/romoperator.hpp"
Expand Down Expand Up @@ -179,7 +180,8 @@ void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, in
postop.SetEGridFunction(E);
postop.SetBGridFunction(B);
postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
Mpi::Print(" Sol. ||E|| = {:.6e} (||RHS|| = {:.6e})\n", E.Norml2(), RHS.Norml2());
Mpi::Print(" Sol. ||E|| = {:.6e} (||RHS|| = {:.6e})\n", linalg::Norml2(A->GetComm(), E),
linalg::Norml2(A->GetComm(), RHS));
if (!iodata.solver.driven.only_port_post)
{
E_elec = postop.GetEFieldEnergy();
Expand Down Expand Up @@ -337,7 +339,8 @@ void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, i
postop.SetEGridFunction(E);
postop.SetBGridFunction(B);
postop.UpdatePorts(spaceop.GetLumpedPortOp(), spaceop.GetWavePortOp(), omega);
Mpi::Print(" Sol. ||E|| = {:.6e}\n", E.Norml2());
// Mpi::Print(" Sol. ||E|| = {:.6e}\n", linalg::Norml2(A->GetComm(), E)); //XX TODO
// PROM
if (!iodata.solver.driven.only_port_post)
{
E_elec = postop.GetEFieldEnergy();
Expand Down
209 changes: 66 additions & 143 deletions palace/drivers/eigensolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
#include "linalg/arpack.hpp"
#include "linalg/complex.hpp"
#include "linalg/divfree.hpp"
#include "linalg/feast.hpp"
#include "linalg/ksp.hpp"
#include "linalg/operator.hpp"
#include "linalg/slepc.hpp"
Expand Down Expand Up @@ -46,9 +45,7 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,

// Define and configure the eigensolver to solve the eigenvalue problem:
// (K + λ C + λ² M) u = 0 or K u = -λ² M u
// with λ = iω. A shift-and-invert strategy is employed to solve for the eigenvalues
// closest to the specified target, σ. In general, the system matrices are complex and
// symmetric.
// with λ = iω. In general, the system matrices are complex and symmetric.
std::unique_ptr<EigenvalueSolver> eigen;
config::EigenSolverData::Type type = iodata.solver.eigenmode.type;
#if defined(PALACE_WITH_ARPACK) && defined(PALACE_WITH_SLEPC)
Expand Down Expand Up @@ -78,21 +75,6 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
if (type == config::EigenSolverData::Type::FEAST)
{
MFEM_ABORT("FEAST eigenvalue solver is currently not supported!");
#if defined(PALACE_WITH_SLEPC)
// Mpi::Print("\nConfiguring FEAST eigenvalue solver\n");
// if (C)
// {
// eigen = std::make_unique<feast::FeastPEPSolver>(
// K->GetComm(), iodata, spaceop, iodata.solver.eigenmode.feast_contour_np,
// iodata.problem.verbose);
// }
// else
// {
// eigen = std::make_unique<feast::FeastEPSSolver>(
// K->GetComm(), iodata, spaceop, iodata.solver.eigenmode.feast_contour_np,
// iodata.problem.verbose);
// }
#endif
}
else if (type == config::EigenSolverData::Type::ARPACK)
{
Expand Down Expand Up @@ -158,70 +140,79 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
Mpi::Print(" Scaling γ = {:.3e}, δ = {:.3e}\n", eigen->GetScalingGamma(),
eigen->GetScalingDelta());

// If desired, use an M-inner product for orthogonalizing the eigenvalue subspace. The
// constructed matrix just references the real SPD part of the mass matrix (no copy is
// performed). Boundary conditions don't need to be eliminated here.
std::unique_ptr<ParOperator> Mr;
if (iodata.solver.eigenmode.mass_orthog)
{
// Mpi::Print(" Basis uses M-inner product\n");
// Mr = std::make_unique<ParOperator>(
// std::make_unique<SumOperator>(M->LocalOperator().Real(), 1.0),
// M->GetFESpace());
// eigen->SetBMat(*Mr);

Mpi::Print(" Basis uses (K + M)-inner product\n");
auto KM = std::make_unique<SumOperator>(M->LocalOperator().Real(), 1.0);
KM->AddOperator(K->LocalOperator().Real(), 1.0);
Mr = std::make_unique<ParOperator>(std::move(KM), M->GetFESpace());
eigen->SetBMat(*Mr);
}

// Construct a divergence-free projector so the eigenvalue solve is performed in the space
// orthogonal to the zero eigenvalues of the stiffness matrix.
std::unique_ptr<DivFreeSolver> divfree;
if (iodata.solver.linear.divfree_max_it > 0)
{
constexpr int divfree_verbose = 0;
divfree = std::make_unique<DivFreeSolver>(
spaceop.GetMaterialOp(), spaceop.GetNDSpace(), spaceop.GetH1Spaces(),
spaceop.GetAuxBdrTDofLists(), iodata.solver.linear.divfree_tol,
iodata.solver.linear.divfree_max_it, divfree_verbose);
eigen->SetDivFreeProjector(*divfree);
}

// Set up the initial space for the eigenvalue solve. Satisfies boundary conditions and is
// projected appropriately.
if (iodata.solver.eigenmode.init_v0)
{
ComplexVector v0;
if (iodata.solver.eigenmode.init_v0_const)
{
Mpi::Print(" Using constant starting vector\n");
spaceop.GetConstantInitialVector(v0);
}
else
{
Mpi::Print(" Using random starting vector\n");
spaceop.GetRandomInitialVector(v0);
}
if (divfree)
{
divfree->Mult(v0);
}
eigen->SetInitialSpace(v0); // Copies the vector

// Debug
// auto Grad = spaceop.GetComplexGradMatrix();
// ComplexVector r0(Grad->Width());
// Grad->MultTranspose(v0, r0);
// r0.Print();
}

// Configure the shift-and-invert strategy is employed to solve for the eigenvalues
// closest to the specified target, σ.
const double target = iodata.solver.eigenmode.target;
const double f_target = iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, target);
std::unique_ptr<ComplexParOperator> A;
std::vector<std::unique_ptr<ParOperator>> P, AuxP;
std::unique_ptr<ComplexKspSolver> ksp;
// #if defined(PALACE_WITH_SLEPC)
// auto *feast = dynamic_cast<feast::FeastEigenSolver *>(eigen.get());
// if (feast)
// {
// // Configure the FEAST integration contour. The linear solvers are set up inside
// the
// // solver.
// if (iodata.solver.eigenmode.feast_contour_np > 1)
// {
// double contour_ub = iodata.solver.eigenmode.feast_contour_ub;
// double f_contour_ub =
// iodata.DimensionalizeValue(IoData::ValueType::FREQUENCY, contour_ub);
// double contour_ar = iodata.solver.eigenmode.feast_contour_ar;
// MFEM_VERIFY(contour_ub > target,
// "FEAST eigensolver requires a specified upper frequency target!");
// MFEM_VERIFY(
// contour_ar >= 0.0 && contour_ar <= 1.0,
// "Contour aspect ratio for FEAST eigenvalue solver must be in range
// [0.0, 1.0]!");
// Mpi::Print(" FEAST search contour: σ_lower = {:.3e} GHz ({:.3e})\n"
// " σ_upper = {:.3e} GHz ({:.3e})\n"
// " AR = {:.1e}\n",
// f_target, target, f_contour_ub, contour_ub, contour_ar);
// if (C)
// {
// // Search for eigenvalues in the range λ = iσₗₒ to iσₕᵢ.
// double h = (contour_ub - target) * contour_ar;
// feast->SetContour(-0.5 * h, target, 0.5 * h, contour_ub, false, true);
// }
// else
// {
// // Linear EVP has eigenvalues μ = -λ² = ω². Search for eigenvalues from μ =
// σₗₒ² to
// // σₕᵢ².
// double h = (contour_ub * contour_ub - target * target) * contour_ar;
// feast->SetContour(target * target, -0.5 * h, contour_ub * contour_ub, 0.5 * h);
// }
// }
// else
// {
// Mpi::Print(" FEAST search target: σ = {:.3e} GHz ({:.3e})\n", f_target, target);
// if (C)
// {
// feast->SetContour(0.0, target, 0.0, target, false, true);
// }
// else
// {
// feast->SetContour(target * target, 0.0, target * target, 0.0);
// }
// }
// }
// else
// #endif
{
Mpi::Print(" Shift-and-invert σ = {:.3e} GHz ({:.3e})\n", f_target, target);
if (C)
{
// Search for eigenvalues closest to λ = iσ.
eigen->SetShiftInvert(0.0, target);
eigen->SetShiftInvert(1i * target);
if (type == config::EigenSolverData::Type::ARPACK)
{
// ARPACK searches based on eigenvalues of the transformed problem. The eigenvalue
Expand All @@ -237,7 +228,7 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
else
{
// Linear EVP has eigenvalues μ = -λ² = ω². Search for eigenvalues closest to μ = σ².
eigen->SetShiftInvert(target * target, 0.0);
eigen->SetShiftInvert(target * target);
if (type == config::EigenSolverData::Type::ARPACK)
{
// ARPACK searches based on eigenvalues of the transformed problem. 1 / (μ - σ²)
Expand Down Expand Up @@ -265,80 +256,12 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
ksp->SetOperator(*A, P, &AuxP);
eigen->SetLinearSolver(*ksp);
}

// If desired, use an M-inner product for orthogonalizing the eigenvalue subspace. The
// constructed matrix just references the real SPD part of the mass matrix (no copy is
// performed).
std::unique_ptr<Operator> Mr;
if (iodata.solver.eigenmode.mass_orthog)
{
// Mpi::Print(" Basis uses M-inner product\n");
// Mr = std::make_unique<SumOperator>(M->Real(), 1.0);
// eigen->SetBMat(*Mr);

Mpi::Print(" Basis uses (K + M)-inner product\n");
auto KM = std::make_unique<SumOperator>(M->Real(), 1.0);
KM->AddOperator(K->Real(), 1.0);
Mr = std::move(KM);
eigen->SetBMat(*Mr);
}

// Construct a divergence-free projector so the eigenvalue solve is performed in the space
// orthogonal to the zero eigenvalues of the stiffness matrix.
std::unique_ptr<DivFreeSolver> divfree;
if (iodata.solver.linear.divfree_max_it > 0)
{
constexpr int divfree_verbose = 0;
divfree = std::make_unique<DivFreeSolver>(
spaceop.GetMaterialOp(), spaceop.GetNDSpace(), spaceop.GetH1Spaces(),
spaceop.GetAuxBdrTDofLists(), iodata.solver.linear.divfree_tol,
iodata.solver.linear.divfree_max_it, divfree_verbose);
eigen->SetDivFreeProjector(*divfree);
}

// Set up the initial space for the eigenvalue solve. Satisfies boundary conditions and is
// projected appropriately.
if (iodata.solver.eigenmode.init_v0)
{
ComplexVector v0;
if (iodata.solver.eigenmode.init_v0_const)
{
Mpi::Print(" Using constant starting vector\n");
spaceop.GetConstantInitialVector(v0);
}
else
{
Mpi::Print(" Using random starting vector\n");
spaceop.GetRandomInitialVector(v0);
}
if (divfree)
{
divfree->Mult(v0);
}
eigen->SetInitialSpace(v0); // Copies the vector

// Debug
// auto Grad = spaceop.GetComplexGradMatrix();
// ComplexVector r0(Grad->Width());
// Grad->MultTranspose(v0, r0);
// r0.Print();
}
timer.construct_time += timer.Lap();

// Eigenvalue problem solve.
Mpi::Print("\n");
int num_conv = eigen->Solve();
#if defined(PALACE_WITH_SLEPC)
// auto *feast = dynamic_cast<feast::FeastEigenSolver *>(eigen.get());
// if (feast)
// {
// SaveMetadata(feast->GetLinearSolver());
// }
// else
#endif
{
SaveMetadata(*ksp);
}
SaveMetadata(*ksp);
timer.solve_time += timer.Lap();

// Postprocess the results.
Expand Down
1 change: 0 additions & 1 deletion palace/linalg/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ target_sources(${TARGET_NAME}
${CMAKE_CURRENT_SOURCE_DIR}/curlcurl.cpp
${CMAKE_CURRENT_SOURCE_DIR}/distrelaxation.cpp
${CMAKE_CURRENT_SOURCE_DIR}/divfree.cpp
${CMAKE_CURRENT_SOURCE_DIR}/feast.cpp
${CMAKE_CURRENT_SOURCE_DIR}/gmg.cpp
${CMAKE_CURRENT_SOURCE_DIR}/jacobi.cpp
${CMAKE_CURRENT_SOURCE_DIR}/ksp.cpp
Expand Down
Loading

0 comments on commit c435ace

Please sign in to comment.