[petsc-users] a question on DMPlexSetAnchors

[Matthew Knepley]

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/20170105/79e36096/attachment-0001.html>