<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><div class="">(moving to petsc-dev) </div><div class=""><br class=""></div>To follow up further on this, Matt is correct as to what's happening now, but periodic coordinates aren't sufficiently supported yet in DMStag, so I will add something.<div class=""><br class=""></div><div class="">The way things are set up now has a conceptual elegance to it, in that to define coordinates, you use another DM which has coordinate information on it, instead of other field information. It's periodic iff the primary DM is. So there is no point on the right boundary, at 2 * pi in the 1D version of this example, because that point would be identical to the point at 0, on the left boundary.</div><div class=""><br class=""></div><div class="">The problem with the current implementation (for DMStag) is that the right boundary of the domain [0, 2*pi) is never stored. There's no way to know the width of the last cell on the right. You need that information for at least two important reasons:</div><div class="">1. to visualize the mesh, where even though the boundary point is the same point on the torus, you are plotting it on the plane and want different representations of the point on the left and right.</div><div class="">2. to use PIC methods (DMSwarm), where we need a way to determine if a particle is in the last cell.</div><div class=""><br class=""></div><div class="">Matt, Mark, Dave, et. al, it'd be very helpful to know if the following seems like a good/bad idea to you, since I assume you resolved this same issue for DMPlex + DMSwarm:</div><div class=""><br class=""></div><div class="">A tempting way to proceed here is to use the existing DMSetPeriodicity(), which allows you to specify that missing piece of information and store it in the DM. This could be called from the DMStagSetUniformCoordinatesXXX() functions, so the user wouldn't have to worry about it in that case. That also makes conceptual sense as that's the stage, after setup, in which you specify the "embedding" part of the DM. A next step would be to make DMLocalizeCoordinates() work for DMStag (and DMDA if possible, while I'm at it). </div><div class=""><br class=""></div><div class=""><br class=""></div><div class=""><br class=""><div><blockquote type="cite" class=""><div class="">Am 25.03.2021 um 00:20 schrieb Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>>:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" style="caret-color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; -webkit-text-stroke-width: 0px; text-decoration: none;" class=""><div dir="ltr" class="">On Wed, Mar 24, 2021 at 7:17 PM Jorti, Zakariae via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" class="">petsc-users@mcs.anl.gov</a>> wrote:<br class=""></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div style="overflow-wrap: break-word;" class=""><div id="gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><p class="">Hi Patrick,</p><p class=""><br class=""></p><p class="">Thanks for your responses. </p><p class="">As for the code, I was not granted permission to share it yet. So, I cannot send it to you for the moment. I apologize for that.</p><p class=""><br class=""></p><p class="">I wanted to let you know that while I was testing my code, I discovered that when the periodic boundary conditions are activated, the coordinates accessed might be incorrect on one side of the boundary. </p><p class="">Let me give you an example in cylindrical coordinates with a 3x3x3 DMStag mesh: </p><div class=""> <br class="webkit-block-placeholder"></div><p class=""><br class=""></p><p class=""><br class=""></p><p class="">PetscInt<span style="white-space: pre-wrap;" class=""> </span> startr,startphi,startz,nr,nphi,nz,d;</p><div class=""><br class="webkit-block-placeholder"></div><p class="">PetscInt<span style="white-space: pre-wrap;" class=""> </span> er,ephi,ez,<span style="font-size: 12pt;" class="">icErmphip[3];</span></p><div class=""><br class="webkit-block-placeholder"></div><p class="">DM dmCoorda, coordDA;</p><p class="">Vec coordaLocal;</p><p class="">PetscScalar ****arrCoord;</p><p class="">PetscScalar surf;</p><p class=""><br class=""></p><p class="">DMStagCreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,3,3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,1,1,DMSTAG_STENCIL_BOX,1,NULL,NULL,NULL,&coordDA);<br class=""></p><p class=""><br class=""></p><p class="">DMSetFromOptions(coordDA);</p><p class="">DMSetUp(coordDA);</p><p class=""><br class=""></p><p class="">DMStagGetCorners(coordDA,&startr,&startphi,&startz,&nr,&nphi,&nz,NULL,NULL,NULL);</p><p class=""><br class=""></p><p class="">DMGetCoordinateDM(coordDA,&dmCoorda);</p><p class="">DMGetCoordinatesLocal(coordDA,&coordaLocal);</p><p class=""> <span class="Apple-converted-space"> </span>DMStagVecGetArrayRead(dmCoorda,coordaLocal,&arrCoord);</p><p class=""><br class=""></p><p class="">for (d=0; d< 3; ++d){</p><p class="">DMStagGetLocationSlot(dmCoorda,UP_LEFT,d,&icErmphip[d]);<br class=""></p><p class="">}</p><p class=""><br class=""></p><p class="">er = 1; ez = 0;</p><p class="">for (ephi=0; ephi< 3; ++ephi){</p><p class=""><span style="font-size: 12pt;" class=""> PetscPrintf(PETSC_COMM</span><span style="font-size: 12pt;" class="">_</span><span style="font-size: 12pt;" class="">WORLD,"Phi_p(%d,%d,%d) = %E</span><span style="font-size: 12pt;" class="">\</span><span style="font-size: 12pt;" class="">n</span><span style="font-size: 12pt;" class="">",</span><span style="font-size: 12pt;" class="">er</span><span style="font-size: 12pt;" class="">,</span><span style="font-size: 12pt;" class="">ephi</span><span style="font-size: 12pt;" class="">,</span><span style="font-size: 12pt;" class="">ez</span><span style="font-size: 12pt;" class="">,(double)</span><span style="font-size: 12pt;" class="">arrCoord[ez][ephi][er][icErmphip[1]</span><span style="font-size: 12pt;" class="">);</span></p><p class="">}</p><p class=""><br class=""></p><p class="">When I execute this example, I get this output:</p><p class=""><span style="font-size: 12pt;" class="">Phi_p(1,0</span><span style="font-size: 12pt;" class="">,0) =</span><span style="font-size: 12pt;" class=""> </span><span style="font-size: 12pt;" class="">2</span><span style="font-size: 12pt;" class="">.094395</span><span style="font-size: 12pt;" class="">E+00</span></p><p class=""><span style="font-size: 12pt;" class=""></span>Phi_p(1,1,0) = 4.188790E+00<br class=""></p><p class="">Phi_p(1,2,0) = 0.000000E+00</p><p class=""><br class=""></p><p class="">Note here that the first two lines correspond to 2<span style="color: rgb(33, 33, 33); font-family: Calibri, Helvetica, sans-serif, serif, EmojiFont; font-size: 12pt;" class="">π / 3 and 4</span><span style="color: rgb(33, 33, 33); font-family: Calibri, Helvetica, sans-serif, serif, EmojiFont; font-size: 12pt;" class="">π / 3 respectively. Thus, nothing is wrong here. </span></p><p class=""><span style="color: rgb(33, 33, 33); font-family: Calibri, Helvetica, sans-serif, serif, EmojiFont; font-size: 12pt;" class="">But the last line should rather give 2</span><span style="color: rgb(33, 33, 33); font-family: Calibri, Helvetica, sans-serif, serif, EmojiFont; font-size: 12pt;" class="">π instead of 0.</span></p><p class=""><span style="color: rgb(33, 33, 33); font-family: Calibri, Helvetica, sans-serif, serif, EmojiFont; font-size: 12pt;" class=""><br class=""></span></p><p class=""><span style="font-size: 12pt;" class="">I understand that degrees of freedom should be the same on both sides of the boundary, but should the coordinates not be preserved? </span></p><div class=""><span style="font-size: 12pt;" class=""></span><br class="webkit-block-placeholder"></div></div></div></blockquote><div class="">I don't think so. The circle has coordinates in [0, 2\pi), so the point at 2\pi is identified with the point at 0 and you must choose</div><div class="">one, so we choose 0.</div><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt </div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div style="overflow-wrap: break-word;" class=""><div id="gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><p class=""><span style="font-size: 12pt;" class="">Thank you.</span></p><p class=""><span style="font-size: 12pt;" class="">Best regards,</span></p><p class=""><span style="font-size: 12pt;" class=""><br class=""></span></p><p class=""><span style="font-size: 12pt;" class="">Zakariae Jorti</span></p></div><hr style="display: inline-block; width: 561.25px;" class=""><div id="gmail-m_-5207082944899504102divRplyFwdMsg" dir="ltr" class=""><font face="Calibri, sans-serif" style="font-size: 11pt;" class=""><b class="">From:</b><span class="Apple-converted-space"> </span>Patrick Sanan <<a href="mailto:patrick.sanan@gmail.com" target="_blank" class="">patrick.sanan@gmail.com</a>><br class=""><b class="">Sent:</b><span class="Apple-converted-space"> </span>Tuesday, March 23, 2021 11:37:04 AM<br class=""><b class="">To:</b><span class="Apple-converted-space"> </span>Jorti, Zakariae<br class=""><b class="">Cc:</b><span class="Apple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov" target="_blank" class="">petsc-users@mcs.anl.gov</a><br class=""><b class="">Subject:</b><span class="Apple-converted-space"> </span>[EXTERNAL] Re: Question about periodic conditions</font><div class=""> </div></div><div class="">Hi Zakariae - sorry about the delay - responses inline below.<div class=""><br class=""></div><div class="">I'd be curious to see your code (which you can send directly to me if you don't want to post it publicly), so I can give you more comments, as DMStag is a new component. <br class=""><div class=""><br class=""></div><div class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">Am 23.03.2021 um 00:54 schrieb Jorti, Zakariae <<a href="mailto:zjorti@lanl.gov" target="_blank" class="">zjorti@lanl.gov</a>>:</div><br class=""><div class=""><div id="gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class="">Hi, </div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><br class=""></div><div style="margin-top: 0px; margin-bottom: 0px;" class="">I implemented a PETSc code to solve Maxwell's equations for the magnetic and electric fields (B and E) in a cylinder: </div><div style="margin-top: 0px; margin-bottom: 0px;" class="">0 < r_min <= r <= r_max; with <span style="font-size: 12pt;" class="">r_max</span><span style="font-size: 12pt;" class=""> > r_min</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class="">phi_min = 0 <span style="font-size: 12pt;" class=""><= r <= phi_max</span><span style="font-size: 12pt;" class=""> <span class="Apple-converted-space"> </span>= 2 </span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">π</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">z_min <= z =< z_max; </span><span style="font-size: 12pt;" class="">with </span><span style="font-size: 12pt;" class="">z_max</span><span style="font-size: 12pt;" class=""> > z_min.</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font size="3" class="">I am using a PETSc staggered grid with the electric field E defined on edge centers and the magnetic field B </font>defined<font size="3" class=""> on face centers. (</font><span style="font-size: 12pt;" class="">dof0 = 0, dof1 = 1,dof2 = 1, dof3 = 0;</span><font size="3" class="">). </font><font size="3" class=""> </font></div></div></div></blockquote><blockquote type="cite" class=""><div class=""><div id="gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-size: 12pt;" class=""></span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">I have two versions of my code: </span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font face="Calibri, Helvetica, sans-serif" size="3" class="">1 - A first version<span class=""> </span></font><font face="Calibri, Helvetica, sans-serif" class="">in which I set the boundary type to DM_BOUNDARY_NONE in the three directions r, phi and z</font></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">2- A second version in which I set </span><span style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class="">the boundary type to DM_BOUNDARY_NONE in the r and z directions</span><span style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class="">, </span><font face="Calibri, Helvetica, sans-serif" size="3" class="">and DM_BOUNDARY_PERIODIC in the phi di</font><font face="Calibri, Helvetica, sans-serif" class="">rection.</font><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""> </span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">When I print the solution vector X, which contains both E and B components, I notice that the vector is shorter with the second version compared to the first one. </span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">Is it normal? </span></div></div></div></blockquote>Yes - with the periodic boundary conditions, there will be fewer points since there won't be the "extra" layer of faces and edges at phi = 2 * pi .</div><div class=""><br class=""></div><div class="">If you consider a 1-d example with 1 dof on vertices and cells, with three elements, the periodic case looks like this, globally,</div><div class=""><br class=""></div><div class=""> <span class="Apple-converted-space"> </span>x ---- x ---- x ----</div><div class=""><br class=""></div><div class="">as opposed to the non-periodic case,</div><div class=""><br class=""></div><div class=""> x ---- x ---- x ---- x</div><div class=""><br class=""></div><div class=""><br class=""></div><div class=""><blockquote type="cite" class=""><div class=""><div id="gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">Besides, I was wondering if I have to change the way I define the value of the solution on the boundary. What I am doing so far in both versions is something like:</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">B_phi [phi = 0] = 1.0;</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif;" class="">B_phi </span><span style="font-family: Calibri, Helvetica, sans-serif;" class="">[phi = 2</span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">π</span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">] = 1.0;</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class="">E_z [r, phi = 0</span><span style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class="">] = 1/r;</span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">E_z [r, phi = 2</span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">π</span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">] = 1/r;</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""><br class=""></span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><font face="Calibri, Helvetica, sans-serif" class="">Assuming that values at phi = 0 should be the same as at phi=2</font><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">π with the periodic boundary conditions, </span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">is it sufficient for example </span><font face="Calibri, Helvetica, sans-serif" size="3" class="">to </font><font face="Calibri, Helvetica, sans-serif" class="">have</font><font face="Calibri, Helvetica, sans-serif" size="3" class=""> only the following boundary conditions</font><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">:</span></div><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class="">B_phi [phi = 0] = 1.0;</span></div><p style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""></span></p><div style="margin-top: 0px; margin-bottom: 0px;" class=""><span style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class="">E_z [r, phi = 0</span><span style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class="">] = 1/r ? </span><span style="font-family: Calibri, Helvetica, sans-serif; font-size: 12pt;" class=""><br class=""></span></div></div></div></blockquote><div class=""><br class=""></div>Yes - this is the intention, since the boundary at phi = 2 * pi is represented by the same entries in the global vector. </div><div class=""><br class=""></div><div class="">Of course, you need to make sure that your continuous problem is well-posed, which in general could change when using different boundary conditions. </div><div class=""><br class=""><blockquote type="cite" class=""><div class=""><div id="gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; text-decoration: none; font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class=""><div class=""><span style="font-size: 12pt; font-family: Calibri, Helvetica, sans-serif;" class="">Thank you.</span></div><div class=""><font face="Calibri, Helvetica, sans-serif" class="">Best regards,</font></div><div class=""><font face="Calibri, Helvetica, sans-serif" class=""><br class=""></font></div><div class=""><font face="Calibri, Helvetica, sans-serif" class="">Zakariae Jorti</font></div></div></div></blockquote></div><br class=""></div></div></div></div></blockquote></div><br clear="all" class=""><div class=""><br class=""></div>--<span class="Apple-converted-space"> </span><br class=""><div dir="ltr" class="gmail_signature"><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a></div></div></div></div></div></div></div></div></div></blockquote></div><br class=""></div></body></html>