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 <[email protected]> wrote: > On Thu, Jan 5, 2017 at 6:35 PM, Rochan Upadhyay <[email protected]> > 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 <[email protected]> >> wrote: >> >>> On Tue, Jan 3, 2017 at 4:02 PM, Rochan Upadhyay <[email protected]> >>> wrote: >>> >>>> I think I sent my previous question (on Dec 28th) to the wrong place >>>> ([email protected]). >>>> >>> >>> 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 >
