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

[Matthew Knepley]

On Wed, Feb 20, 2019 at 3:38 PM Jordan Wagner <j.wagner at rice.edu> wrote:

> 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> 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...
>
> You are right. They do not give negative indices (I was just spacing out),
they do not give indices at all. The sizes in
the global section for those points are smaller than the local 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
>
> That is a documentation failure. I did not copy over the constrained
indices into the global section because they
do not really mean anything there. In the local section, you can find the
dofs that get constrained, but in the global
section they are eliminated, so the numbering would not make sense.
Unfortunately, I did not think of a tricky way to
avoid printing them out. I guess I should set negative values and have the
printing check that. So as you say, the sizes
are right, and the constrained indices are meaningless.

  Thanks,

    Matt

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

-- 
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/87f8005c/attachment-0001.html>