<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Jan 6, 2017 at 8:52 AM, Rochan Upadhyay <span dir="ltr"><<a href="mailto:u.rochan@gmail.com" target="_blank">u.rochan@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Constraints come from so-called cohomology conditions. In practical applications,<br><div class="gmail_extra">they arise when you couple field models (e.g. Maxwell's equations) with lumped<br></div><div class="gmail_extra">models (e.g. circuit equations). They are described in this paper :<br><a href="http://gmsh.info/doc/preprints/gmsh_homology_preprint.pdf" target="_blank">http://gmsh.info/doc/<wbr>preprints/gmsh_homology_<wbr>preprint.pdf</a><br></div><div class="gmail_extra"></div></div></blockquote><div><br></div><div>This looks interesting, but I wish they were more explicit about what was actually being solved.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra">In their matrix in page 12 all rows and columns involving the terms <E1,*>, <E2,*>,<br><*,E1> and,  <*,E2> are non-local. That is because the "cohomology" basis functions <br>E1 and E2, are sums of basis functions defined on all points contained in a group of cells.</div></div></blockquote><div><br></div><div>Okay, then to me it looks like you have</div><div><br></div><div>  M + L = / A 0 \</div><div>               \ 0  I /</div><div><br></div><div>where M is a sparse, block diagonal matrix (maybe you do not have the I), and L is low-rank. You</div><div>can certainly lay this out with a Section by having 4 fields. It will almost certainly be that the</div><div>Jacobian layout is wrong due to the non-locality, but you can turn off checking so that you can insert</div><div>new nonzeros using MatSetOption(M, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE).<br></div><div>Does this make sense?</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra">I guess this structure will kill the performance of most existing preconditioners but I <br></div><div class="gmail_extra">would like to initially look at smallish problems. </div></div></blockquote><div><br></div><div>Yes, it will kill performance unless we treat the matrix as M + L, which you can do using the MatLRC</div><div>type. However, we can postpone that until you have everything working and want to get bigger. Also,</div><div>integral constraints can sometimes be handled using fast methods for integral equations.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Jan 5, 2017 at 8:40 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span class="gmail-m_-7161042593548520915gmail-">On Thu, Jan 5, 2017 at 6:35 PM, Rochan Upadhyay <span dir="ltr"><<a href="mailto:u.rochan@gmail.com" target="_blank">u.rochan@gmail.com</a>></span> wrote:<br><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div dir="ltr"><div><div>Thanks for prompt reply. I don't need hanging nodes or Dirichlet conditions which can<br></div></div>be easily done by adding constraint DoFs in the Section as you mention. <br><div>My requirement is the following:<br><div>>>> Constraints among Fields:</div><div>>>> I would recommend just putting the constraint in as an equation. In your case the effect can</div><div>>>> be non-local, so this seems like the best strategy.</div><div>The constraint dof is described by an equation. In fact I have easily<br></div><div>set up residuals for the system. My (perceived) difficulties are in the Jacobian. My additional <br></div><div>Dof is a scalar quantity that is not physically tied to any specific point but needs to be solved tightly coupled<br>to a FEM system. In order to use the global section (default section for the FEM system)<br></div><div>to fill up the Mats and Vecs, I have artificially appended this extra dof to a particular point.<br></div><div>Now in the Jacobian matrix there will be one extra row and column that, once filled, should be dense <br>(rather block dense) due to the non-local dependence of this extra Dof on field values at some other points. <br></div></div></div></blockquote><div><br></div></span><div>Now, if you want good performance, you have to describe the constraint in terms of the topology. All our DMs</div><div>are setup for local equations. Nonlocal equations are not correctly preallocated. You can</div><div><br></div><div>  a) Just turn off checking for proper preallocation, MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)<br></div><div><br></div><div>  b) Do the preallocation yourself</div><div><br></div><div>If instead, the pattern "fits inside" a common pattern described by these</div><div><br></div><div>  <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMPlexGetAdjacencyUseClosure.html" target="_blank">http://www.mcs.anl.gov/petsc<wbr>/petsc-current/docs/manualpage<wbr>s/DM/DMPlexGetAdjacencyUseClos<wbr>ure.html</a></div><div>  <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMPlexSetAdjacencyUseCone.html" target="_blank">http://www.mcs.anl.gov/petsc<wbr>/petsc-current/docs/manualpage<wbr>s/DM/DMPlexSetAdjacencyUseCone<wbr>.html</a></div><div><br></div><div>you can just use that.</div><div><br></div><div>What creates your constraints?</div><div><br></div><div>   Matt</div><div><div class="gmail-m_-7161042593548520915gmail-h5"><div><br></div><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div dir="ltr"><div><div>My question is once the DM has allocated non-zeros for the matrix (based on the given section) would it be<br>possible to add non-zeros in non-standard locations (namely a few dense sub-rows and sub-columns) in a way<br></div><div>that does not destroy performance. Does using the built in routine DMSetDefaultConstraint (or for that<br>matter the DMPlexSetAnchors) create another (separate) constraint matrix that presumably does an efficient job<br>of incorporating these additional non-zeros ? Or does this Constraint matrix only come in during the DMLocalToGLobal<br></div><div>(& vice versa) calls as mentioned in the documentation ?  <br></div>I appreciate your reading through my rather verbose mail, especially considering the numerous other queries that<br><div>you receive each day.<br></div><div>Thanks.   <br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Jan 4, 2017 at 5:59 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span>On Tue, Jan 3, 2017 at 4:02 PM, Rochan Upadhyay <span dir="ltr"><<a href="mailto:u.rochan@gmail.com" target="_blank">u.rochan@gmail.com</a>></span> wrote:<br><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div dir="ltr">I think I sent my previous question (on Dec 28th) to the wrong place <br>(<a href="mailto:petsc-users-request@mcs.anl.gov" target="_blank">petsc-users-request@mcs.anl.g<wbr>ov</a>).</div></blockquote><div><br></div></span><div>Yes, this is the correct mailing list.</div><span><div> </div><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div dir="ltr"><div>To repeat,<br><br>I am having bit of a difficulty in understanding the introduction of</div><div>constraints in DMPlex. From a quick study of the User Manual I gather <br>that it is easiest done using DMPlexSetAnchors ? The description of this <br>routine says that there is an anchorIS that specifies the anchor points (rows in the </div><div>matrix). This is okay and easily understood.</div></div></blockquote><div><br></div></span><div>I think this is not the right mechanism for you.</div><div><br></div><div>Anchors:</div><div><br></div><div>This is intended for constraints in the discretization, such as hanging nodes, which are</div><div>purely local, and intended to take place across the entire domain. That determines the</div><div>interface.</div><div><br></div><div>Dirichlet Boundary Conditions:</div><div><br></div><div>For these, I would recommend using the Constraint interface in PetscSection, which</div><div>eliminates these unknowns from the global system, but includes the values in the local</div><div>vectors used in assembly.</div><div><br></div><div>You can also just alter your equations for constrained unknowns.</div><div><br></div><div>Constraints among Fields:</div><div><br></div><div>I would recommend just putting the constraint in as an equation. In your case the effect can</div><div>be non-local, so this seems like the best strategy.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><span><div> </div><blockquote style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex" class="gmail_quote"><div dir="ltr"><div>There is also an anchorSection which is described as a map from constraint points </div><div>(columns ?)  to the anchor points listed in the anchorIS. Should this not be a map between <br></div><div>solution indices (i.e. indices appearing in the vectors and matrices) ?<br></div><div><br></div><div>For example I am completely unable to set up a simple constraint matrix for the following (say):<br></div><div><br>Point 1, Field A, B<br></div><div>Point 2-10 Field A<br></div><div>At point 1, Field B depends on Field A at points 1-10 <br><br></div><div>When I set it up it appears to create a matrix where field A depends on field A values at points 1-10.<br></div><div><br></div><div>How does the mapping work in this case ? Will the DMPlexSetAnchors() routine work <br>for this simple scenario ? <br> </div><div>If not, is the only recourse to create the constraint matrix oneself </div><div>using DMSetDefaultConstraints ?<br> <br></div><div>Also documentation for DMSetDefaultConstraints is incomplete.<br></div><div>The function accepts three arguments (dm, section and Mat) but<br></div><div>what the section is is not described at all. <br></div><div><br></div><div>I don't know if my question makes any sense. If it does not then it is</div><div>only a reflection of my utter confusion regarding the routine DMPlexSetAnchors :-(</div><div><br></div><div>Regards,<br></div>Rochan</div>
</blockquote></span></div><span class="gmail-m_-7161042593548520915gmail-m_6249290054596694524gmail-m_1765618532389985798HOEnZb"><font color="#888888"><br><br clear="all"><span class="gmail-HOEnZb"><font color="#888888"><span class="gmail-m_-7161042593548520915gmail-m_6249290054596694524gmail-HOEnZb"><font color="#888888"><div><br></div>-- <br><div class="gmail-m_-7161042593548520915gmail-m_6249290054596694524gmail-m_1765618532389985798m_4890121365330834689gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</font></span></font></span></font></span></div></div><span class="gmail-HOEnZb"><font color="#888888">
</font></span></blockquote></div><span class="gmail-HOEnZb"><font color="#888888"><br></font></span></div><span class="gmail-HOEnZb"><font color="#888888">
</font></span></blockquote></div></div></div><span class="gmail-HOEnZb"><font color="#888888"><div><div class="gmail-m_-7161042593548520915gmail-h5"><br><br clear="all"><div><br></div>-- <br><div class="gmail-m_-7161042593548520915gmail-m_6249290054596694524gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div></font></span></div></div>
</blockquote></div><br></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>