[petsc-users] Structured Finite Difference Method With Periodic BC Help

[Barry Smith]

> On May 6, 2016, at 5:09 PM, Paul Stephen Urbanczyk <gomer at stanford.edu> wrote:
> 
> Hello Barry,
> 
> Thank you for responding. I think I'm following you. I set up a little 1-D test problem, and used the built-in DMDASetUniformCoordinates function. Then, by turning on and off the periodic BCs, I see the change in behavior. I think I was making a mistake with the right boundary indexing.
> 
> However, this raises another question: shouldn't the coordinates behave differently than other field variables? In other words, shouldn't the periodic ghost points at the left end have coordinates (0-h) and (0-2h), rather than (1-h) and (1-2h)? And, at the right side, shouldn't the x-coordinates be 1 and 1+h?

   I don't know. But I believe they should not; periodic boundary conditions can be interpreted as having the domain [0,1) connected at the two ends so there would never be coordinates less than 0 or 1 or larger. Another way, I think, of saying the same thing is that if you take a point x in domain and add some positive distance d then the coordinate of the result is give by fractionalpart(x + d) something similar can be done for subtracting some distance. I remember vaguely there are precise mathematical ways to define this kind of thing but I don't think it matters. 

  Barry

> 
> Apologies for what may be somewhat remedial questions, but I'm somewhat new to PETSc and have not actually programmed a periodic case before.
> 
> Thanks for all of your help!
> 
> -Paul
> 
> ________________________________________
> From: Barry Smith <bsmith at mcs.anl.gov>
> Sent: Tuesday, May 3, 2016 7:26:38 PM
> To: Paul Stephen Urbanczyk
> Cc: petsc-users at mcs.anl.gov
> Subject: Re: [petsc-users] Structured Finite Difference Method With Periodic BC Help
> 
>> On May 3, 2016, at 7:38 PM, Paul Urbanczyk <gomer at stanford.edu> wrote:
>> 
>> Hello,
>> 
>> I'm trying to implement a CFD code using finite differences on a structured mesh with PETSc.
>> 
>> I'm having a bit of trouble understanding how to properly use the periodic boundary condition with distributed arrays, and need some clarification.
>> 
>> If I set the boundaries in the x-direction to be DM_BOUNDARY_PERIODIC, and set a stencil width of two, there should be two ghost cells in the x-direction at each end of the mesh, which looks like it's happening just fine.
>> 
>> However, it seems that the assumption being made by PETSc when filling in these values is that the mesh is a cell-centered finite difference, rather than a vertex-centered finite difference. Is there a way to change this behavior? In other words, I need the first ghost cell on each side to correspond to the opposite side's first interior point, rather than the opposite boundary point.
> 
>   DMDA doesn't really have a concept of if the values are cell-centered or vertex centered. I think this issue you are facing is only an issue of labeling, not of real substance; in some sense it is up to you to interpret the meaning of the values. In your case with vertex centered values here is how I see it in a picture. Consider a domain from [0,1) periodic with n grid points in the global vector
> 
>    x =  0                                                        1-2h 1-h     1
>    i  =  0   1   2                                                       n-1
> 
> So you divide the domain into n sections and label the vertices (left end of the sections) starting with zero, note that the last section has no right hand side index because the value at x=1 is the same as the value at x=0  now if you have two processors and a stencil width of two the local domains look like
> 
>    x = 1-2h 1-h   0
>    i =   n-2   n-1  0 1  2  ....
> 
> 
>    x =                                                       1-2h 1-h 1=0*  1+h
>    i  =                                                        n-2   n-1  0     1
> 
> This is what the DMDA will deliver in the local vector after a DMGlobalToLocalBegin/End
> 
> In periodic coordinates the location x=1 is the same as the location x=0
>> 
>> If there is not a way to change this behavior, then I need to set my own ghost cells; however, I'm having trouble implementing this...
>> 
>> If I change the boundaries to DM_BOUNDARY_GHOSTED, with a stencil width of two, I have the necessary ghost cells at either end of the mesh. I then try to do the following:
>> 
>> 1) Scatter the global vector to local vectors on each rank
>> 
>> 2) Get a local array referencing the local vectors
>> 
>> 3) Calculate ghost values and fill them in the appropriate local arrays
>> 
>> 4) Restore the local vectors from the arrays
>> 
>> 5) Scatter the local vector info back to the global vector
>> 
>> However, if I then re-scatter and look at the local vectors, the ghost cell values are zero. It seems as though the ghost values are lost when scattered back to the global vector.
> 
>   This is nuts; in theory the basic DMDA does what one needs to handle periodic boundary conditions. No need to try to implement it yourself.
> 
>   Of course we could always have bugs so if you have a problem with my reasoning send a simple 1d code that demonstrates the issue.
> 
>   Barry
> 
>> 
>> What am I doing wrong here?
>> 
>> Thanks for your help!
>> 
>> -Paul
>> 
>> 
> 
>