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

Matthew Knepley knepley at gmail.com
Sun Feb 24 16:50:22 CST 2019 

On Sun, Feb 24, 2019 at 4:28 PM Appel, Thibaut <t.appel17 at imperial.ac.uk>
wrote:

> Hi Matthew,
>
> With the following data declaration and routines calls:
>

I am not a Fortran90 expert, but your declaration looks different than the
ones I have seen work:


https://bitbucket.org/petsc/petsc/src/master/src/dm/impls/plex/examples/tutorials/ex1f90.F90

  Thanks,

     Matt


>
>   DM :: da
>   ISLocalToGlobalMapping :: ltog
>   PetscInt, DIMENSION(:), POINTER :: id_ltog
>
>   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)
>
>   CALL DMGetLocalToGlobalMapping(da,ltog,ierr)
>   CHKERRA(ierr)
>
>   CALL ISLocalToGlobalMappingGetIndicesF90(ltog,id_ltog,ierr)
>   CHKERRA(ierr)
>
>   CALL ISLocalToGlobalMappingRestoreIndicesF90(ltog,id_ltog,ierr)
>   CHKERRA(ierr)
>
>   CALL DMDestroy(da,ierr)
>   CHKERRA(ierr)
>
>
> I get, with the most recent Intel Fortran compiler and a fresh “git pull”
> PETsc:
>
> forrtl: severe (408): fort: (7): Attempt to use pointer ID_LTOG when it is
> not associated with a target
>
>
> Thibaut
>
> On 22 Feb 2019, at 15:13, Matthew Knepley <knepley at gmail.com> wrote:
>
> 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/
>
>
>
> --
> 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/
>
>
>

-- 
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/20190224/48611fb8/attachment.html>

More information about the petsc-users mailing list