[petsc-users] a question on DMPlexSetAnchors

[Rochan Upadhyay]

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

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

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/manualpages/DM/
> DMPlexGetAdjacencyUseClosure.html
>   http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170106/2828640c/attachment.html>