[petsc-users] a question on DMPlexSetAnchors

[Matthew Knepley]

On Fri, Jan 6, 2017 at 8:52 AM, Rochan Upadhyay <u.rochan at gmail.com> wrote:

> Constraints come from so-called cohomology conditions. In practical
> applications,
> they arise when you couple field models (e.g. Maxwell's equations) with
> lumped
> models (e.g. circuit equations). They are described in this paper :
> http://gmsh.info/doc/preprints/gmsh_homology_preprint.pdf
>

This looks interesting, but I wish they were more explicit about what was
actually being solved.


> In their matrix in page 12 all rows and columns involving the terms
> <E1,*>, <E2,*>,
> <*,E1> and,  <*,E2> are non-local. That is because the "cohomology" basis
> functions
> E1 and E2, are sums of basis functions defined on all points contained in
> a group of cells.
>

Okay, then to me it looks like you have

  M + L = / A 0 \
               \ 0  I /

where M is a sparse, block diagonal matrix (maybe you do not have the I),
and L is low-rank. You
can certainly lay this out with a Section by having 4 fields. It will
almost certainly be that the
Jacobian layout is wrong due to the non-locality, but you can turn off
checking so that you can insert
new nonzeros using MatSetOption(M, MAT_NEW_NONZERO_LOCATION_ERR,
PETSC_FALSE).
Does this make sense?

I guess this structure will kill the performance of most existing
> preconditioners but I
> would like to initially look at smallish problems.
>

Yes, it will kill performance unless we treat the matrix as M + L, which
you can do using the MatLRC
type. However, we can postpone that until you have everything working and
want to get bigger. Also,
integral constraints can sometimes be handled using fast methods for
integral equations.

  Thanks,

     Matt


> On Thu, Jan 5, 2017 at 8:40 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Thu, Jan 5, 2017 at 6:35 PM, Rochan Upadhyay <u.rochan at gmail.com>
>> wrote:
>>
>>> Thanks for prompt reply. I don't need hanging nodes or Dirichlet
>>> conditions which can
>>> be easily done by adding constraint DoFs in the Section as you mention.
>>> My requirement is the following:
>>> >>> Constraints among Fields:
>>> >>> I would recommend just putting the constraint in as an equation. In
>>> your case the effect can
>>> >>> be non-local, so this seems like the best strategy.
>>> The constraint dof is described by an equation. In fact I have easily
>>> set up residuals for the system. My (perceived) difficulties are in the
>>> Jacobian. My additional
>>> Dof is a scalar quantity that is not physically tied to any specific
>>> point but needs to be solved tightly coupled
>>> to a FEM system. In order to use the global section (default section for
>>> the FEM system)
>>> to fill up the Mats and Vecs, I have artificially appended this extra
>>> dof to a particular point.
>>> Now in the Jacobian matrix there will be one extra row and column that,
>>> once filled, should be dense
>>> (rather block dense) due to the non-local dependence of this extra Dof
>>> on field values at some other points.
>>>
>>
>> Now, if you want good performance, you have to describe the constraint in
>> terms of the topology. All our DMs
>> are setup for local equations. Nonlocal equations are not correctly
>> preallocated. You can
>>
>>   a) Just turn off checking for proper preallocation, MatSetOption(A,
>> MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)
>>
>>   b) Do the preallocation yourself
>>
>> If instead, the pattern "fits inside" a common pattern described by these
>>
>>   http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpage
>> s/DM/DMPlexGetAdjacencyUseClosure.html
>>   http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpage
>> s/DM/DMPlexSetAdjacencyUseCone.html
>>
>> you can just use that.
>>
>> What creates your constraints?
>>
>>    Matt
>>
>> My question is once the DM has allocated non-zeros for the matrix (based
>>> on the given section) would it be
>>> possible to add non-zeros in non-standard locations (namely a few dense
>>> sub-rows and sub-columns) in a way
>>> that does not destroy performance. Does using the built in routine
>>> DMSetDefaultConstraint (or for that
>>> matter the DMPlexSetAnchors) create another (separate) constraint matrix
>>> that presumably does an efficient job
>>> of incorporating these additional non-zeros ? Or does this Constraint
>>> matrix only come in during the DMLocalToGLobal
>>> (& vice versa) calls as mentioned in the documentation ?
>>> I appreciate your reading through my rather verbose mail, especially
>>> considering the numerous other queries that
>>> you receive each day.
>>> Thanks.
>>>
>>> On Wed, Jan 4, 2017 at 5:59 PM, Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>>> On Tue, Jan 3, 2017 at 4:02 PM, Rochan Upadhyay <u.rochan at gmail.com>
>>>> wrote:
>>>>
>>>>> I think I sent my previous question (on Dec 28th) to the wrong place
>>>>> (petsc-users-request at mcs.anl.gov).
>>>>>
>>>>
>>>> Yes, this is the correct mailing list.
>>>>
>>>>
>>>>> To repeat,
>>>>>
>>>>> I am having bit of a difficulty in understanding the introduction of
>>>>> constraints in DMPlex. From a quick study of the User Manual I gather
>>>>> that it is easiest done using DMPlexSetAnchors ? The description of
>>>>> this
>>>>> routine says that there is an anchorIS that specifies the anchor
>>>>> points (rows in the
>>>>> matrix). This is okay and easily understood.
>>>>>
>>>>
>>>> I think this is not the right mechanism for you.
>>>>
>>>> Anchors:
>>>>
>>>> This is intended for constraints in the discretization, such as hanging
>>>> nodes, which are
>>>> purely local, and intended to take place across the entire domain. That
>>>> determines the
>>>> interface.
>>>>
>>>> Dirichlet Boundary Conditions:
>>>>
>>>> For these, I would recommend using the Constraint interface in
>>>> PetscSection, which
>>>> eliminates these unknowns from the global system, but includes the
>>>> values in the local
>>>> vectors used in assembly.
>>>>
>>>> You can also just alter your equations for constrained unknowns.
>>>>
>>>> Constraints among Fields:
>>>>
>>>> I would recommend just putting the constraint in as an equation. In
>>>> your case the effect can
>>>> be non-local, so this seems like the best strategy.
>>>>
>>>>   Thanks,
>>>>
>>>>      Matt
>>>>
>>>>
>>>>> There is also an anchorSection which is described as a map from
>>>>> constraint points
>>>>> (columns ?)  to the anchor points listed in the anchorIS. Should this
>>>>> not be a map between
>>>>> solution indices (i.e. indices appearing in the vectors and matrices) ?
>>>>>
>>>>> For example I am completely unable to set up a simple constraint
>>>>> matrix for the following (say):
>>>>>
>>>>> Point 1, Field A, B
>>>>> Point 2-10 Field A
>>>>> At point 1, Field B depends on Field A at points 1-10
>>>>>
>>>>> When I set it up it appears to create a matrix where field A depends
>>>>> on field A values at points 1-10.
>>>>>
>>>>> How does the mapping work in this case ? Will the DMPlexSetAnchors()
>>>>> routine work
>>>>> for this simple scenario ?
>>>>>
>>>>> If not, is the only recourse to create the constraint matrix oneself
>>>>> using DMSetDefaultConstraints ?
>>>>>
>>>>> Also documentation for DMSetDefaultConstraints is incomplete.
>>>>> The function accepts three arguments (dm, section and Mat) but
>>>>> what the section is is not described at all.
>>>>>
>>>>> I don't know if my question makes any sense. If it does not then it is
>>>>> only a reflection of my utter confusion regarding the routine
>>>>> DMPlexSetAnchors :-(
>>>>>
>>>>> Regards,
>>>>> Rochan
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> 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
>>>>
>>>
>>>
>>
>>
>> --
>> 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
>>
>
>


-- 
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/20170106/31a32e51/attachment.html>