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

[Paul Stephen Urbanczyk]

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?

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
>
>