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

[Jorti, Zakariae]

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?

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210324/67804765/attachment-0001.html>