[Dune] Runtime error with parallel linear solver

[Tobias Ritschel]

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/20150613/c8151782/attachment.htm>