Direct solver for the block matrix system #4699
Replies: 1 comment 6 replies
-
|
You can use |
Beta Was this translation helpful? Give feedback.
-
|
No |
Beta Was this translation helpful? Give feedback.
-
|
If you use |
Beta Was this translation helpful? Give feedback.
-
|
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? |
Beta Was this translation helpful? Give feedback.
-
|
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; |
Beta Was this translation helpful? Give feedback.
-
|
Thank you @pazner. This helps me a lot. :-) |
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
Uh oh!
There was an error while loading. Please reload this page.
-
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.
Beta Was this translation helpful? Give feedback.
All reactions