[petsc-users] Parallel assembly FEM matrix

[Lawrence Mitchell]

> On 6 Jun 2016, at 01:43, Nestor Cerpa <ncerpagi at umn.edu> wrote:
> 
> Hello,
> 
> I am learning to use Petsc in order to modify my sequential FEM code to work in parallel. First, I am trying to understand how a FEM assembly loop would work in parallel using Petsc but I run into some problems.
> 
> I made a simple example to explain the issue :
> - I read a simple mesh (generated by gmsh) with 4 quadrilaterls and generate the corresponding DMPlex
> - I define a scalar field on vertices
> - Vertices on Bottom and Top boundaries are contrained (so I have 3 unknowns)
> - I create a DM Matrix
> - For the purpose of my example during the assembly I define an “elementary matrix” containing only ones
> 
> When I run the program with one process I obtain what the matrix is supposed to be :
> 
> Mat Object: 1 MPI processes
>   type: seqaij
> row 0: (0, 2.)  (2, 2.)
> row 1: (1, 2.)  (2, 2.)
> row 2: (0, 2.)  (1, 2.)  (2, 4.)
> 
> But running with two processes gives me this :
> 
> Mat Object: 2 MPI processes
>   type: mpiaij
> row 0: (0, 2.)  (2, 0.)
> row 1: (1, 2.)  (2, 2.)
> row 2: (0, 0.)  (1, 2.)  (2, 2.)
> 
> It seems that in the way I am doing the assembly, Process 0 is not giving the contribution of its local assembly to process 1.
> 
> 
> For the assembly I am doing something like this (the complete example is attached) :
> 
> 
>  !*************************************
>  !--- Matrix Assembly
>  !*************************************
> 
>  call DMCreateMatrix(dm,A,ierr)
>  call DMGetLocalToGlobalMapping(dm,ltog,ierr)
>  call MatSetLocalToGlobalMapping(A,ltog,ltog,ierr);
>  call ISLocalToGlobalMappingDestroy(ltog,ierr);
> 
> 
> do icell = !--- Loop over cells
> 
>            do ic = ! Loop over closure (somehow only vertices)
>                ipoin = closure((ic-1)*2+1)
>                call DMPlexGetPointLocal(dm,ipoin,pStart,pEnd,ierr)
>                !
>                !  Something to get location of point data in a local Vec
>                !
>           enddo
> 
>           elstf(:,:) = 1.0
> 
>          do i = !Loop over rows corresponding to dofs which are not constrained
>                call MatSetValuesLocal(A,1,local_ldof(i)-1,ndof*nodes,local_ldof(:)-1, elstf(i,:),ADD_VALUES,ierr)
>          enddo
> 
> enddo
> 
>   call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
>   call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
> 
>  !*************************************
>  !--- End Assembly
>  !*************************************
> 
> 
> 
> There is something I am definitely missing because I don’t understand why it does not work. I was thinking maybe I need a function like DMLocallToGlobal…() but for matrices.


You can either do this "by hand" using DMPlexGetTransitiveClosure to get the points. Or else you use DMPlexMatSetClosure to scatter the element matrix to the appropriate points in the matrix.

Note there is also a matching DMPlexVecSetClosure to set values in the residual vector, and DMPlexVecGetClosure to get the values in the closure of a point.

So in pseudo-code you'd probably do:

for c in cells:
   DMPlexVecGetClosure(dm, NULL, coeff, c, &csize, &values)
   ! build element matrix (using values)
   ...
   DMPlexVecRestoreClosure(dm, NULL, coeff, c, &size, &values)
   DMPlexMatSetClosure(dm, NULL, NULL, mat, c, element_matrix, ADD_VALUES)

Cheers,

Lawrence

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160606/040924cd/attachment.pgp>