On Wed, Mar 31, 2021 at 6:40 AM Patrick Sanan <patrick.sa...@gmail.com> wrote:
> Am 31.03.2021 um 12:11 schrieb Matthew Knepley <knep...@gmail.com>: > > On Wed, Mar 31, 2021 at 3:21 AM Patrick Sanan <patrick.sa...@gmail.com> > 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/DMStagPopulateLocalToGlobalInjective.html > 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 <knep...@gmail.com>: >> >> On Wed, Mar 24, 2021 at 7:17 PM Jorti, Zakariae via petsc-users < >> petsc-us...@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.sa...@gmail.com> >>> *Sent:* Tuesday, March 23, 2021 11:37:04 AM >>> *To:* Jorti, Zakariae >>> *Cc:* petsc-us...@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 <zjo...@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/> >> >> >> > > -- > 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/>