Direct solver for the block matrix system

[Hemanta822]

Hi All,

Is there any way to implement a direct solver (UMFPackSolver) in the block operator matrix system (Darcy code - example 5). In example 5, UMFPackSolver has been only used for a diagonal block of a preconditioner. Please let me know your thoughts.

[jandrej]

You can use HypreParMatrixFromBlocks to generate the monolithic HypreParMatrix.

[jandrej]

No

[pazner]

If you use BlockMatrix instead of BlockOperator (at least for the non-PA version), then BlockMatrix::CreateMonolithic will give you what you want.

[Hemanta822]

Thank you @pazner. I modified BlockOperator to BlockMatrix. It runs the code, but I see big drop in the H1 convergence rate for the displacement (u) at finer refinements. Do I need to modify anything in right side block?

[pazner]

No, you shouldn't have to make any further modifications. These changes to ex5 work for me:

diff --git a/examples/ex5.cpp b/examples/ex5.cpp
index 1a8b25094..1a0522fb4 100644
--- a/examples/ex5.cpp
+++ b/examples/ex5.cpp
@@ -180,126 +180,34 @@ int main(int argc, char *argv[])
    bVarf->Assemble();
    if (!pa) { bVarf->Finalize(); }
 
-   BlockOperator darcyOp(block_offsets);
-
-   TransposeOperator *Bt = NULL;
+   BlockMatrix darcyOp(block_offsets);
+   SparseMatrix *Bt = nullptr;
 
    if (pa)
    {
-      Bt = new TransposeOperator(bVarf);
-
-      darcyOp.SetBlock(0,0, mVarf);
-      darcyOp.SetBlock(0,1, Bt, -1.0);
-      darcyOp.SetBlock(1,0, bVarf, -1.0);
+      MFEM_ABORT("");
    }
    else
    {
       SparseMatrix &M(mVarf->SpMat());
       SparseMatrix &B(bVarf->SpMat());
       B *= -1.;
-      Bt = new TransposeOperator(&B);
+      Bt = Transpose(B);
 
       darcyOp.SetBlock(0,0, &M);
       darcyOp.SetBlock(0,1, Bt);
       darcyOp.SetBlock(1,0, &B);
    }
 
-   // 10. Construct the operators for preconditioner
-   //
-   //                 P = [ diag(M)         0         ]
-   //                     [  0       B diag(M)^-1 B^T ]
-   //
-   //     Here we use Symmetric Gauss-Seidel to approximate the inverse of the
-   //     pressure Schur Complement
-   SparseMatrix *MinvBt = NULL;
-   Vector Md(mVarf->Height());
-
-   BlockDiagonalPreconditioner darcyPrec(block_offsets);
-   Solver *invM, *invS;
-   SparseMatrix *S = NULL;
-
-   if (pa)
-   {
-      mVarf->AssembleDiagonal(Md);
-      auto Md_host = Md.HostRead();
-      Vector invMd(mVarf->Height());
-      for (int i=0; i<mVarf->Height(); ++i)
-      {
-         invMd(i) = 1.0 / Md_host[i];
-      }
-
-      Vector BMBt_diag(bVarf->Height());
-      bVarf->AssembleDiagonal_ADAt(invMd, BMBt_diag);
-
-      Array<int> ess_tdof_list;  // empty
-
-      invM = new OperatorJacobiSmoother(Md, ess_tdof_list);
-      invS = new OperatorJacobiSmoother(BMBt_diag, ess_tdof_list);
-   }
-   else
-   {
-      SparseMatrix &M(mVarf->SpMat());
-      M.GetDiag(Md);
-      Md.HostReadWrite();
-
-      SparseMatrix &B(bVarf->SpMat());
-      MinvBt = Transpose(B);
-
-      for (int i = 0; i < Md.Size(); i++)
-      {
-         MinvBt->ScaleRow(i, 1./Md(i));
-      }
-
-      S = Mult(B, *MinvBt);
-
-      invM = new DSmoother(M);
-
-#ifndef MFEM_USE_SUITESPARSE
-      invS = new GSSmoother(*S);
-#else
-      invS = new UMFPackSolver(*S);
-#endif
-   }
-
-   invM->iterative_mode = false;
-   invS->iterative_mode = false;
-
-   darcyPrec.SetDiagonalBlock(0, invM);
-   darcyPrec.SetDiagonalBlock(1, invS);
+   std::unique_ptr<SparseMatrix> monolithic(darcyOp.CreateMonolithic());
+   UMFPackSolver solver(*monolithic);
 
-   // 11. Solve the linear system with MINRES.
-   //     Check the norm of the unpreconditioned residual.
-   int maxIter(1000);
-   real_t rtol(1.e-6);
-   real_t atol(1.e-10);
-
-   chrono.Clear();
    chrono.Start();
-   MINRESSolver solver;
-   solver.SetAbsTol(atol);
-   solver.SetRelTol(rtol);
-   solver.SetMaxIter(maxIter);
-   solver.SetOperator(darcyOp);
-   solver.SetPreconditioner(darcyPrec);
-   solver.SetPrintLevel(1);
-   x = 0.0;
    solver.Mult(rhs, x);
    if (device.IsEnabled()) { x.HostRead(); }
    chrono.Stop();
 
-   if (solver.GetConverged())
-   {
-      std::cout << "MINRES converged in " << solver.GetNumIterations()
-                << " iterations with a residual norm of "
-                << solver.GetFinalNorm() << ".\n";
-   }
-   else
-   {
-      std::cout << "MINRES did not converge in " << solver.GetNumIterations()
-                << " iterations. Residual norm is " << solver.GetFinalNorm()
-                << ".\n";
-   }
-   std::cout << "MINRES solver took " << chrono.RealTime() << "s.\n";
+   std::cout << "Direct solver took " << chrono.RealTime() << "s.\n";
 
    // 12. Create the grid functions u and p. Compute the L2 error norms.
    GridFunction u, p;
@@ -372,11 +280,7 @@ int main(int argc, char *argv[])
    // 17. Free the used memory.
    delete fform;
    delete gform;
-   delete invM;
-   delete invS;
-   delete S;
    delete Bt;
-   delete MinvBt;
    delete mVarf;
    delete bVarf;
    delete W_space;

[Hemanta822]

Thank you @pazner. This helps me a lot. :-)