[Dune] Parallel solvers in DUNE

[Markus Blatt]

Hi Max,

welcome to Dune.
On Fri, Jun 18, 2021 at 11:01:23AM +0000, Firmbach, Max wrote:
>I have some questions regarding parallel solvers in DUNE. I want to solve a linear systems of equations in parallel
>
>and faced some difficulties trying to get a nonoverlapping version to work ...
>
>My approach is to split my mesh up with the loadbalance command and afterwards assemble the stiffness matrix
>
>on each local partition using the "interiorBorder (hope that is correct)" partition type. Afterwards the boundary
>
>conditions are set.
>

I assume your are not using DUNE for the assembly. What kind of grid/diiscretization are you using? Where are the degrees of
freedom attached?

Since you are using the InteriorBorder, I assume that your degrees of freedom are attached to the vertices.

Note that our naming is a bit misleading. It rather a way for dealing with DUNE grids that do not support overlap but
use ghost cells. E.g. for allow Alugrid one could not identify vertices of different ghost cells that are actually the
same vertex. Those have a different index.

Just look into the code to get an idea what we have to do to accomdate this.

If you have a grid without any overlap/ghost (and do not need powerful preconditioners) you could implement your own
ScalarProduct/Preconditioner/Operator in an additive Schwarz sense. Oliver has an example in his fabulous book chapter 7.4.4
https://www.springer.com/de/book/9783030597016

If your grid is able to have one layer of overlap, then you can simply use your normal assembly and the Overlapping versions.

Otherwise I would recommend, to algraically extend the linear system such that:
- Each row is assigned to one process (marked as owner in the index set), all other are marked as copy
- Each row owned is consistent, i.e. Has the same sparsity and entries as in a serial run. I.e. you will need
  add entries for unknowns at the border here.
  (For a real partitioning you will need to add up the values from your assembly on InteriorBorder here)
- Set up the Index set accordingly.

There should be code in e.g. Dumux that does that. Seems like we never found time to add this to dune-istl. Such a pity.

>
>Considering the preconditioner and solver, I define the following things:
>
>
>  using Communication = OwnerOverlapCopyCommunication<int, int>;
>  using SSOR                    = SeqSSOR<operatorType, blockVector, blockVector>;
>  using Preconditioner   = NonoverlappingBlockPreconditioner<Communication, SSOR>;
>  using ScalarProduct     = NonoverlappingSchwarzScalarProduct<blockVector, Communication>;
>  using LinearOperator  = NonoverlappingSchwarzOperator<operatorType, blockVector, blockVector, Communication>;
>
>  Communication  comm; comm.remoteIndices().rebuild<true>();
>  SSOR                    ssor(stiffnessMatrix, 3, 0.7);
>  Preconditioner  preconditioner(ssor, comm);
>  ScalarProduct    scalarProduct(comm);
>  LinearOperator linearOperator(stiffnessMatrix, comm);
>
>  CGSolver<blockVector> solver(linearOperator, scalarProduct, preconditioner, 1e-8, 1000, (mpiRank==0) ? 2 : 0);
>
>  InverseOperatorResult statistics;
>  solver.apply(displacementVector, loadVector, statistics);
>
>
>For a single process I get the right solution, but for more processes I get a mapping segmentation fault. 

I am not sure what that means. More details, please?


Markus

-- 

Dr. Markus Blatt - HPC-Simulation-Software & Services http://www.dr-blatt.de
Pedettistr. 38, 85072 Eichstätt, Germany,  USt-Id: DE279960836
Tel.: +49 (0) 160 97590858