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 <[email protected]> wrote: > On Fri, Jan 6, 2017 at 8:52 AM, Rochan Upadhyay <[email protected]> > 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 <[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/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 <[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 >>> >> >> > > > -- > 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 >
