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

[David Scott]

Thanks for the replies. I obviously need to think about the problem 
more. I was aware that the solution was not fully determined and I had 
created a MatNullSpace to deal with the arbitrary constant.

I have contacted the original author of the code (I am just modifying it 
to use PETSc) to find out more about the formulation of the problem. I 
should point out that the the actual problem I want to deal with does 
not have a zero RHS.

This systems of equations is only part of the complete system and only 
the gradient of the solution is required elsewhere so the arbitrary 
constant has no physical significance.

David

On 10/01/2013 19:18, Sanjay Govindjee wrote:
> Perhaps more simply, consider the 1D version of your problem:
>
> u''=0 on [0,1]   with "Boundary Conditions"  u'(0)=a and u'(1) =b
>
> The solution is u(x) = a x + c  (for any! value of c).
> This is not going to solve happily :(
>
>
> On 1/10/13 8:27 AM, Jed Brown wrote:
>> ... and you describe a "boundary condition" as "second derivative is
>> zero", which is not a boundary condition, making your problem
>> ill-posed. (Indeed, consider the family of problems in which you
>> extend the domain in that patch and apply _any_ boundary conditions in
>> the extended domain. All of those solutions are also solutions of your
>> problem with "second derivative is zero" on your "boundary".)
>>
>>
>> On Thu, Jan 10, 2013 at 9:44 AM, David Scott <d.scott at ed.ac.uk
>> <mailto:d.scott at ed.ac.uk>> wrote:
>>
>>     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>
>>         <mailto: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.
>>
>>
>


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