[petsc-dev] [petsc-users] [EXTERNAL] Re: Question about periodic conditions

[Patrick Sanan]

(moving to petsc-dev) 

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.

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.

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:
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.
2.  to use PIC methods (DMSwarm), where we need a way to determine if a particle is in the last cell.

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:

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). 



> Am 25.03.2021 um 00:20 schrieb Matthew Knepley <knepley at gmail.com>:
> 
> On Wed, Mar 24, 2021 at 7:17 PM Jorti, Zakariae via petsc-users <petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>> wrote:
> Hi Patrick,
> 
> 
> 
> Thanks for your responses. 
> 
> 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.
> 
> 
> 
> 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. 
> 
> Let me give you an example in cylindrical coordinates with a 3x3x3 DMStag mesh:  
> 
>  
> 
> 
> 
> 
> PetscInt      startr,startphi,startz,nr,nphi,nz,d;
> 
> 
> PetscInt      er,ephi,ez,icErmphip[3];
> 
> 
> DM                dmCoorda, coordDA;
> 
> Vec               coordaLocal;
> 
> PetscScalar       ****arrCoord;
> 
> PetscScalar       surf;
> 
> 
> 
> 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);
> 
> 
> 
> DMSetFromOptions(coordDA);
> 
> DMSetUp(coordDA);
> 
> 
> 
> DMStagGetCorners(coordDA,&startr,&startphi,&startz,&nr,&nphi,&nz,NULL,NULL,NULL);
> 
> 
> 
> DMGetCoordinateDM(coordDA,&dmCoorda);
> 
> DMGetCoordinatesLocal(coordDA,&coordaLocal);
> 
>   DMStagVecGetArrayRead(dmCoorda,coordaLocal,&arrCoord);
> 
> 
> 
> for (d=0; d< 3; ++d){
> 
> DMStagGetLocationSlot(dmCoorda,UP_LEFT,d,&icErmphip[d]);
> 
> }
> 
> 
> 
> er = 1; ez = 0;
> 
> for (ephi=0; ephi< 3; ++ephi){
> 
>  PetscPrintf(PETSC_COMM_WORLD,"Phi_p(%d,%d,%d) = %E\n",er,ephi,ez,(double)arrCoord[ez][ephi][er][icErmphip[1]);
> 
> }
> 
> 
> 
> When I execute this example, I get this output:
> 
> Phi_p(1,0,0) = 2.094395E+00
> 
> Phi_p(1,1,0) = 4.188790E+00
> 
> Phi_p(1,2,0) = 0.000000E+00
> 
> 
> 
> Note here that the first two lines correspond to 2π / 3 and 4π / 3 respectively. Thus, nothing is wrong here. 
> 
> But the last line should rather give 2π instead of 0.
> 
> 
> 
> I understand that degrees of freedom should be the same on both sides of the boundary, but should the coordinates not be preserved? 
> 
> 
> 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
> one, so we choose 0.
> 
>   Thanks,
> 
>      Matt 
> Thank you.
> 
> Best regards,
> 
> 
> 
> Zakariae Jorti
> 
> From: Patrick Sanan <patrick.sanan at gmail.com <mailto:patrick.sanan at gmail.com>>
> Sent: Tuesday, March 23, 2021 11:37:04 AM
> To: Jorti, Zakariae
> Cc: petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>
> Subject: [EXTERNAL] Re: Question about periodic conditions
>  
> Hi Zakariae - sorry about the delay - responses inline below.
> 
> 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. 
> 
> 
>> Am 23.03.2021 um 00:54 schrieb Jorti, Zakariae <zjorti at lanl.gov <mailto:zjorti at lanl.gov>>:
>> 
>> Hi, 
>> 
>> I implemented a PETSc code to solve Maxwell's equations for the magnetic and electric fields (B and E) in a cylinder: 
>> 0 < r_min <= r <= r_max;                with r_max > r_min
>> phi_min = 0 <= r <= phi_max  = 2 π
>> z_min <= z =< z_max;                    with z_max > z_min.
>> 
>> I am using a PETSc staggered grid with the electric field E defined on edge centers and the magnetic field B defined on face centers. (dof0 = 0, dof1 = 1,dof2 = 1, dof3 = 0;).  
>> 
>> I have two versions of my code: 
>> 1 - A first version in which I set the boundary type to DM_BOUNDARY_NONE in the three directions r, phi and z
>> 2- A second version in which I set the boundary type to DM_BOUNDARY_NONE in the r and z directions, and DM_BOUNDARY_PERIODIC in the phi direction. 
>> 
>> 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. 
>> Is it normal?  
> 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 .
> 
> If you consider a 1-d example with 1 dof on vertices and cells, with three elements, the periodic case looks like this, globally,
> 
>     x ---- x ---- x ----
> 
> as opposed to the non-periodic case,
> 
>    x ---- x ---- x ---- x
> 
> 
>> 
>> 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:
>> B_phi [phi = 0] = 1.0;
>> B_phi [phi = 2π] = 1.0;
>> E_z [r, phi = 0] = 1/r;
>> E_z [r, phi = 2π] = 1/r;
>> 
>> Assuming that values at phi = 0 should be the same as at phi=2π with the periodic boundary conditions, is it sufficient for example to have only the following boundary conditions:
>> B_phi [phi = 0] = 1.0;
>> E_z [r, phi = 0] = 1/r ? 
> 
> Yes - this is the intention, since the boundary at phi = 2 * pi is represented by the same entries in the global vector. 
> 
> Of course, you need to make sure that your continuous problem is well-posed, which in general could change when using different boundary conditions. 
> 
>> Thank you.
>> Best regards,
>> 
>> Zakariae Jorti
> 
> 
> 
> -- 
> 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
> 
> https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20210331/2c0aa0a9/attachment-0001.html>