[petsc-users] ksp/dmmg laplace examples

[Darach Golden]

Hi,

I am planning on implementing a basic 3D Poisson solver using petsc --
preferably using the DMMG interface to multigrid with finite
difference.  Initially I'm using Dirichlet boundary conditions.  To
get started I've been looking at the following multigrid examples in
1D, 2D and 3D respectively

1) $PETSC_DIR/src/ksp/ksp/examples/tutorials/ex25.c
2) $PETSC_DIR/src/ksp/ksp/examples/tutorials/ex29.c
3) $PETSC_DIR/src/ksp/ksp/examples/tutorials/ex22.c

In each of these cases I've set any constants in the code so that the
problem is reduced to the basic laplacian/poisson

I have some questions about the finite difference formulation.  The
questions may very well be naive so apologies in advance for that. If
they are, could you point me off to some standard
reference/examples/slides?

In each of the examples above, the linear system to be solved seems to
contain entries for the total number of nodes (including boundary
nodes), with rows associated with boundary nodes containing values on
the diagonal only (in the Dirichlet case).  Is this correct?

It's these diagonal matrix values that are confusing me:

1) ex25
--------------------------------------------------------
 Partial differential equation

   d  (1 + e*sine(2*pi*k*x)) d u = 1, 0 < x < 1,
   --                        ---
   dx                        dx
with boundary conditions

   u = 0 for x = 0, x = 1
--------------------------------------------------------

In ComputeMatrix() the diagonal values in the matrix associated with
boundary nodes seem to be 2.0, and the rows associated with internal
nodes are divided by h. while in ComputeRHS(), the RHS appears to be
multiplied by h, Why is the value 2.0 in the diagonal?

2) ex29
--------------------------------------------------------
Inhomogeneous Laplacian in 2D. Modeled by the partial differential
equation

   div \rho grad u = f,  0 < x,y < 1,

with forcing function

   f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}

with Dirichlet boundary conditions

   u = f(x,y) for x = 0, x = 1, y = 0, y = 1
--------------------------------------------------------

The RHS seems to be multiplied by hx * hy, while the values on the
diagonal appear to be 2.0*rho*((hx/hy) + (hy/hx)) where I take rho=1.
Why is this?

Also, is the forcing function in ComputeRHS() 
      f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}?

3) ex22
--------------------------------------------------------
Laplacian in 3D. Modeled by the partial differential equation

   - Laplacian u = 1,0 < x,y,z < 1,

with boundary conditions

   u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
--------------------------------------------------------

The RHS seems to be multiplied by hx*hy*hz, while the values on the
diagonals (associated with boundary nodes) appear to be

2.0*((hx*hy/hz) + (hx*hz/hy) + (hy*hz/hx)).

Why?


Darach