Mfem: Direct MUMPS solver via PETSc for Complex OperatorHandle(ex22p) or BlockOperator(ex5p)

[hongbo-yao]

Hi:

I'm solving the complex linear system (ex22p) by Direct MUMPS solver via PETSc. My way is to change the complex linear system into 2x2 BlockOperator, and create a PetscParMatrix based on this BlockOperator. Then create PetscLinearSolver and call MUMPS to solve this linear system. The code:

`ComplexHypreParMatrix * complex_mat = dynamic_cast(A.Ptr());
HypreParMatrix Mat_r = complex_mat->real();
HypreParMatrix Mat_i = complex_mat->imag();

Array block_trueOffsets(3); // number of variables + 1
block_trueOffsets[0] = 0;
block_trueOffsets[1] = pfes->TrueVSize();
block_trueOffsets[2] = pfes->TrueVSize();
block_trueOffsets.PartialSum();
BlockOperator* Op = new BlockOperator(block_trueOffsets);
Op->SetDiagonalBlock(0, &Mat_r, 1.0);
Op->SetDiagonalBlock(1, &Mat_r, -1.0);
Op->SetBlock (0, 1, &Mat_i, -1.0);
Op->SetBlock (1, 0, &Mat_i, -1.0);

const char *petscrc_file = "petsc_opts";
MFEMInitializePetsc(NULL, NULL, petscrc_file, NULL);

PetscParMatrix* _A = new PetscParMatrix(pmesh->GetComm(), Op, Operator::PETSC_MATAIJ);
PetscLinearSolver* solver = new PetscLinearSolver(pmesh->GetComm());
PetscPreconditioner* prec = new PetscPreconditioner(_A,"solver_");
solver->SetOperator(_A);
solver->SetPreconditioner(*prec);
solver->SetAbsTol(1.e-10);
solver->SetRelTol(1.e-6);
solver->SetMaxIter(500);
solver->SetPrintLevel(2);
solver->Mult(B, X);`

Unfortunately, it prints error:
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for possible LU and Cholesky solvers
[0]PETSC ERROR: MatSolverType mumps does not support matrix type nest
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.12.0, Sep, 29, 2019
[0]PETSC ERROR: Unknown Name on a arch-darwin-c-opt named MacBook-Pro by yaohb Tue Apr 21 20:00:03 2020
[0]PETSC ERROR: Configure options --download-mumps --download-scalapack --with-debugging=no
[0]PETSC ERROR: #1 MatGetFactor() line 4419 in /petsc-3.12.0/src/mat/interface/matrix.c
[0]PETSC ERROR: #2 PCSetUp_LU() line 88 in /petsc-3.12.0/src/ksp/pc/impls/factor/lu/lu.c
[0]PETSC ERROR: #3 PCSetUp() line 894 in /petsc-3.12.0/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #4 KSPSetUp() line 377 in /petsc-3.12.0/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #5 KSPSolve() line 707 in /petsc-3.12.0/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #6 void mfem::PetscLinearSolver::MultKernel(const mfem::Vector &, mfem::Vector &, bool) const() line 2574 in linalg/petsc.cpp

MFEM abort: Error in PETSc. See stacktrace above.
... in function: void mfem::PetscLinearSolver::MultKernel(const mfem::Vector &, mfem::Vector &, bool) const
... in file: linalg/petsc.cpp:2574

application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0

This question is similar to #282 but it failed.

---

So:

1. I have an idea, that is to create a HypreParMatrix from this 2x2 blocks using #1422, it may work?

1. Can you give me some advice on solving ex22p through parallel direct solver?

Thanks!
Hongbo

[psocratis]

Hi @hongbo-yao, you can follow ex25p that solves Maxwell equations using superLU.
If you want to use MUMPS through petsc you can do something like

cpp HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix(); const char *petscrc_file = "petscrc"; MFEMInitializePetsc(NULL, NULL, petscrc_file, NULL); PetscLinearSolver * petsc = new PetscLinearSolver(MPI_COMM_WORLD); // Convert to PetscParMatrix petsc->SetOperator(PetscParMatrix(A, Operator::PETSC_MATAIJ)); petsc->Mult(B, X); MFEMFinalizePetsc();
The file petscrc just contains the petsc options

-pc_type lu
-pc_factor_mat_solver_type mumps

Hope this helps,
Socratis

[hongbo-yao]

Hi @psocratis, Thanks a lot! MUMPS works... I will try SuperLU follow ex25p ASAP.

Hongbo

[psocratis]

Hi @hongbo-yao , I am glad it works. MUMPS is expected to be slower since it's a direct solver with much higher complexity. I expect that SuperLU will give similar timings but it's worth trying it.

For the AMS I am assuming you create the AMS preconditioner for the diagonal blocks that correspond to the real part of the matrix. Is that the case? What exactly is the problem you are solving? There are ways to extract all the components of the AMG solver and modify accordingly (see in hypre.cpp). I experimented a bit in the past with that for block operators coming from Least Squares formulation. You might find the following helpful; https://github.com/mfem/mfem/blob/gmg-solver/examples/waves/LS_blkAMS.cpp

Socratis

[hongbo-yao]

Thanks for your quick reply, @psocratis !

I may need more time to think about it on my own. And I will learn more about AMG solver and consider to modify it follow your advice.

For other problems, maybe I will ask you for help. I plan to close this issue for now.

Again, Thanks for your help!
Best wishes in this tense period!
Hongbo