[petsc-users] a question on DMPlexSetAnchors

Rochan Upadhyay u.rochan at gmail.com
Fri Jan 6 12:22:29 CST 2017 

Yes, the option MatSetOption(M, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)
seems to be the path of least resistance. Especially as it is something I
am doing out of my
own curiosity and not part of anything larger.
I might have to bug you again very soon on how to optimize or move forward
based on how
things turn out after a first run.
Thanks.
Regards,
Rochan


On Fri, Jan 6, 2017 at 10:15 AM, Matthew Knepley <knepley at gmail.com> wrote:

> 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/0b108a1f/attachment.html>

More information about the petsc-users mailing list