On Wed, Mar 31, 2021 at 6:40 AM Patrick Sanan <[email protected]> wrote:
Am 31.03.2021 um 12:11 schrieb Matthew Knepley <[email protected]>:
On Wed, Mar 31, 2021 at 3:21 AM Patrick Sanan <[email protected]> wrote:
(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).
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 is overkill here and would not give you any added functionality.
Ah, cool, I was wondering what the coordinate field function was for. Does the DMSetPeriodicity() stuff get used anymore, now that you have a different approach? Quickly looking it seems like it gets passed around as you create DMs from existing ones, and it's perhaps used in some Plex output functionality.
We use it to identify that the mesh is periodic and in what directions, and use the length if we have to figure out the coordinates.
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.
Computational work or implementation work?
implementation.
I already have logic in DMStag to support something related, which was required for L2G with INSERT_VALUES in the periodic, 1-rank case, where multiple local points can respond to the same global point. So building the local coordinates with the far point and then doing L2G doesn't sound bad to implement, at least.
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMSTAG/DMStagPo...
That may work. It is a problem in Plex because that strategy destroys the topological queries. However, it sounds like you do not have that problem. Thanks, Matt
THanks,
Matt
Am 25.03.2021 um 00:20 schrieb Matthew Knepley <[email protected]>:
On Wed, Mar 24, 2021 at 7:17 PM Jorti, Zakariae via petsc-users < [email protected]> 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 <[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]>:
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/>
-- 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/>
-- 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/>