On 2/19/19 7:40 AM, Matthew Knepley wrote:
On Tue, Feb 19, 2019 at 1:13 AM Jordan Wagner via petsc-users
<[email protected] <mailto:[email protected]>> 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/>
