[petsc-users] About DMDA (and extracting its ordering)

[Thibaut Appel]

I reckon that what corresponds best to my situation is to create my own 
MPIAIJ matrix, recover the DMDA local to global mapping, preallocate 
myself and fill my values with MatSetValues.

However, what is the correct way to call 
ISLocalToGlobalMappingGetIndices in Fortran? It seems there's something 
missing in the documentation as my code just won't work.

Some examples include a PetscOffset argument which is not present in the 
documentation and there is no manual page for 
ISLocalToGlobalMappingGetIndicesF90: the only information I could find 
is on /src/vec/is/utils/f90-custom/zisltogf90.c. Are you supposed to 
allocate the index array yourself? I tried declaring it as a pointer 
too. I keep getting an error

[0]PETSC ERROR: #1 User provided function() line 0 in User file
application called MPI_Abort(MPI_COMM_SELF, -42924) - process 0

Minimal example:

PROGRAM test_dmda

#include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petscdmda.h>

   USE PetscDMDA

   IMPLICIT NONE

   PetscErrorCode :: ierr
   PetscInt       :: irank, icx, icy, lx, ly, n_ltog
   DM             :: da
   ISLocalToGlobalMapping :: ltog

   PetscInt, DIMENSION(:), ALLOCATABLE :: id_ltog
   !PetscInt, DIMENSION(:), POINTER :: id_ltog

   PetscInt, PARAMETER :: nx=24, ny=9, neq=4
   CHARACTER(LEN=100) :: err_msg

   CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)

   CALL MPI_COMM_RANK(PETSC_COMM_WORLD,irank,ierr)

   CALL 
DMDACreate2D(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,(nx+1),(ny+1), 
&
PETSC_DECIDE,PETSC_DECIDE,neq,0,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
   CHKERRA(ierr)

   CALL DMSetUp(da,ierr)
   CHKERRA(ierr)

   CALL 
DMDAGetCorners(da,icx,icy,PETSC_NULL_INTEGER,lx,ly,PETSC_NULL_INTEGER,ierr)
   CHKERRA(ierr)

   WRITE(*,'(*(1X,A,I3))') 'I am process #', irank, ' - my lower left 
corner is icx=', icx, &
                           'and icy=', icy, ' - the width is lx=', lx, ' 
and ly=', ly

   CALL DMGetLocalToGlobalMapping(da,ltog,ierr)
   CHKERRA(ierr)

   CALL ISLocalToGlobalMappingGetSize(ltog,n_ltog,ierr)
   CHKERRA(ierr)

   ALLOCATE(id_ltog(1:n_ltog),STAT=ierr,ERRMSG=err_msg)

   CALL ISLocalToGlobalMappingGetIndices(ltog,id_ltog,ierr)
   CHKERRA(ierr)

   CALL ISLocalToGlobalMappingRestoreIndices(ltog,id_ltog,ierr)
   CHKERRA(ierr)

   DEALLOCATE(id_ltog,STAT=ierr,ERRMSG=err_msg)

   CALL DMDestroy(da,ierr)
   CHKERRA(ierr)

   CALL PetscFinalize(ierr)

END PROGRAM test_dmda


On 21/02/2019 17:49, Matthew Knepley wrote:
> On Thu, Feb 21, 2019 at 11:16 AM Thibaut Appel via petsc-users 
> <petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>> wrote:
>
>     Dear PETSc developers/users,
>
>     I’m solving linear PDEs on a regular grid with high-order finite
>     differences, assembling an MPIAIJ matrix to solve linear systems
>     or eigenvalue problems. I’ve been using vertex major, natural
>     ordering for the parallelism with PetscSplitOwnership (yielding
>     rectangular slices of the physical domain) and wanted to move to
>     DMDA to have a more square-ish domain decomposition and minimize
>     communication between processes.
>
>     However, my application is memory critical, and I have
>     finely-tuned matrix preallocation routines for allocating memory
>     “optimally”. It seems the memory of a DMDA matrix is allocated
>     along the value of the stencil width of DMDACreate and the manual
>     says about it
>
>     “These DMDA stencils have nothing directly to do with any finite
>     difference stencils one might chose to use for a discretization”
>
>     And despite reading the manual pages there must be something I do
>     not understand in the DM topology, what is that "stencil width"
>     for then? I will not use ghost values for my FD-method, right?
>
>
> What this is saying is, "You might be using some stencil that is not 
> STAR or BOX, but we are preallocating according to one of those".
> If you really care about how much memory is preallocated, which it 
> seems you do, then you might be able to use
>
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMDA/DMDASetBlockFills.html
>
> to tell use exactly how to preallocate.
>
>     I was then wondering if I could just create a MPIAIJ matrix, and
>     with a PETSc routine get the global indices of the domain for each
>     process: in other words, an equivalent of PetscSplitOwnership that
>     gives me the DMDA unknown ordering. So I can feed and loop on that
>     in my preallocation and assembly routines.
>
>
> You can make an MPIAIJ matrix yourself of course. It should have the 
> same division of rows as the DMDA division of dofs. Also, 
> MatSetValuesStencil() will not work for a custom matrix.
>
>   Thanks,
>
>      Matt
>
>     Thanks very much,
>
>     Thibaut
>
>
>
> -- 
> What most experimenters take for granted before they begin their 
> experiments is infinitely more interesting than any results to which 
> their experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/ 
> <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190222/3b60dc7d/attachment.html>