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

Fri Feb 22 09:13:49 CST 2019

On Fri, Feb 22, 2019 at 9:10 AM Thibaut Appel via petsc-users <
petsc-users at mcs.anl.gov> wrote:

> 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
>
1) There was a bug in that file, but I don't think it should have affected
you. I will put in a PR

2)  You are supposed to pass in a pointer, unallocated. We manage all the
allocation, which is why you call Restore. What happens when you do that?

  Thanks,

     Matt

> [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> 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/>
>
>

-- 
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/73218c20/attachment.html>