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 <[email protected]>
Sent: Tuesday, March 23, 2021 11:37:04 AM
To: Jorti, Zakariae
Cc: [email protected]
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
<[email protected]<mailto:[email protected]>>:
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
