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

[Matthew Knepley]

On Wed, Mar 24, 2021 at 7:17 PM Jorti, Zakariae via petsc-users <
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>
> *Sent:* Tuesday, March 23, 2021 11:37:04 AM
> *To:* Jorti, Zakariae
> *Cc:* 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>:
>
> 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-users/attachments/20210324/edf6a82c/attachment.html>