[petsc-dev] [petsc-users] [EXTERNAL] Re: Question about periodic conditions
Matthew Knepley
knepley at gmail.com
Wed Mar 31 05:11:30 CDT 2021
On Wed, Mar 31, 2021 at 3:21 AM Patrick Sanan <patrick.sanan at 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.
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.
THanks,
Matt
> Am 25.03.2021 um 00:20 schrieb Matthew Knepley <knepley at gmail.com>:
>
> 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/>
>
>
>
--
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-dev/attachments/20210331/0cc8d781/attachment.html>
More information about the petsc-dev
mailing list