[petsc-users] Parallel assembly FEM matrix

[Matthew Knepley]

On Tue, Jun 7, 2016 at 6:38 AM, Nestor Cerpa <ncerpagi at umn.edu> wrote:

> 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 ?
>

MatSetValuesLocal() just uses a mapping from local indices to global
indices, and then calls MatSetValues(). It does not
assemble a different type of matrix.

 Thanks,

     Matt


> 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 <cerpa003 at umn.edu>
>
>


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160607/5a1f56cf/attachment-0001.html>