[Dune] Runtime error with parallel linear solver

[Tobias Ritschel]

I was a bit too quick to ask. I found out how to synchronize the vectors
after calling apply on the solver object.

I stumbled upon this blog
<http://marcoagnese.blogspot.dk/2014/06/socis-2014-dune-parallel-communication.html>
which has the example parallel.cc that seems to be the same that is
discussed in the tutorial on "parallel communication interface
<http://www.dune-project.org/doc/comm/communication.pdf>". It also links to
the research article "C++ Components Describing Parallel Domain
Decompositionand Communication
<http://conan.iwr.uni-heidelberg.de/people/peter/pdf/BlattBastian_IJPEDS_2009.pdf>"
which seems to be a newer version of the tutorial.

Just thought I'd let you know.

On 13 June 2015 at 23:18, Tobias Ritschel <tobiasritschel at gmail.com> wrote:

> Hi again
>
> I managed to create an index set and rebuild the RemoteIndices such that
> the code compiles and run without errors.
>
> I am not sure exactly how to use the library though. The system has size
> 2*Sim.nc x 2*Sim.nc and is created with opm autodiff and hence has the
> SparseMatrix
> <http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html> format
> from Eigen. On each processor I create a Communication object and create an
> index set similar to Is from Figure 1 in communicator.pdf
> <https://www.dune-project.org/doc/comm/communication.pdf> with the
> following code.
>
> int rank;
> int procs;
> MPI_Comm_rank(MPI_COMM_WORLD, &rank);
> MPI_Comm_size(MPI_COMM_WORLD, &procs);
>
> /* Create communication structure */
> Comm comm(MPI_COMM_WORLD);
>  typedef Dune::OwnerOverlapCopyAttributeSet::AttributeSet AS;
> typedef Dune::ParallelLocalIndex<AS> LocalIndex;
> typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet
> PIndexSet;
> typedef Dune::RemoteIndices<PIndexSet> RIndices;
>
> PIndexSet &IS = comm.indexSet();
> int elem_per_proc = 2*Sim.nc/procs; /* Truncated */
> int offset = 0;
>
> IS.beginResize();
>
> int start = elem_per_proc*rank;
>
> /* Last processor takes the rest */
> int end = elem_per_proc*(rank+1);
> if(rank == procs-1)
> end = 2*Sim.nc;
>
> /* Add the previous element as a ghost point if not on first processor */
> if(rank > 0){
> IS.add(start-1, LocalIndex(0, AS::overlap));
> offset = 1;
> }
>
> /* Indices owned by this process */
> for(int i = start; i < end; ++i){
> IS.add(i, LocalIndex(i-start+offset, AS::owner));
> }
>
> /* Add the following element as a ghost point if not on last processor */
> if(rank < procs-1)
> IS.add(end, LocalIndex(end-start+offset, AS::overlap));
>
> IS.endResize();
>
> RIndices &RI = comm.remoteIndices();
> RI.rebuild<true>();
>
> When I need to solve the linear system I first create an OPM DuneMatrix
> <http://www.opm-project.org/documentation/master/opm-autodiff/html/class_opm_1_1_dune_matrix.php> which
> I then use to construct an object of the type on each processor
>
> Dune::OverlappingSchwarzOperator<Opm::DuneMatrix, Vector, Vector, Comm>
>
> Then I construct the sequential and parallel preconditioners and the
> scalarproduct with the types
>
> typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
> typedef Dune::BlockVector<Dune::FieldVector<double,1>> Vector;
> typedef Dune::OverlappingSchwarzOperator<Opm::DuneMatrix, Vector, Vector,
> Comm> Matrix;
>
> typedef Dune::OverlappingSchwarzScalarProduct<Vector, Comm> ScalarProduct;
> typedef Dune::SeqILU0<Opm::DuneMatrix, Vector, Vector> SeqPrec;
> typedef Dune::BlockPreconditioner<Vector, Vector, Comm, SeqPrec> ParPrec;
>
> Lastly I call apply. I also tried to run all code but the part creating
> the communicator and the index set on the processor with rank 0 but then
> the program hangs when apply is called on the solver object.
>
> I guess I am missing some synchronisation or communication between the
> processors? For instance, when apply is called, I guess it only updates the
> part that is owned by that processor. By the way, I realize that speedup
> attained from this approach may be limited since the system matrix is
> created sequentially. It is more to show the possibility of using Dune ISTL
> in parallel.
>
> On 11 June 2015 at 11:54, Tobias Ritschel <tobiasritschel at gmail.com>
> wrote:
>
>> Thanks, I will have a look at the documentation. I solve for both
>> pressure and saturation in each cell, so I I'll create the index set.
>>
>> On 11 June 2015 at 11:36, Markus Blatt <markus at dr-blatt.de> wrote:
>>
>>> On Thu, Jun 11, 2015 at 10:38:10AM +0200, Tobias Ritschel wrote:
>>> > Ah, this is clearly my error. I am using the constructor
>>> >
>>> > OwnerOverlapCopyCommunication
>>> > <
>>> http://www.dune-project.org/doc/doxygen/dune-istl-html/class_dune_1_1_owner_overlap_copy_communication.html#a93f446939bd4d6623bdf2e544f183f8c
>>> >
>>> > (MPI_Comm
>>> > comm_, SolverCategory::Category
>>> > <
>>> http://www.dune-project.org/doc/doxygen/dune-istl-html/struct_dune_1_1_solver_category.html#ae061380ac961490be6ee353cf0dc1733
>>> >
>>> > cat_=SolverCategory::overlapping,
>>> > bool freecomm_=false)
>>> >
>>> > and used the created object without an index set. Is the information
>>> > automatically set up when constructing an IndexInfoFromGrid object, or
>>> > should they be set with AddLocalIndex and AddRemoteIndex? And is this
>>> the
>>> > most convenient way to set up the indices?
>>>
>>> Don't use IndexInfoGrid. I assume that is a left over from ancient
>>> dune-disc times. You can access the underlying index set by
>>> OwnerOverlapCopyCommunication::indexSet() and use that to add the
>>> indices with the correct tag. Later on get the remote indices with
>>> OwnerOverlapCopyCommunication::remoteIndices() and call rebuild on it,
>>>
>>> You will find more information on the index sets and how they are used
>>> in dune-istl in the "Description of the parallel communication
>>> interface" on https://dune-project.org/doc/index.html. You need one
>>> index per degree of freedom and use the attribute owner to mark a
>>> partitioning. In the easiest case you need to mark connected degrees
>>> of freedom that are not owner on the process as copy.
>>>
>>> If you are using dune-cornerpoint with one unknown per cell, then you
>>> can use ParallelISTLInformation::copyValuesTo from opm-core to do the
>>> job.
>>>
>>> Otherwise there is some manual work to do. E.g. compute the owner
>>> region using grid communication.
>>>
>>> Markus
>>>
>>> >
>>> > On 11 June 2015 at 09:59, Markus Blatt <markus at dr-blatt.de> wrote:
>>> >
>>> > > Hi,
>>> > >
>>> > > On Wed, Jun 10, 2015 at 09:00:36PM +0200, Tobias Ritschel wrote:
>>> > > >
>>> > > > terminate called after throwing an instance of
>>> > > > 'Dune::InterfaceBuilder::RemotexIndicesStateError'
>>> > > >
>>> > > > As far as I can see, it is called when I call the member function
>>> apply
>>> > > on
>>> > > > the RestartedGMRES solver object.
>>> > > >
>>> > >
>>> > > This only happens if the information in the RemoteIndices is not in
>>> > > sync with the underlying IndexSets. This can be case if these were
>>> > > changed, but there was not call to RemoteIndices::rebuild().
>>> > > Thi is also reflected in the error message: "RemoteIndices is not in
>>> > > sync with the index set. Call RemoteIndices::rebuild first!"
>>> > >
>>> > > So the question is: How do you set up the index sets and the remote
>>> > > information that you are using in OwnerOverlapCopyCommunication?
>>> > >
>>> > > Cheers,
>>> > >
>>> > > Markus
>>> > >
>>> > >
>>> > > Dr. Markus Blatt - HPC-Simulation-Software & Services
>>> > > http://www.dr-blatt.de
>>> > > Hans-Bunte-Str. 8-10, 69123 Heidelberg, Germany
>>> > > Tel.: +49 (0) 160 97590858
>>> > >
>>> > > _______________________________________________
>>> > > Dune mailing list
>>> > > Dune at dune-project.org
>>> > > http://lists.dune-project.org/mailman/listinfo/dune
>>> > >
>>> > >
>>>
>>> --
>>> Do you need more support with DUNE or HPC in general?
>>>
>>> Dr. Markus Blatt - HPC-Simulation-Software & Services
>>> http://www.dr-blatt.de
>>> Hans-Bunte-Str. 8-10, 69123 Heidelberg, Germany
>>> Tel.: +49 (0) 160 97590858
>>>
>>> _______________________________________________
>>> Dune mailing list
>>> Dune at dune-project.org
>>> http://lists.dune-project.org/mailman/listinfo/dune
>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20150614/9ea3459b/attachment.htm>