<div dir="ltr"><div dir="ltr">On Wed, Mar 31, 2021 at 3:21 AM Patrick Sanan <<a href="mailto:patrick.sanan@gmail.com">patrick.sanan@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="overflow-wrap: break-word;"><div>(moving to petsc-dev) </div><div><br></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><br></div><div>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><br></div><div>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>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>2.  to use PIC methods (DMSwarm), where we need a way to determine if a particle is in the last cell.</div><div><br></div><div>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><br></div><div>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></blockquote><div><br></div><div>That is what I added that stuff for. In Plex, in order to generalize to situations not in a box, we went to a formulation that uses a DG coordinate field instead. I think that</div><div>is overkill here and would not give you any added functionality.</div><div><br></div><div>Jed may argue that he wants you to retain the far point and use L2G to eliminate it, but that sounds like a lot more work.</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 style="overflow-wrap: break-word;"><div><div><blockquote type="cite"><div>Am 25.03.2021 um 00:20 schrieb Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>:</div><br><div><div dir="ltr" style="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;text-decoration:none"><div dir="ltr">On Wed, Mar 24, 2021 at 7:17 PM Jorti, Zakariae via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div id="gmail-m_3112071783902229590gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif"><p>Hi Patrick,</p><p><br></p><p>Thanks for your responses. </p><p>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><br></p><p>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>Let me give you an example in cylindrical coordinates with a 3x3x3 DMStag mesh:  </p><div> <br></div><p><br></p><p><br></p><p>PetscInt<span style="white-space:pre-wrap"> </span>     startr,startphi,startz,nr,nphi,nz,d;</p><div><br></div><p>PetscInt<span style="white-space:pre-wrap"> </span>     er,ephi,ez,<span style="font-size:12pt">icErmphip[3];</span></p><div><br></div><p>DM                dmCoorda, coordDA;</p><p>Vec               coordaLocal;</p><p>PetscScalar       ****arrCoord;</p><p>PetscScalar       surf;</p><p><br></p><p>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></p><p><br></p><p>DMSetFromOptions(coordDA);</p><p>DMSetUp(coordDA);</p><p><br></p><p>DMStagGetCorners(coordDA,&startr,&startphi,&startz,&nr,&nphi,&nz,NULL,NULL,NULL);</p><p><br></p><p>DMGetCoordinateDM(coordDA,&dmCoorda);</p><p>DMGetCoordinatesLocal(coordDA,&coordaLocal);</p><p> <span> </span>DMStagVecGetArrayRead(dmCoorda,coordaLocal,&arrCoord);</p><p><br></p><p>for (d=0; d< 3; ++d){</p><p>DMStagGetLocationSlot(dmCoorda,UP_LEFT,d,&icErmphip[d]);<br></p><p>}</p><p><br></p><p>er = 1; ez = 0;</p><p>for (ephi=0; ephi< 3; ++ephi){</p><p><span style="font-size:12pt"> PetscPrintf(PETSC_COMM</span><span style="font-size:12pt">_</span><span style="font-size:12pt">WORLD,"Phi_p(%d,%d,%d) = %E</span><span style="font-size:12pt">\</span><span style="font-size:12pt">n</span><span style="font-size:12pt">",</span><span style="font-size:12pt">er</span><span style="font-size:12pt">,</span><span style="font-size:12pt">ephi</span><span style="font-size:12pt">,</span><span style="font-size:12pt">ez</span><span style="font-size:12pt">,(double)</span><span style="font-size:12pt">arrCoord[ez][ephi][er][icErmphip[1]</span><span style="font-size:12pt">);</span></p><p>}</p><p><br></p><p>When I execute this example, I get this output:</p><p><span style="font-size:12pt">Phi_p(1,0</span><span style="font-size:12pt">,0) =</span><span style="font-size:12pt"> </span><span style="font-size:12pt">2</span><span style="font-size:12pt">.094395</span><span style="font-size:12pt">E+00</span></p><p><span style="font-size:12pt"></span>Phi_p(1,1,0) = 4.188790E+00<br></p><p>Phi_p(1,2,0) = 0.000000E+00</p><p><br></p><p>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">π / 3 and 4</span><span style="color:rgb(33,33,33);font-family:Calibri,Helvetica,sans-serif,serif,EmojiFont;font-size:12pt">π / 3 respectively. Thus, nothing is wrong here. </span></p><p><span style="color:rgb(33,33,33);font-family:Calibri,Helvetica,sans-serif,serif,EmojiFont;font-size:12pt">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">π instead of 0.</span></p><p><span style="color:rgb(33,33,33);font-family:Calibri,Helvetica,sans-serif,serif,EmojiFont;font-size:12pt"><br></span></p><p><span style="font-size:12pt">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><span style="font-size:12pt"></span><br></div></div></div></blockquote><div>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>one, so we choose 0.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div id="gmail-m_3112071783902229590gmail-m_-5207082944899504102divtagdefaultwrapper" dir="ltr" style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif"><p><span style="font-size:12pt">Thank you.</span></p><p><span style="font-size:12pt">Best regards,</span></p><p><span style="font-size:12pt"><br></span></p><p><span style="font-size:12pt">Zakariae Jorti</span></p></div><hr style="display:inline-block;width:561.25px"><div id="gmail-m_3112071783902229590gmail-m_-5207082944899504102divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" style="font-size:11pt"><b>From:</b><span> </span>Patrick Sanan <<a href="mailto:patrick.sanan@gmail.com" target="_blank">patrick.sanan@gmail.com</a>><br><b>Sent:</b><span> </span>Tuesday, March 23, 2021 11:37:04 AM<br><b>To:</b><span> </span>Jorti, Zakariae<br><b>Cc:</b><span> </span><a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a><br><b>Subject:</b><span> </span>[EXTERNAL] Re: Question about periodic conditions</font><div> </div></div><div>Hi Zakariae - sorry about the delay - responses inline below.<div><br></div><div>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><div><br></div><div><div><br><blockquote type="cite"><div>Am 23.03.2021 um 00:54 schrieb Jorti, Zakariae <<a href="mailto:zjorti@lanl.gov" target="_blank">zjorti@lanl.gov</a>>:</div><br><div><div id="gmail-m_3112071783902229590gmail-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"><div style="margin-top:0px;margin-bottom:0px">Hi, </div><div style="margin-top:0px;margin-bottom:0px"><br></div><div style="margin-top:0px;margin-bottom:0px">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">0 < r_min <= r <= r_max;                with <span style="font-size:12pt">r_max</span><span style="font-size:12pt"> > r_min</span></div><div style="margin-top:0px;margin-bottom:0px">phi_min = 0 <span style="font-size:12pt"><= r <= phi_max</span><span style="font-size:12pt"> <span> </span>= 2 </span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">π</span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">z_min <= z =< z_max;                    </span><span style="font-size:12pt">with </span><span style="font-size:12pt">z_max</span><span style="font-size:12pt"> > z_min.</span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-size:12pt"><br></span></div><div style="margin-top:0px;margin-bottom:0px"><font size="3">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"> on face centers. (</font><span style="font-size:12pt">dof0 = 0, dof1 = 1,dof2 = 1, dof3 = 0;</span><font size="3">). </font><font size="3"> </font></div></div></div></blockquote><blockquote type="cite"><div><div id="gmail-m_3112071783902229590gmail-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"><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><br></span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-size:12pt"></span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">I have two versions of my code: </span></div><div style="margin-top:0px;margin-bottom:0px"><font face="Calibri, Helvetica, sans-serif" size="3">1 - A first version<span> </span></font><font face="Calibri, Helvetica, sans-serif">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"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">2- A second version in which I set </span><span style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif">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">, </span><font face="Calibri, Helvetica, sans-serif" size="3">and DM_BOUNDARY_PERIODIC in the phi di</font><font face="Calibri, Helvetica, sans-serif">rection.</font><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"> </span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><br></span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">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"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">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><br></div><div>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><br></div><div>   <span> </span>x ---- x ---- x ----</div><div><br></div><div>as opposed to the non-periodic case,</div><div><br></div><div>   x ---- x ---- x ---- x</div><div><br></div><div><br></div><div><blockquote type="cite"><div><div id="gmail-m_3112071783902229590gmail-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"><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><br></span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">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"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">B_phi [phi = 0] = 1.0;</span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif">B_phi </span><span style="font-family:Calibri,Helvetica,sans-serif">[phi = 2</span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">π</span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">] = 1.0;</span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif">E_z [r, phi = 0</span><span style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif">] = 1/r;</span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><br></span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">E_z [r, phi = 2</span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">π</span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">] = 1/r;</span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><br></span></div><div style="margin-top:0px;margin-bottom:0px"><font face="Calibri, Helvetica, sans-serif">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">π with the periodic boundary conditions, </span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">is it sufficient for example </span><font face="Calibri, Helvetica, sans-serif" size="3">to </font><font face="Calibri, Helvetica, sans-serif">have</font><font face="Calibri, Helvetica, sans-serif" size="3"> only the following boundary conditions</font><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">:</span></div><div style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt">B_phi [phi = 0] = 1.0;</span></div><p style="margin-top:0px;margin-bottom:0px"><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"></span></p><div style="margin-top:0px;margin-bottom:0px"><span style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif">E_z [r, phi = 0</span><span style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif">] = 1/r ? </span><span style="font-family:Calibri,Helvetica,sans-serif;font-size:12pt"><br></span></div></div></div></blockquote><div><br></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><br></div><div>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><br><blockquote type="cite"><div><div id="gmail-m_3112071783902229590gmail-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"><div><span style="font-size:12pt;font-family:Calibri,Helvetica,sans-serif">Thank you.</span></div><div><font face="Calibri, Helvetica, sans-serif">Best regards,</font></div><div><font face="Calibri, Helvetica, sans-serif"><br></font></div><div><font face="Calibri, Helvetica, sans-serif">Zakariae Jorti</font></div></div></div></blockquote></div><br></div></div></div></div></blockquote></div><br clear="all"><div><br></div>--<span> </span><br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a></div></div></div></div></div></div></div></div></div></blockquote></div><br></div></div></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>