diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp index 39c7f0471..013704977 100644 --- a/palace/drivers/eigensolver.cpp +++ b/palace/drivers/eigensolver.cpp @@ -148,13 +148,13 @@ EigenSolver::Solve(const std::vector> &mesh) const std::unique_ptr KM; if (iodata.solver.eigenmode.mass_orthog) { - // Mpi::Print(" Basis uses M-inner product\n"); - // KM = spaceop.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get()); - // eigen->SetBMat(*KM); - - Mpi::Print(" Basis uses (K + M)-inner product\n"); - KM = spaceop.GetInnerProductMatrix(1.0, 1.0, K.get(), M.get()); + Mpi::Print(" Basis uses M-inner product\n"); + KM = spaceop.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get()); eigen->SetBMat(*KM); + + // Mpi::Print(" Basis uses (K + M)-inner product\n"); + // KM = spaceop.GetInnerProductMatrix(1.0, 1.0, K.get(), M.get()); + // eigen->SetBMat(*KM); } // Construct a divergence-free projector so the eigenvalue solve is performed in the space @@ -266,20 +266,21 @@ EigenSolver::Solve(const std::vector> &mesh) const } SaveMetadata(*ksp); - // Calculate and record the error indicators. - Mpi::Print("Computing solution error estimates\n\n"); + // Calculate and record the error indicators and postprocess the results. + BlockTimer bt2(Timer::POSTPRO); + Mpi::Print("Computing solution error estimates and postprocessing computed modes\n\n"); CurlFluxErrorEstimator estimator( spaceop.GetMaterialOp(), spaceop.GetNDSpace(), iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0); ErrorIndicator indicator; - for (int i = 0; i < iodata.solver.eigenmode.n; i++) + if (!KM) { - eigen->GetEigenvector(i, E); - estimator.AddErrorIndicator(E, indicator); + // Normalize the finalized eigenvectors with respect to mass matrix (unit electric field + // energy) even if they are not computed to be orthogonal with respect to it. + KM = spaceop.GetInnerProductMatrix(0.0, 1.0, nullptr, M.get()); + eigen->SetBMat(*KM); + eigen->RescaleEigenvectors(num_conv); } - - // Postprocess the results. - BlockTimer bt2(Timer::POSTPRO); for (int i = 0; i < num_conv; i++) { // Get the eigenvalue and relative error. @@ -300,6 +301,11 @@ EigenSolver::Solve(const std::vector> &mesh) const // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in // PostOperator for all postprocessing operations. eigen->GetEigenvector(i, E); + if (i < iodata.solver.eigenmode.n) + { + // Only consider the desired number of modes for the error indicator. + estimator.AddErrorIndicator(E, indicator); + } Curl.Mult(E.Real(), B.Real()); Curl.Mult(E.Imag(), B.Imag()); B *= -1.0 / (1i * omega); diff --git a/palace/linalg/arpack.cpp b/palace/linalg/arpack.cpp index 96a3bcd58..94e42f2fc 100644 --- a/palace/linalg/arpack.cpp +++ b/palace/linalg/arpack.cpp @@ -225,6 +225,7 @@ void ArpackEigenvalueSolver::SetNumModes(int num_eig, int num_vec) eig.reset(); perm.reset(); res.reset(); + xscale.reset(); } if (ncv > 0 && num_vec != ncv) { @@ -435,6 +436,23 @@ void ArpackEigenvalueSolver::GetEigenvector(int i, ComplexVector &x) const MFEM_VERIFY(x.Size() == n, "Invalid size mismatch for provided eigenvector!"); const int &j = perm.get()[i]; x.Set(V.get() + j * n, n); + if (xscale.get()[j] > 0.0) + { + x *= xscale.get()[j]; + } +} + +double ArpackEigenvalueSolver::GetEigenvectorNorm(const ComplexVector &x, + ComplexVector &Bx) const +{ + if (opB) + { + return linalg::Norml2(comm, x, *opB, Bx); + } + else + { + return linalg::Norml2(comm, x); + } } double ArpackEigenvalueSolver::GetError(int i, EigenvalueSolver::ErrorType type) const @@ -454,6 +472,18 @@ double ArpackEigenvalueSolver::GetError(int i, EigenvalueSolver::ErrorType type) return 0.0; } +void ArpackEigenvalueSolver::RescaleEigenvectors(int num_eig) +{ + res = std::make_unique(num_eig); + xscale = std::make_unique(num_eig); + for (int i = 0; i < num_eig; i++) + { + x1.Set(V.get() + i * n, n); + xscale.get()[i] = 1.0 / GetEigenvectorNorm(x1, y1); + res.get()[i] = GetResidualNorm(eig.get()[i], x1, y1) / linalg::Norml2(comm, x1); + } +} + // EPS specific methods ArpackEPSSolver::ArpackEPSSolver(MPI_Comm comm, int print) @@ -483,9 +513,9 @@ void ArpackEPSSolver::SetOperators(const ComplexOperator &K, const ComplexOperat } // Set up workspace. - x.SetSize(opK->Height()); - y.SetSize(opK->Height()); - z.SetSize(opK->Height()); + x1.SetSize(opK->Height()); + y1.SetSize(opK->Height()); + z1.SetSize(opK->Height()); n = opK->Height(); } @@ -493,7 +523,7 @@ int ArpackEPSSolver::Solve() { // Set some defaults (default maximum iterations from SLEPc). CheckParameters(); - HYPRE_BigInt N = linalg::GlobalSize(comm, z); + HYPRE_BigInt N = linalg::GlobalSize(comm, z1); if (ncv > N) { ncv = mfem::internal::to_int(N); @@ -525,21 +555,13 @@ int ArpackEPSSolver::Solve() { eig = std::make_unique[]>(nev + 1); perm = std::make_unique(nev); - res = std::make_unique(nev); } // Solve the generalized eigenvalue problem. int num_conv = SolveInternal(n, r.get(), V.get(), eig.get(), perm.get()); // Compute the eigenpair residuals: || (K - λ M) x ||₂ for eigenvalue λ. - for (int i = 0; i < nev; i++) - { - const std::complex l = eig.get()[i]; - x.Set(V.get() + i * n, n); - opK->Mult(x, y); - opM->AddMult(x, y, -l); - res.get()[i] = linalg::Norml2(comm, y); - } + RescaleEigenvectors(nev); // Reset for next solve. info = 0; @@ -553,37 +575,46 @@ void ArpackEPSSolver::ApplyOp(const std::complex *px, // y = M⁻¹ K x . // Case 2: Shift-and-invert spectral transformation (opInv = (K - σ M)⁻¹) // y = (K - σ M)⁻¹ M x . - x.Set(px, n); + x1.Set(px, n); if (!sinvert) { - opK->Mult(x, z); - opInv->Mult(z, y); - y *= 1.0 / gamma; + opK->Mult(x1, z1); + opInv->Mult(z1, y1); + y1 *= 1.0 / gamma; } else { - opM->Mult(x, z); - opInv->Mult(z, y); - y *= gamma; + opM->Mult(x1, z1); + opInv->Mult(z1, y1); + y1 *= gamma; } if (opProj) { - // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(comm, y)); - opProj->Mult(y); - // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(comm, y)); + // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(comm, y1)); + opProj->Mult(y1); + // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(comm, y1)); } - y.Get(py, n); + y1.Get(py, n); } void ArpackEPSSolver::ApplyOpB(const std::complex *px, std::complex *py) const { MFEM_VERIFY(opB, "No B operator for weighted inner product in ARPACK solve!"); - x.Set(px, n); - opB->Mult(x.Real(), y.Real()); - opB->Mult(x.Imag(), y.Imag()); - y *= delta * gamma; - y.Get(py, n); + x1.Set(px, n); + opB->Mult(x1.Real(), y1.Real()); + opB->Mult(x1.Imag(), y1.Imag()); + y1 *= delta * gamma; + y1.Get(py, n); +} + +double ArpackEPSSolver::GetResidualNorm(std::complex l, const ComplexVector &x, + ComplexVector &r) const +{ + // Compute the i-th eigenpair residual: || (K - λ M) x ||₂ for eigenvalue λ. + opK->Mult(x, r); + opM->AddMult(x, r, -l); + return linalg::Norml2(comm, r); } double ArpackEPSSolver::GetBackwardScaling(std::complex l) const @@ -638,7 +669,7 @@ void ArpackPEPSolver::SetOperators(const ComplexOperator &K, const ComplexOperat x2.SetSize(opK->Height()); y1.SetSize(opK->Height()); y2.SetSize(opK->Height()); - z.SetSize(opK->Height()); + z1.SetSize(opK->Height()); n = opK->Height(); } @@ -647,7 +678,7 @@ int ArpackPEPSolver::Solve() // Set some defaults (from SLEPc ARPACK interface). The problem size is the size of the // 2x2 block linearized problem. CheckParameters(); - HYPRE_BigInt N = linalg::GlobalSize(comm, z); + HYPRE_BigInt N = linalg::GlobalSize(comm, z1); if (ncv > 2 * N) { ncv = mfem::internal::to_int(2 * N); @@ -682,25 +713,20 @@ int ArpackPEPSolver::Solve() if (!eig) { eig = std::make_unique[]>(nev + 1); - perm = std::make_unique(nev + 1); - res = std::make_unique(nev + 1); + perm = std::make_unique(nev); } // Solve the linearized eigenvalue problem. int num_conv = SolveInternal(2 * n, s.get(), W.get(), eig.get(), perm.get()); // Extract the eigenvector from the linearized problem and compute the eigenpair - // residuals: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for eigenvalue λ. + // residuals: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for eigenvalue λ. For the + // linearized problem, select the most accurate x for y = [x₁; x₂]. Or, just take x = x₁. for (int i = 0; i < nev; i++) { - const std::complex &l = eig.get()[i]; - ExtractEigenvector(l, W.get() + i * 2 * n, V.get() + i * n); - x1.Set(V.get() + i * n, n); - opK->Mult(x1, y1); - opC->AddMult(x1, y1, l); - opM->AddMult(x1, y1, l * l); - res.get()[i] = linalg::Norml2(comm, y1); + std::copy(W.get() + i * 2 * n, W.get() + (i * 2 + 1) * n, V.get() + i * n); } + RescaleEigenvectors(nev); // Reset for next solve. info = 0; @@ -729,9 +755,9 @@ void ArpackPEPSolver::ApplyOp(const std::complex *px, // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(comm, y1)); } - opK->Mult(x1, z); - opC->AddMult(x2, z, std::complex(gamma, 0.0)); - opInv->Mult(z, y2); + opK->Mult(x1, z1); + opC->AddMult(x2, z1, std::complex(gamma, 0.0)); + opInv->Mult(z1, y2); y2 *= -1.0 / (gamma * gamma); if (opProj) { @@ -743,9 +769,9 @@ void ArpackPEPSolver::ApplyOp(const std::complex *px, else { y2.AXPBYPCZ(sigma, x1, gamma, x2, 0.0); // Just temporarily - opM->Mult(y2, z); - opC->AddMult(x1, z, std::complex(1.0, 0.0)); - opInv->Mult(z, y1); + opM->Mult(y2, z1); + opC->AddMult(x1, z1, std::complex(1.0, 0.0)); + opInv->Mult(z1, y1); y1 *= -gamma; if (opProj) { @@ -782,6 +808,17 @@ void ArpackPEPSolver::ApplyOpB(const std::complex *px, y2.Get(py + n, n); } +double ArpackPEPSolver::GetResidualNorm(std::complex l, const ComplexVector &x, + ComplexVector &r) const +{ + // Compute the i-th eigenpair residual: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for + // eigenvalue λ. + opK->Mult(x, r); + opC->AddMult(x, r, l); + opM->AddMult(x, r, l * l); + return linalg::Norml2(comm, r); +} + double ArpackPEPSolver::GetBackwardScaling(std::complex l) const { // Make sure not to use norms from scaling as this can be confusing if they are different. @@ -802,24 +839,6 @@ double ArpackPEPSolver::GetBackwardScaling(std::complex l) const return normK + t * normC + t * t * normM; } -void ArpackPEPSolver::ExtractEigenvector(std::complex l, - const std::complex *py, - std::complex *px) const -{ - // Select the most accurate x for y = [x₁; x₂] from the linearized eigenvalue problem. Or, - // just take x = x₁. - x1.Set(py, n); - if (opB) - { - linalg::Normalize(comm, x1, *opB, y1); - } - else - { - linalg::Normalize(comm, x1); - } - x1.Get(px, n); -} - } // namespace palace::arpack #endif diff --git a/palace/linalg/arpack.hpp b/palace/linalg/arpack.hpp index 25771b735..13fc6c114 100644 --- a/palace/linalg/arpack.hpp +++ b/palace/linalg/arpack.hpp @@ -66,8 +66,8 @@ class ArpackEigenvalueSolver : public EigenvalueSolver // Storage for Arnoldi basis vectors. std::unique_ptr[]> V; - // Storage for computed residual norms. - std::unique_ptr res; + // Storage for computed residual norms and eigenvector scalings. + std::unique_ptr res, xscale; // On input used to define optional initial guess, on output stores final residual // vector. @@ -86,6 +86,9 @@ class ArpackEigenvalueSolver : public EigenvalueSolver // which case identity is used. const Operator *opB; + // Workspace vector for operator applications. + mutable ComplexVector x1, y1, z1; + // Perform the ARPACK RCI loop. int SolveInternal(int n, std::complex *r, std::complex *V, std::complex *eig, int *perm); @@ -97,6 +100,13 @@ class ArpackEigenvalueSolver : public EigenvalueSolver virtual void ApplyOp(const std::complex *px, std::complex *py) const = 0; virtual void ApplyOpB(const std::complex *px, std::complex *py) const = 0; + // Helper routine for computing the eigenvector normalization. + double GetEigenvectorNorm(const ComplexVector &x, ComplexVector &Bx) const; + + // Helper routine for computing the eigenpair residual. + virtual double GetResidualNorm(std::complex l, const ComplexVector &x, + ComplexVector &r) const = 0; + // Helper routine for computing the backward error. virtual double GetBackwardScaling(std::complex l) const = 0; @@ -155,11 +165,17 @@ class ArpackEigenvalueSolver : public EigenvalueSolver // Get the corresponding eigenvalue. std::complex GetEigenvalue(int i) const override; - // Get the corresponding eigenvector. + // Get the corresponding eigenvector. Eigenvectors are normalized such that ||x||₂ = 1, + // unless the B-matrix is set for weighted inner products. void GetEigenvector(int i, ComplexVector &x) const override; // Get the corresponding eigenpair error. double GetError(int i, ErrorType type) const override; + + // Re-normalize the given number of eigenvectors, for example if the matrix B for weighted + // inner products has changed. This does not perform re-orthogonalization with respect to + // the new matrix, only normalization. + void RescaleEigenvectors(int num_eig) override; }; // Generalized eigenvalue problem solver: K x = λ M x . @@ -172,13 +188,13 @@ class ArpackEPSSolver : public ArpackEigenvalueSolver // Operator norms for scaling. mutable double normK, normM; - // Workspace vector for operator applications. - mutable ComplexVector x, y, z; - protected: void ApplyOp(const std::complex *px, std::complex *py) const override; void ApplyOpB(const std::complex *px, std::complex *py) const override; + double GetResidualNorm(std::complex l, const ComplexVector &x, + ComplexVector &r) const override; + double GetBackwardScaling(std::complex l) const override; const char *GetName() const override { return "EPS"; } @@ -204,16 +220,15 @@ class ArpackPEPSolver : public ArpackEigenvalueSolver mutable double normK, normC, normM; // Workspace vectors for operator applications. - mutable ComplexVector x1, x2, y1, y2, z; - - // Do eigenvector extraction from the linearized problem to the actual eigenvectors. - void ExtractEigenvector(std::complex l, const std::complex *py, - std::complex *px) const; + mutable ComplexVector x2, y2; protected: void ApplyOp(const std::complex *px, std::complex *py) const override; void ApplyOpB(const std::complex *px, std::complex *py) const override; + double GetResidualNorm(std::complex l, const ComplexVector &x, + ComplexVector &r) const override; + double GetBackwardScaling(std::complex l) const override; const char *GetName() const override { return "PEP"; } diff --git a/palace/linalg/eps.hpp b/palace/linalg/eps.hpp index 1754b19d8..8aa9d07ae 100644 --- a/palace/linalg/eps.hpp +++ b/palace/linalg/eps.hpp @@ -100,11 +100,17 @@ class EigenvalueSolver // Get the corresponding eigenvalue. virtual std::complex GetEigenvalue(int i) const = 0; - // Get the corresponding eigenvector. + // Get the corresponding eigenvector. Eigenvectors are normalized such that ||x||₂ = 1, + // unless the B-matrix is set for weighted inner products. virtual void GetEigenvector(int i, ComplexVector &x) const = 0; // Get the corresponding eigenpair error. virtual double GetError(int i, ErrorType type) const = 0; + + // Re-normalize the given number of eigenvectors, for example if the matrix B for weighted + // inner products has changed. This does not perform re-orthogonalization with respect to + // the new matrix, only normalization. + virtual void RescaleEigenvectors(int num_eig) = 0; }; } // namespace palace diff --git a/palace/linalg/slepc.cpp b/palace/linalg/slepc.cpp index f2b63ba3d..1c796edd0 100644 --- a/palace/linalg/slepc.cpp +++ b/palace/linalg/slepc.cpp @@ -378,6 +378,19 @@ void SlepcEigenvalueSolver::Customize() } } +PetscReal SlepcEigenvalueSolver::GetEigenvectorNorm(const ComplexVector &x, + ComplexVector &Bx) const +{ + if (opB) + { + return linalg::Norml2(GetComm(), x, *opB, Bx); + } + else + { + return linalg::Norml2(GetComm(), x); + } +} + PetscReal SlepcEigenvalueSolver::GetError(int i, EigenvalueSolver::ErrorType type) const { switch (type) @@ -392,6 +405,20 @@ PetscReal SlepcEigenvalueSolver::GetError(int i, EigenvalueSolver::ErrorType typ return 0.0; } +void SlepcEigenvalueSolver::RescaleEigenvectors(int num_eig) +{ + res = std::make_unique(num_eig); + xscale = std::make_unique(num_eig); + for (int i = 0; i < num_eig; i++) + { + xscale.get()[i] = 0.0; + GetEigenvector(i, x1); + xscale.get()[i] = 1.0 / GetEigenvectorNorm(x1, y1); + res.get()[i] = + GetResidualNorm(GetEigenvalue(i), x1, y1) / linalg::Norml2(GetComm(), x1); + } +} + // EPS specific methods SlepcEPSSolverBase::SlepcEPSSolverBase(MPI_Comm comm, int print, const std::string &prefix) @@ -594,11 +621,7 @@ int SlepcEPSSolverBase::Solve() } // Compute and store the eigenpair residuals. - res = std::make_unique(num_conv); - for (int i = 0; i < num_conv; i++) - { - res.get()[i] = GetResidualNorm(i); - } + RescaleEigenvectors(num_conv); return (int)num_conv; } @@ -624,6 +647,10 @@ void SlepcEPSSolverBase::GetEigenvector(int i, ComplexVector &x) const PalacePetscCall(VecGetArrayRead(v0, &pv0)); x.Set(pv0, n); PalacePetscCall(VecRestoreArrayRead(v0, &pv0)); + if (xscale.get()[i] > 0.0) + { + x *= xscale.get()[i]; + } } BV SlepcEPSSolverBase::GetBV() const @@ -694,8 +721,8 @@ void SlepcEPSSolver::SetOperators(const ComplexOperator &K, const ComplexOperato { PalacePetscCall(MatCreateVecs(A0, nullptr, &v0)); } - x.SetSize(opK->Height()); - y.SetSize(opK->Height()); + x1.SetSize(opK->Height()); + y1.SetSize(opK->Height()); // Configure linear solver for generalized problem or spectral transformation. This also // allows use of the divergence-free projector as a linear solve side-effect. @@ -718,14 +745,13 @@ void SlepcEPSSolver::SetBMat(const Operator &B) PalacePetscCall(BVSetMatrix(bv, B0, PETSC_FALSE)); } -PetscReal SlepcEPSSolver::GetResidualNorm(int i) const +PetscReal SlepcEPSSolver::GetResidualNorm(PetscScalar l, const ComplexVector &x, + ComplexVector &r) const { // Compute the i-th eigenpair residual: || (K - λ M) x ||₂ for eigenvalue λ. - PetscScalar l = GetEigenvalue(i); - GetEigenvector(i, x); - opK->Mult(x, y); - opM->AddMult(x, y, -l); - return linalg::Norml2(GetComm(), y); + opK->Mult(x, r); + opM->AddMult(x, r, -l); + return linalg::Norml2(GetComm(), r); } PetscReal SlepcEPSSolver::GetBackwardScaling(PetscScalar l) const @@ -862,27 +888,21 @@ void SlepcPEPLinearSolver::GetEigenvector(int i, ComplexVector &x) const PalacePetscCall(VecGetArrayRead(v0, &pv0)); x.Set(pv0, n / 2); PalacePetscCall(VecRestoreArrayRead(v0, &pv0)); - - if (opB) - { - linalg::Normalize(GetComm(), x, *opB, y1); - } - else + if (xscale.get()[i] > 0.0) { - linalg::Normalize(GetComm(), x); + x *= xscale.get()[i]; } } -PetscReal SlepcPEPLinearSolver::GetResidualNorm(int i) const +PetscReal SlepcPEPLinearSolver::GetResidualNorm(PetscScalar l, const ComplexVector &x, + ComplexVector &r) const { // Compute the i-th eigenpair residual: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for // eigenvalue λ. - PetscScalar l = GetEigenvalue(i); - GetEigenvector(i, x1); - opK->Mult(x1, y1); - opC->AddMult(x1, y1, l); - opM->AddMult(x1, y1, l * l); - return linalg::Norml2(GetComm(), y1); + opK->Mult(x, r); + opC->AddMult(x, r, l); + opM->AddMult(x, r, l * l); + return linalg::Norml2(GetComm(), r); } PetscReal SlepcPEPLinearSolver::GetBackwardScaling(PetscScalar l) const @@ -1102,11 +1122,7 @@ int SlepcPEPSolverBase::Solve() } // Compute and store the eigenpair residuals. - res = std::make_unique(num_conv); - for (int i = 0; i < num_conv; i++) - { - res.get()[i] = GetResidualNorm(i); - } + RescaleEigenvectors(num_conv); return (int)num_conv; } @@ -1132,6 +1148,10 @@ void SlepcPEPSolverBase::GetEigenvector(int i, ComplexVector &x) const PalacePetscCall(VecGetArrayRead(v0, &pv0)); x.Set(pv0, n); PalacePetscCall(VecRestoreArrayRead(v0, &pv0)); + if (xscale.get()[i] > 0.0) + { + x *= xscale.get()[i]; + } } BV SlepcPEPSolverBase::GetBV() const @@ -1211,8 +1231,8 @@ void SlepcPEPSolver::SetOperators(const ComplexOperator &K, const ComplexOperato { PalacePetscCall(MatCreateVecs(A0, nullptr, &v0)); } - x.SetSize(opK->Height()); - y.SetSize(opK->Height()); + x1.SetSize(opK->Height()); + y1.SetSize(opK->Height()); // Configure linear solver. if (first) @@ -1234,16 +1254,15 @@ void SlepcPEPSolver::SetBMat(const Operator &B) PalacePetscCall(BVSetMatrix(bv, B0, PETSC_FALSE)); } -PetscReal SlepcPEPSolver::GetResidualNorm(int i) const +PetscReal SlepcPEPSolver::GetResidualNorm(PetscScalar l, const ComplexVector &x, + ComplexVector &r) const { // Compute the i-th eigenpair residual: || P(λ) x ||₂ = || (K + λ C + λ² M) x ||₂ for // eigenvalue λ. - PetscScalar l = GetEigenvalue(i); - GetEigenvector(i, x); - opK->Mult(x, y); - opC->AddMult(x, y, l); - opM->AddMult(x, y, l * l); - return linalg::Norml2(GetComm(), y); + opK->Mult(x, r); + opC->AddMult(x, r, l); + opM->AddMult(x, r, l * l); + return linalg::Norml2(GetComm(), r); } PetscReal SlepcPEPSolver::GetBackwardScaling(PetscScalar l) const @@ -1280,15 +1299,15 @@ PetscErrorCode __mat_apply_EPS_A0(Mat A, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opK->Mult(ctx->x, ctx->y); - ctx->y *= ctx->delta; + ctx->opK->Mult(ctx->x1, ctx->y1); + ctx->y1 *= ctx->delta; PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1306,15 +1325,15 @@ PetscErrorCode __mat_apply_EPS_A1(Mat A, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opM->Mult(ctx->x, ctx->y); - ctx->y *= ctx->delta * ctx->gamma; + ctx->opM->Mult(ctx->x1, ctx->y1); + ctx->y1 *= ctx->delta * ctx->gamma; PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1332,16 +1351,16 @@ PetscErrorCode __mat_apply_EPS_B(Mat A, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opB->Mult(ctx->x.Real(), ctx->y.Real()); - ctx->opB->Mult(ctx->x.Imag(), ctx->y.Imag()); - ctx->y *= ctx->delta * ctx->gamma; + ctx->opB->Mult(ctx->x1.Real(), ctx->y1.Real()); + ctx->opB->Mult(ctx->x1.Imag(), ctx->y1.Imag()); + ctx->y1 *= ctx->delta * ctx->gamma; PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1362,28 +1381,28 @@ PetscErrorCode __pc_apply_EPS(PC pc, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opInv->Mult(ctx->x, ctx->y); + ctx->opInv->Mult(ctx->x1, ctx->y1); if (!ctx->sinvert) { - ctx->y *= 1.0 / (ctx->delta * ctx->gamma); + ctx->y1 *= 1.0 / (ctx->delta * ctx->gamma); } else { - ctx->y *= 1.0 / ctx->delta; + ctx->y1 *= 1.0 / ctx->delta; } if (ctx->opProj) { - // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y)); - ctx->opProj->Mult(ctx->y); - // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y)); + // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1)); + ctx->opProj->Mult(ctx->y1); + // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1)); } PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1570,14 +1589,14 @@ PetscErrorCode __mat_apply_PEP_A0(Mat A, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opK->Mult(ctx->x, ctx->y); + ctx->opK->Mult(ctx->x1, ctx->y1); PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1596,14 +1615,14 @@ PetscErrorCode __mat_apply_PEP_A1(Mat A, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opC->Mult(ctx->x, ctx->y); + ctx->opC->Mult(ctx->x1, ctx->y1); PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1622,14 +1641,14 @@ PetscErrorCode __mat_apply_PEP_A2(Mat A, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opM->Mult(ctx->x, ctx->y); + ctx->opM->Mult(ctx->x1, ctx->y1); PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1647,16 +1666,16 @@ PetscErrorCode __mat_apply_PEP_B(Mat A, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opB->Mult(ctx->x.Real(), ctx->y.Real()); - ctx->opB->Mult(ctx->x.Imag(), ctx->y.Imag()); - ctx->y *= ctx->delta * ctx->gamma; + ctx->opB->Mult(ctx->x1.Real(), ctx->y1.Real()); + ctx->opB->Mult(ctx->x1.Imag(), ctx->y1.Imag()); + ctx->y1 *= ctx->delta * ctx->gamma; PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); @@ -1678,28 +1697,28 @@ PetscErrorCode __pc_apply_PEP(PC pc, Vec x, Vec y) const PetscScalar *px; PetscCall(VecGetArrayRead(x, &px)); - ctx->x.Set(px, n); + ctx->x1.Set(px, n); PetscCall(VecRestoreArrayRead(x, &px)); - ctx->opInv->Mult(ctx->x, ctx->y); + ctx->opInv->Mult(ctx->x1, ctx->y1); if (!ctx->sinvert) { - ctx->y *= 1.0 / (ctx->delta * ctx->gamma * ctx->gamma); + ctx->y1 *= 1.0 / (ctx->delta * ctx->gamma * ctx->gamma); } else { - ctx->y *= 1.0 / ctx->delta; + ctx->y1 *= 1.0 / ctx->delta; } if (ctx->opProj) { - // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y)); - ctx->opProj->Mult(ctx->y); - // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y)); + // Mpi::Print(" Before projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1)); + ctx->opProj->Mult(ctx->y1); + // Mpi::Print(" After projection: {:e}\n", linalg::Norml2(ctx->GetComm(), ctx->y1)); } PetscScalar *py; PetscCall(VecGetArrayWrite(y, &py)); - ctx->y.Get(py, n); + ctx->y1.Get(py, n); PetscCall(VecRestoreArrayWrite(y, &py)); PetscFunctionReturn(0); diff --git a/palace/linalg/slepc.hpp b/palace/linalg/slepc.hpp index a6f3a2ba6..912462a1d 100644 --- a/palace/linalg/slepc.hpp +++ b/palace/linalg/slepc.hpp @@ -75,6 +75,9 @@ class SlepcEigenvalueSolver : public EigenvalueSolver JD }; + // Workspace vector for operator applications. + mutable ComplexVector x1, y1; + protected: // Control print level for debugging. int print; @@ -86,8 +89,8 @@ class SlepcEigenvalueSolver : public EigenvalueSolver PetscScalar sigma; bool sinvert, region; - // Storage for computed residual norms. - std::unique_ptr res; + // Storage for computed residual norms and eigenvector normalizations. + std::unique_ptr res, xscale; // Reference to linear solver used for operator action for M⁻¹ (with no spectral // transformation) or (K - σ M)⁻¹ (generalized EVP with shift-and- invert) or P(σ)⁻¹ @@ -112,8 +115,12 @@ class SlepcEigenvalueSolver : public EigenvalueSolver // Customize object with command line options set. virtual void Customize(); + // Helper routine for computing the eigenvector normalization. + PetscReal GetEigenvectorNorm(const ComplexVector &x, ComplexVector &Bx) const; + // Helper routine for computing the eigenpair residual. - virtual PetscReal GetResidualNorm(int i) const = 0; + virtual PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x, + ComplexVector &r) const = 0; // Helper routine for computing the backward error. virtual PetscReal GetBackwardScaling(PetscScalar l) const = 0; @@ -161,6 +168,11 @@ class SlepcEigenvalueSolver : public EigenvalueSolver // Get the corresponding eigenpair error. PetscReal GetError(int i, ErrorType type) const override; + // Re-normalize the given number of eigenvectors, for example if the matrix B for weighted + // inner products has changed. This does not perform re-orthogonalization with respect to + // the new matrix, only normalization. + void RescaleEigenvectors(int num_eig) override; + // Get the basis vectors object. virtual BV GetBV() const = 0; @@ -248,15 +260,13 @@ class SlepcEPSSolver : public SlepcEPSSolverBase // References to matrices defining the generalized eigenvalue problem (not owned). const ComplexOperator *opK, *opM; - // Workspace vector for operator applications. - mutable ComplexVector x, y; - private: // Operator norms for scaling. mutable PetscReal normK, normM; protected: - PetscReal GetResidualNorm(int i) const override; + PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x, + ComplexVector &r) const override; PetscReal GetBackwardScaling(PetscScalar l) const override; @@ -287,14 +297,15 @@ class SlepcPEPLinearSolver : public SlepcEPSSolverBase const ComplexOperator *opK, *opC, *opM; // Workspace vectors for operator applications. - mutable ComplexVector x1, x2, y1, y2; + mutable ComplexVector x2, y2; private: // Operator norms for scaling. mutable PetscReal normK, normC, normM; protected: - PetscReal GetResidualNorm(int i) const override; + PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x, + ComplexVector &r) const override; PetscReal GetBackwardScaling(PetscScalar l) const override; @@ -383,15 +394,13 @@ class SlepcPEPSolver : public SlepcPEPSolverBase // (not owned). const ComplexOperator *opK, *opC, *opM; - // Workspace vector for operator applications. - mutable ComplexVector x, y; - private: // Operator norms for scaling. mutable PetscReal normK, normC, normM; protected: - PetscReal GetResidualNorm(int i) const override; + PetscReal GetResidualNorm(PetscScalar l, const ComplexVector &x, + ComplexVector &r) const override; PetscReal GetBackwardScaling(PetscScalar l) const override; diff --git a/test/examples/ref/cavity/impedance/domain-E.csv b/test/examples/ref/cavity/impedance/domain-E.csv index b71518786..4813e8c9d 100644 --- a/test/examples/ref/cavity/impedance/domain-E.csv +++ b/test/examples/ref/cavity/impedance/domain-E.csv @@ -1,16 +1,16 @@ m, E_elec (J), E_mag (J), E_cap (J), E_ind (J), E_elec[1] (J), p_elec[1], E_mag[1] (J), p_mag[1] - 1.000000000e+00, +9.925512834e-04, +9.925513803e-04, +0.000000000e+00, +0.000000000e+00, +9.925512834e-04, +1.000000000e+00, +9.925513803e-04, +1.000000000e+00 - 2.000000000e+00, +4.729600954e-04, +4.729601412e-04, +0.000000000e+00, +0.000000000e+00, +4.729600954e-04, +1.000000000e+00, +4.729601412e-04, +1.000000000e+00 - 3.000000000e+00, +6.129365465e-04, +6.129366057e-04, +0.000000000e+00, +0.000000000e+00, +6.129365465e-04, +1.000000000e+00, +6.129366057e-04, +1.000000000e+00 - 4.000000000e+00, +8.142171330e-04, +8.142172140e-04, +0.000000000e+00, +0.000000000e+00, +8.142171330e-04, +1.000000000e+00, +8.142172140e-04, +1.000000000e+00 - 5.000000000e+00, +5.798891145e-04, +5.798891697e-04, +0.000000000e+00, +0.000000000e+00, +5.798891145e-04, +1.000000000e+00, +5.798891697e-04, +1.000000000e+00 - 6.000000000e+00, +5.592910111e-04, +5.592910638e-04, +0.000000000e+00, +0.000000000e+00, +5.592910111e-04, +1.000000000e+00, +5.592910638e-04, +1.000000000e+00 - 7.000000000e+00, +4.722364212e-04, +4.722364651e-04, +0.000000000e+00, +0.000000000e+00, +4.722364212e-04, +1.000000000e+00, +4.722364651e-04, +1.000000000e+00 - 8.000000000e+00, +6.116104938e-04, +6.116105484e-04, +0.000000000e+00, +0.000000000e+00, +6.116104938e-04, +1.000000000e+00, +6.116105484e-04, +1.000000000e+00 - 9.000000000e+00, +9.570970594e-04, +9.570971471e-04, +0.000000000e+00, +0.000000000e+00, +9.570970594e-04, +1.000000000e+00, +9.570971471e-04, +1.000000000e+00 - 1.000000000e+01, +1.045672238e-03, +1.045672332e-03, +0.000000000e+00, +0.000000000e+00, +1.045672238e-03, +1.000000000e+00, +1.045672332e-03, +1.000000000e+00 - 1.100000000e+01, +6.789660950e-04, +6.789661588e-04, +0.000000000e+00, +0.000000000e+00, +6.789660950e-04, +1.000000000e+00, +6.789661588e-04, +1.000000000e+00 - 1.200000000e+01, +5.477391878e-04, +5.477392355e-04, +0.000000000e+00, +0.000000000e+00, +5.477391878e-04, +1.000000000e+00, +5.477392355e-04, +1.000000000e+00 - 1.300000000e+01, +9.011021676e-04, +9.011022535e-04, +0.000000000e+00, +0.000000000e+00, +9.011021676e-04, +1.000000000e+00, +9.011022535e-04, +1.000000000e+00 - 1.400000000e+01, +8.989335892e-04, +8.989336730e-04, +0.000000000e+00, +0.000000000e+00, +8.989335892e-04, +1.000000000e+00, +8.989336730e-04, +1.000000000e+00 - 1.500000000e+01, +5.797360828e-04, +5.797361362e-04, +0.000000000e+00, +0.000000000e+00, +5.797360828e-04, +1.000000000e+00, +5.797361362e-04, +1.000000000e+00 + 1.000000000e+00, +9.139656208e-02, +9.139657101e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657101e-02, +1.000000000e+00 + 2.000000000e+00, +9.139656208e-02, +9.139657093e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657093e-02, +1.000000000e+00 + 3.000000000e+00, +9.139656208e-02, +9.139657091e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657091e-02, +1.000000000e+00 + 4.000000000e+00, +9.139656208e-02, +9.139657118e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657118e-02, +1.000000000e+00 + 5.000000000e+00, +9.139656208e-02, +9.139657079e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657079e-02, +1.000000000e+00 + 6.000000000e+00, +9.139656208e-02, +9.139657070e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657070e-02, +1.000000000e+00 + 7.000000000e+00, +9.139656208e-02, +9.139657058e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657058e-02, +1.000000000e+00 + 8.000000000e+00, +9.139656208e-02, +9.139657024e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657024e-02, +1.000000000e+00 + 9.000000000e+00, +9.139656208e-02, +9.139657046e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657046e-02, +1.000000000e+00 + 1.000000000e+01, +9.139656208e-02, +9.139657036e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657036e-02, +1.000000000e+00 + 1.100000000e+01, +9.139656208e-02, +9.139657066e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657066e-02, +1.000000000e+00 + 1.200000000e+01, +9.139656208e-02, +9.139657004e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657004e-02, +1.000000000e+00 + 1.300000000e+01, +9.139656208e-02, +9.139657079e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657079e-02, +1.000000000e+00 + 1.400000000e+01, +9.139656208e-02, +9.139657060e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657060e-02, +1.000000000e+00 + 1.500000000e+01, +9.139656208e-02, +9.139657051e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139657051e-02, +1.000000000e+00 diff --git a/test/examples/ref/cavity/impedance/error-indicators.csv b/test/examples/ref/cavity/impedance/error-indicators.csv index a16c3fcad..86328603d 100644 --- a/test/examples/ref/cavity/impedance/error-indicators.csv +++ b/test/examples/ref/cavity/impedance/error-indicators.csv @@ -1,2 +1,2 @@ Norm, Minimum, Maximum, Mean - +1.786937592e-03, +1.175148950e-04, +3.263378882e-04, +2.040337637e-04 + +9.719333717e-04, +2.308437043e-05, +1.479269298e-04, +1.078490659e-04 diff --git a/test/examples/ref/cavity/pec/domain-E.csv b/test/examples/ref/cavity/pec/domain-E.csv index a1c556d08..810dc0cd9 100644 --- a/test/examples/ref/cavity/pec/domain-E.csv +++ b/test/examples/ref/cavity/pec/domain-E.csv @@ -1,16 +1,16 @@ m, E_elec (J), E_mag (J), E_cap (J), E_ind (J), E_elec[1] (J), p_elec[1], E_mag[1] (J), p_mag[1] - 1.000000000e+00, +9.925512863e-04, +9.925513657e-04, +0.000000000e+00, +0.000000000e+00, +9.925512863e-04, +1.000000000e+00, +9.925513657e-04, +1.000000000e+00 - 2.000000000e+00, +4.728817282e-04, +4.728817661e-04, +0.000000000e+00, +0.000000000e+00, +4.728817282e-04, +1.000000000e+00, +4.728817661e-04, +1.000000000e+00 - 3.000000000e+00, +6.131076891e-04, +6.131077382e-04, +0.000000000e+00, +0.000000000e+00, +6.131076891e-04, +1.000000000e+00, +6.131077382e-04, +1.000000000e+00 - 4.000000000e+00, +8.142171353e-04, +8.142172004e-04, +0.000000000e+00, +0.000000000e+00, +8.142171353e-04, +1.000000000e+00, +8.142172004e-04, +1.000000000e+00 - 5.000000000e+00, +5.799958853e-04, +5.799959317e-04, +0.000000000e+00, +0.000000000e+00, +5.799958853e-04, +1.000000000e+00, +5.799959317e-04, +1.000000000e+00 - 6.000000000e+00, +5.593259764e-04, +5.593260210e-04, +0.000000000e+00, +0.000000000e+00, +5.593259764e-04, +1.000000000e+00, +5.593260210e-04, +1.000000000e+00 - 7.000000000e+00, +4.726815423e-04, +4.726815802e-04, +0.000000000e+00, +0.000000000e+00, +4.726815423e-04, +1.000000000e+00, +4.726815802e-04, +1.000000000e+00 - 8.000000000e+00, +6.133889633e-04, +6.133890123e-04, +0.000000000e+00, +0.000000000e+00, +6.133889633e-04, +1.000000000e+00, +6.133890123e-04, +1.000000000e+00 - 9.000000000e+00, +9.570973399e-04, +9.570974169e-04, +0.000000000e+00, +0.000000000e+00, +9.570973399e-04, +1.000000000e+00, +9.570974169e-04, +1.000000000e+00 - 1.000000000e+01, +1.045672121e-03, +1.045672205e-03, +0.000000000e+00, +0.000000000e+00, +1.045672121e-03, +1.000000000e+00, +1.045672205e-03, +1.000000000e+00 - 1.100000000e+01, +6.789660964e-04, +6.789661506e-04, +0.000000000e+00, +0.000000000e+00, +6.789660964e-04, +1.000000000e+00, +6.789661506e-04, +1.000000000e+00 - 1.200000000e+01, +5.477372102e-04, +5.477372540e-04, +0.000000000e+00, +0.000000000e+00, +5.477372102e-04, +1.000000000e+00, +5.477372540e-04, +1.000000000e+00 - 1.300000000e+01, +9.011033625e-04, +9.011034348e-04, +0.000000000e+00, +0.000000000e+00, +9.011033625e-04, +1.000000000e+00, +9.011034348e-04, +1.000000000e+00 - 1.400000000e+01, +8.989337903e-04, +8.989338623e-04, +0.000000000e+00, +0.000000000e+00, +8.989337903e-04, +1.000000000e+00, +8.989338623e-04, +1.000000000e+00 - 1.500000000e+01, +5.800575034e-04, +5.800575498e-04, +0.000000000e+00, +0.000000000e+00, +5.800575034e-04, +1.000000000e+00, +5.800575498e-04, +1.000000000e+00 + 1.000000000e+00, +9.139656208e-02, +9.139656939e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656939e-02, +1.000000000e+00 + 2.000000000e+00, +9.139656208e-02, +9.139656940e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656940e-02, +1.000000000e+00 + 3.000000000e+00, +9.139656208e-02, +9.139656940e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656940e-02, +1.000000000e+00 + 4.000000000e+00, +9.139656208e-02, +9.139656940e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656940e-02, +1.000000000e+00 + 5.000000000e+00, +9.139656208e-02, +9.139656940e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656940e-02, +1.000000000e+00 + 6.000000000e+00, +9.139656208e-02, +9.139656937e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656937e-02, +1.000000000e+00 + 7.000000000e+00, +9.139656208e-02, +9.139656940e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656940e-02, +1.000000000e+00 + 8.000000000e+00, +9.139656208e-02, +9.139656939e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656939e-02, +1.000000000e+00 + 9.000000000e+00, +9.139656208e-02, +9.139656944e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656944e-02, +1.000000000e+00 + 1.000000000e+01, +9.139656208e-02, +9.139656941e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656941e-02, +1.000000000e+00 + 1.100000000e+01, +9.139656208e-02, +9.139656939e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656939e-02, +1.000000000e+00 + 1.200000000e+01, +9.139656208e-02, +9.139656939e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656939e-02, +1.000000000e+00 + 1.300000000e+01, +9.139656208e-02, +9.139656942e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656942e-02, +1.000000000e+00 + 1.400000000e+01, +9.139656208e-02, +9.139656940e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656940e-02, +1.000000000e+00 + 1.500000000e+01, +9.139656208e-02, +9.139656940e-02, +0.000000000e+00, +0.000000000e+00, +9.139656208e-02, +1.000000000e+00, +9.139656940e-02, +1.000000000e+00 diff --git a/test/examples/ref/cavity/pec/error-indicators.csv b/test/examples/ref/cavity/pec/error-indicators.csv index 40ae37224..d673e73ad 100644 --- a/test/examples/ref/cavity/pec/error-indicators.csv +++ b/test/examples/ref/cavity/pec/error-indicators.csv @@ -1,2 +1,2 @@ Norm, Minimum, Maximum, Mean - +1.786988675e-03, +1.174739444e-04, +3.265272544e-04, +2.040325936e-04 + +9.719345476e-04, +2.308434728e-05, +1.479271253e-04, +1.078491905e-04