[petsc-users] Dealing with 2nd Derivatives on Boundaries

[David Scott]

All right, I'll say it differently. I wish to solve
   div.grad phi = 0
with the boundary conditions that I have described.

On 10/01/2013 15:28, Jed Brown wrote:
> Second derivative is not a boundary condition for Poisson; that is the
> equation satisfied in the interior. Unless you are intentionally
> attempting to apply a certain kind of outflow boundary condition (i.e.,
> you're NOT solving Laplace) then there is a problem with your
> formulation. I suggest you revisit the continuum problem and establish
> that it is well-posed before concerning yourself with implementation
> details.
>
>
> On Thu, Jan 10, 2013 at 9:03 AM, David Scott <d.scott at ed.ac.uk
> <mailto:d.scott at ed.ac.uk>> wrote:
>
>     Hello,
>
>     I am solving Poisson's equation (actually Laplace's equation in this
>     simple test case) on a 3D structured grid. The boundary condition in
>     the first dimension is periodic. In the others there are Von Neumann
>     conditions except for one surface where the second derivative is
>     zero. I have specified DMDA_BOUNDARY_NONE in these two dimensions
>     and deal with the boundary conditions by constructing an appropriate
>     matrix. Here is an extract from the Fortran code:
>
>                    if (j==0) then
>                      ! Von Neumann boundary conditions on y=0 boundary.
>                      v(1) = 1
>                        col(MatStencil_i, 1) = i
>                        col(MatStencil_j, 1) = j
>                        col(MatStencil_k, 1) = k
>                      v(2) = -1
>                        col(MatStencil_i, 2) = i
>                        col(MatStencil_j, 2) = j+1
>                        col(MatStencil_k, 2) = k
>                      call MatSetValuesStencil(B, 1, row, 2, col, v,
>     INSERT_VALUES, ierr)
>                    else if (j==maxl) then
>                      ! Boundary condition on y=maxl boundary.
>                      v(1) = 1
>                        col(MatStencil_i, 1) = i
>                        col(MatStencil_j, 1) = j
>                        col(MatStencil_k, 1) = k
>                      v(2) = -2
>                        col(MatStencil_i, 2) = i
>                        col(MatStencil_j, 2) = j-1
>                        col(MatStencil_k, 2) = k
>                      v(3) = 1
>                        col(MatStencil_i, 3) = i
>                        col(MatStencil_j, 3) = j-2
>                        col(MatStencil_k, 3) = k
>                      call MatSetValuesStencil(B, 1, row, 3, col, v,
>     INSERT_VALUES, ierr)
>                    else if (k==0) then
>
>
>     Here the second clause deals with the second derivative on the boundary.
>
>     In order for this code to work I have to set the stencil width to 2
>     even though 'j-2' refers to an interior, non-halo
>     point in the grid. This leads to larger halo swaps than would be
>     required if a stencil width of 1 could be used.
>
>     Is there a better way to encode the problem?
>
>     David
>
>     --
>     The University of Edinburgh is a charitable body, registered in
>     Scotland, with registration number SC005336.
>
>


-- 
Dr. D. M. Scott
Applications Consultant
Edinburgh Parallel Computing Centre
Tel. 0131 650 5921

The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.