[petsc-dev] [petsc-users] [EXTERNAL] Re: Question about periodic conditions
Matthew Knepley
knepley at gmail.com
Wed Mar 31 06:08:03 CDT 2021
On Wed, Mar 31, 2021 at 6:40 AM Patrick Sanan <patrick.sanan at gmail.com>
wrote:
> Am 31.03.2021 um 12:11 schrieb Matthew Knepley <knepley at gmail.com>:
>
> 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.
>
> 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 <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/>
>
>
>
--
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/319595d9/attachment-0001.html>
More information about the petsc-dev
mailing list