[petsc-users] Parallel assembly FEM matrix

[Nestor Cerpa]

Hi Lawrence, 

Thank you. Now it works!

I’ve looked at the source code for DMPlexMatSetClosure() and it uses MatSetValues() so I am wondering in which situation should we use MatSetValuesLocal() ? and how to scatter the local matrices assembled with this ?


Nestor


> On Jun 6, 2016, at 2:38 AM, Lawrence Mitchell <lawrence.mitchell at imperial.ac.uk> wrote:
> 
> 
>> 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
> 



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

Nestor Cerpa Gilvonio
Postdoctoral research associate

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

E-mail : ncerpagi at umn.edu <mailto:cerpa003 at umn.edu>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160607/8ac1946e/attachment.html>