[petsc-users] Including constrained dofs in a matrix/vector created with DMCreateXXXX()

Wed Feb 20 14:38:32 CST 2019

On 2/19/19 7:40 AM, Matthew Knepley wrote:
> On Tue, Feb 19, 2019 at 1:13 AM Jordan Wagner via petsc-users 
> <petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>> wrote:
>
>     Hi,
>
>     Over the past few months, I have implemented dmplex/section in my
>     preexisting finite element code. We already had our own
>     implementations for things like the finite element space, boundary
>     conditions, projections, etc, so I am mostly using dmplex/section
>     at its lowest level; no DT/DS stuff for the moment. Everything is
>     working fine at this level, and I am overall very happy with how
>     it ties into the code that we already have.
>
>     I am currently doing my own preallocation for the system Jacobian
>     but am wanting to use DMCreateMatrix() instead, as it seems
>     pointless to keep using my own crappy function when this one is
>     just sitting there waiting. However, I dislike how I am currently
>     implementing this and was hoping to get some suggestions or some
>     advice. My problem is handling nicely the Dirichlet dofs. My code
>     takes the approach of first assembling the entire Jacobian matrix
>     and load vector as if no constraints are imposed. We then have a
>     function that applies the essential conditions after assembly and
>     extracts a submatrix and subvector similar to the way ex49
>     <https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex49.c.html>
>     is doing.
>
>     Since DMCreateMatrix() automatically partitions out the
>     constraints specified in the section, I find myself having to
>     create two nearly identical sections: one that has constraints and
>     one that does not, and then I clone the DM and set the default
>     section of each respectively. I can then use DMCreateMatrix() on
>     the one dm/section with no constraints to preallocate the larger
>     matrix that I am assembling. Then I use the dm/section with
>     constraints to do all of my boundary condition and projection
>     operations on the previously created matrix and vector.
>     Essentially, I am using an unconstrained section only for
>     allocation and a constrained section for most everything else.
>
>     This process actually works fine for me, but its pretty ugly, and
>     I am sure there is a better way. I am wondering if I am missing
>     something that could keep me from having to go through this
>     process of creating two sections that differ only in the
>     constraints. It seems to me if there were an option for
>     DMCreateMatrix()/DMCreateVector() to either include or not include
>     the dofs, that would completely solve my problem. That way, I can
>     use the same section for both creation and allocation.
>
>     So, the main question I have is just whether there is a function
>     or option that I am missing like DMCreateXXXX() but with the
>     ability to retain the constrained rows and columns rather than
>     compressing them out.
>
> I do not completely understand the above, but here is how Section 
> works. A _local_ Section contains all dofs and is intended for 
> assembly using local Vecs. To communicate with the solver, you make a 
> _global_ Section using PetscSectionCreateGlobalSection(). This has an 
> option to include or exclude constraints. For example, I make a global 
> Section with constraints when I output for visualization. Maybe you 
> can use that to do what you want.
>
Thanks for the prompt reply and sorry for the unclear initial question. 
I think this helps in my understanding of how this tool was intended to 
be used.

> What I do not understand is why you do not just use a global Section 
> without constraints to assemble your matrix, rather than eliminating 
> the rows/columns later. The global Section gives back negative indices 
> for unknown that are constrained or not owned. These are ignored by 
> MatSetValues() and VecSetValues(). This is how the default Plex 
> assembly makes operators without constrained variables.
>
I have rarely used a global section directly because I have never been 
able to  make sense of exactly what it was giving me. For instance as 
you mentioned (and I have seen written elsewhere), the global section 
should give negative indices for unowned or constrained dofs. This seems 
to work fine for me with the unowned points---I clearly get a negative 
size and offset from the section on these points. But my problem is for 
points that have constraints. I don't see where I would get negative 
indices from. The sizes and offsets for these points are not negative.  
Of course, from the constraint size and constraint indices on these 
points, I can manually find out which should be negative and go from 
there. However, that leads me to the next problem I am having with 
global section...

For some reason, when I am creating the global section from the local 
one, my constrained component indices are getting messed up. The sizes 
are right but the order seems to be random and/or nonsensical. For 
instance, here is a snippet from my code being run in serial

Local section:

   (17615) dim  2 offset 794 constrained 0
   (17616) dim  2 offset 796 constrained 0 1
   (17617) dim  2 offset 798 constrained 0 1
   (17618) dim  2 offset 800 constrained 0 1
   (17619) dim  2 offset 802 constrained 0 1
   (17620) dim  2 offset 804 constrained 0 1
   (17621) dim  2 offset 806 constrained 0 1

Global section:

   (17615) dim  2 offset  80 constrained 1
   (17616) dim  2 offset  81 constrained 0 0
   (17617) dim  2 offset  81 constrained 0 0
   (17618) dim  2 offset  81 constrained 0 0
   (17619) dim  2 offset  81 constrained 0 0
   (17620) dim  2 offset  81 constrained 0 0
   (17621) dim  2 offset  81 constrained 0 1

So using the global section, I am unable to even manually tell which 
components are constrained. Is this normal for the global section or is 
it something going wrong in its creation? or am I completely 
misunderstanding something?

Thanks very much.

>   Thanks,
>
>      Matt
>
>     Thanks a lot, as always.
>
>
>
> -- 
> 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/20190220/b238b4d7/attachment.html>