diff --git a/palace/linalg/mumps.hpp b/palace/linalg/mumps.hpp index f84308939..57811b2c3 100644 --- a/palace/linalg/mumps.hpp +++ b/palace/linalg/mumps.hpp @@ -29,9 +29,9 @@ class MumpsSolver : public mfem::MUMPSSolver iodata.problem.type == ProblemType::ELECTROSTATIC || iodata.problem.type == ProblemType::MAGNETOSTATIC) ? mfem::MUMPSSolver::SYMMETRIC_POSITIVE_DEFINITE - : iodata.solver.linear.complex_coarse_solve - ? mfem::MUMPSSolver::UNSYMMETRIC - : mfem::MUMPSSolver::SYMMETRIC_INDEFINITE, + : iodata.boundaries.periodic.wave_vector == std::array{0.0, 0.0, 0.0} + ? mfem::MUMPSSolver::SYMMETRIC_INDEFINITE + : mfem::MUMPSSolver::UNSYMMETRIC, iodata.solver.linear.sym_factorization, (iodata.solver.linear.strumpack_compression_type == SparseCompression::BLR) ? iodata.solver.linear.strumpack_lr_tol diff --git a/palace/linalg/solver.cpp b/palace/linalg/solver.cpp index 357b524a4..8c2539f75 100644 --- a/palace/linalg/solver.cpp +++ b/palace/linalg/solver.cpp @@ -55,8 +55,9 @@ void MfemWrapperSolver::SetOperator(const ComplexOperator &op) { if (complex_matrix) { - // A = [Ar, -Ai] - // [Ai, Ar] + // A = [Ar, Ai] + // [Ai, -Ar] + // We solve A [xr; -xi] = [br; bi] mfem::Array2D blocks(2, 2); mfem::Array2D block_coeffs(2, 2); blocks(0, 0) = hAr; @@ -64,9 +65,9 @@ void MfemWrapperSolver::SetOperator(const ComplexOperator &op) blocks(1, 0) = hAi; blocks(1, 1) = hAr; block_coeffs(0, 0) = 1.0; - block_coeffs(0, 1) = -1.0; + block_coeffs(0, 1) = 1.0; block_coeffs(1, 0) = 1.0; - block_coeffs(1, 1) = 1.0; + block_coeffs(1, 1) = -1.0; A.reset(mfem::HypreParMatrixFromBlocks(blocks, &block_coeffs)); } else @@ -168,6 +169,8 @@ void MfemWrapperSolver::Mult(const ComplexVector &x, Y.ReadWrite(); yr.MakeRef(Y, 0, Ny); yi.MakeRef(Y, Ny, Ny); + // [yr; yi] is the complex conjugate of the solution + yi *= -1.0; y.Real() = yr; y.Imag() = yi; }