i am trying to solve a poisson equation:

d^2/dx^2(PHI) + d^2/dy^2(PHI) = S

where S is a modified source term that satisfy the consequence of Green Thereom. the boundary conditions are the Neumann type, where at the boundary:

the normal derivative of PHI=0, i.e. d(PHI)/dn=0

By utilizing ADI, solution is obtained. but the problem is when i substitute the solution to check if the boundary conditions are satisfied, it turns out that :

d(PHI)/dn <> 0 (not equal to)

the question are: what can be the possible cause? how can this problem be tackled?

any advice is appreciated. regards, yfyap

When you have only Neumann type boundary conditions, the solution is not unique, i.e. any PHI+C will satisfy the Poisson equation.

As to how to resolve this -- boundary conditions of the mixed types (Dirichlet+Neumann+Robin) should be used. At least, analytical solution will be unique.

For discussion, think of this as the heat conduction equation with PHI the temperature. If you have a steady state (and you do, D PHI / Dt = 0), while adding energy to the system (S > 0) or removing it (S < 0), then S has to be balanced by energy crossing the boundary. If the normal derivitive of temperature on the boundary is everywhere 0, no energy crosses the boundary, and a steady state is not possible.

I'm sure this can be stated completely mathematically. Your problem is that your boundary condition is not compatible with your equation.

There are several possibilities.

If the normal derivative is zero everywhere on the boundary, then the Greens theorem would imply that the volume (or area in 2D) integral of the source term also be zero. Convergence to machine zero is possible only if discrete version of the Green's theorem is satisfied.

When you say "solution is obtained", what does that mean ? ADI is not a direct solver, it is iterative. What criterion do you use to stop the iterations and conclude that the solution is accurate enough. It is quite likely that the residual could still be quite high.

When you build tridiagonal matrices (for sweeping along each direction), do you build into it, the (zero gradient) boundary conditions.

ADI sometimes has convergence problems if you have a singular matrix and based on your BCs, it is evident that you have one.

thanks for the precious advice. iteration stops when

{ PHI[i][j][t+1]-PHI[i][j][t] } < {permitted error}

when i said "solution is obtained", what i mean is that if i reduce the {permitted error}, let's say by a factor of 10e-3, it gives identical solution.

in the previous posting i forgot to mention that since any PHI+C is solution, then i fixed PHI[0][0]=constant, then, the solution should be unique, right? the BCs are actually of the Robin type.

i have checked that:

int{ d^2/dx^2(PHI) + d^2/dy^2(PHI) }dA - int{S}dA < 10-e16

note: int = integration of

can this implies that "discrete version of the Green's theorem is satisfied."?

i am checking if the residual is large and the imposition of the BCs.

thanks for the precious advice.

What kind of discetization method do you use (FDM, FVM, FVEM or FEM). Myself or others people could give you better answers.

i used finite difference method.

thanks. both positive and negative values of S exist for my problem, then, there are sink and source within the boundary. i have checked that: int{ d^2/dx^2(PHI) + d^2/dy^2(PHI) }dA - int{S}dA < 10-e16 note: int = integration of

in the previous posting i forgot to mention that since any PHI+C is solution, then i fixed PHI[0][0]=constant, then, the solution should be unique, right? the BCs are actually of the Robin type.

check your code everywhere carefully!

thanks. both positive and negative values of S exist for my problem, then, there are sinc and source within the boundary. i have checked that:

int{ d^2/dx^2(PHI) + d^2/dy^2(PHI) }dA - int{S}dA < 10-e16

note: int = integration of

in the previous posting i forgot to mention that since any PHI+C is solution, then i fixed PHI[0][0]=constant, then, the solution should be unique, right? the BCs are actually of the Robin type.

thanks. doing the checking...

1) Do you have any unknowns on the domain boundary?

2) How do you impose the boundary condition?

Quite likeky you have to check out how do you impose and implement the boundary conditions with your FD approach.

I wonder what scheme CDS or upwind you use to implement the boundary conditions. If you use CDS, you will impose the zero gradient condition across the wall but this does not garantee when you do the check out that the gradient is zero on the wall. This doesn't happen if you use upwind.

Any case, you should take a look at how you implement the boundary conditions on your domain.

thanks for the advice.

thanks for your concerns. i applied some sort of linearized scheme so that there are some terms in the source term evaluated in the preceding station. other than that, PHI is the only unknown at the boundary. to impose the boundary conditions, i used the reflection method. it is found that i get the same equations at the boundary if i used three-point-one-sided finite differece formula. it should be ok, right? regards, yfyap