[petsc-users] Parallel assembly FEM matrix

[Nestor Cerpa]

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. 
I hope someone can help me with this and I am sorry if there is an obvious answer…
Thanks in advance. 

Nestor






---------------------------------------------------
---------------------------------------------------

Nestor Cerpa Gilvonio
Postdoctoral researcher

Department of Earth Sciences
University of Minnesota
310 Pillsbury Drive SE
Minnepolis, MN 55455

E-mail : ncerpagi at umn.edu

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160605/ffcc3c01/attachment-0003.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ParAssembly.F90
Type: application/octet-stream
Size: 8108 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160605/ffcc3c01/attachment-0002.obj>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160605/ffcc3c01/attachment-0004.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Q1_4cells.msh
Type: application/octet-stream
Size: 441 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160605/ffcc3c01/attachment-0003.obj>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160605/ffcc3c01/attachment-0005.html>