Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6534 Discussions

Solving diffusion type equation using Poisson Solver



I am trying to solve a 3D diffusion type equation with periodic in X and Y and Neumann boundary condition  in Z direction using the MKL Poisson Solver and facing couple of problems.

First of all the BCTYPE if i use 'PPPPNN' and put  

bd_az[i + j * (nx+1)]= 0.0
bd_bz[i + j * (nx+1)]= 0.0

in the Z boundary, is it considering the Zero Neumann condition acurately at the boundaries?

and the other question is if i write the diffusion equation in terms of poisson equation then my RHS or the 'f' will be time dependent as du/dt. So, in that case will there be any conflict between the time scheme (u_old , u_new)?

Or is there any other way to use the MKL Poisson solver for the time dependent equations with the above mentioned boundary conditions ?

Please explain.




0 Kudos
5 Replies

Hi Swagnik, 


1. Yes if you assign zeros in the boundary arrays under Neumann Boundary condition, then it will be the case of zero Neumann BC which means that the edge of the boundary is isolated and no flux of heat enters or exits the boundary. 


2. Will you describe more of the case or provide your test case for better understanding? 


Thank you

Black Belt

Here is one point which you may have overlooked, and can cause numerical difficulties. If the "diffusion type equation" is homogeneous, and so are the boundary conditions (periodic; Neumann), if φ(x,y,z,t) is a solution, so is φ(x,y,z,t) + C, where C is a constant. To overcome this problem, you may have to pin one node to a specific value of φ.



I would like to address your question about whether we can use a time derivative in the Poisson equation solver.  Let us analyze the repercussions:  In other words you want to solve the following heat equation problem:

Find u(x,y,z,t) on 0 <= t <= T and (x,y,z) in some closed rectangular region Ω = [x0,x1]x[y0,y1]x[z0,z1] such that

du/dt = -Δ​u

with periodic boundary conditions in x and y on and neumann boundary conditions in z with an initial state

u(x,y,z, t=0) = u0(x,y,z)  for all (x,y,z) in Ω​.


This is fundamentally a different class of equations (parabolic equations) than the poisson equation (elliptic equation) and as such, requires different tools and techniques for solving.  Typical spatial discretizations of the elliptic equations reduce the differential equation to a linear algebra system of the type

A U = B     

where A is a sparse (hopefully invertible, see mecej4's comment about difficult issues with certain BC's) matrix and U,B are vectors.  B comes from a discretization of the right hand side and U is the vector of unknown values which then can be converted back to an approximate solution on the region Ω.

Now, if we apply the same spatial discretization techniques to the heat equation above, we end up with a (intermediate) vector ODE system with continuous dependence on time but that is discretized in space:

Find U(t) for 0<= t <= T such that

M dU(t)/dt = A(t) U(t)

with initial conditions

M U(0) = U0.

The mass matrix M is used here and represents the discretization of space of the solution  (when using finite element or finite volume techniques, it may be something different than the identity, but when using finite difference approaches it could be thought of as the identity matrix I, or when applying the Fourier transform or other spectral techniques to convert spatial derivatives to algebraic operations, it is also typically the identity, I).  Note that we have explicitly denoted the remaining dependence on time.  You can now use your favorite ODE solver technique to solve this.  For most people the first thing to try is forward or backward Euler temporal discretization. I will proceed with the unconditionally stable Backward Euler technique to discretize in time with time step, dt>0, and using superscript to denote vector solution at time t^{k},  U^k := U(t^k) := U(0 + k*dt)   (Recalling that U is itself a vector representing spatial unknowns): 

M [ (U^{k+1} - U^{k}) / dt ] = A(t^{k+1}) U^{k+1}

M U(0) = U0.

This can further be simplified (moving the unknowns U^{k+1} to left and knowns U^{k} to the right) to the update equation

(M - dt A(t^{k+1}) )U^{k+1} = M U^{k}

with initial conditions

M U^0 = U0.

So if A depends on time (it doesn't in the simplest heat equation case above) then A must be recomputed at each time step.  Otherwise, you must somehow assemble the sparse matrix  M - dt A   and then you have a series of linear algebra equations to solve as you propagate through time.  

The Poisson solver does not do this for you as it is built only for spatial discretizations and elliptic equations.  It uses a Fourier Transform approach to spatial discretization to convert the derivatives into algebra which then is solved using an LU factorization solver.  It only assembles the A U = B system. 

[As a note, the above material on spatial and temporal discretizations of ODEs and PDEs is a very brief overview of standard techniques covered in a first semester graduate course in the numerical analysis of partial differential equations or as a second or third semester of undergraduate numerical analysis. ]


Thanks, a lot for the answers. But now i am really in a trouble . I just want to solve Poisson's equation in 3d with the same above condition i.e.

periodic in x and Y and non-periodic and neumann in Z. My equation is Δ2 u = Δ·p . where p is a vector (px(i,j,k), py(i,j,k), pz(i,j,k)).

Now whwn i am applying BCTYPE as 'PPPPDD' and applying boundary condition i am getting the exact  ans but when i am applying the neumann condition it shows and error as:

The problem is degenerate due to rounding errors. The approximate solution
that provides the minimal Euclidean norm of the solution will be computed.

Double precision 3D Poisson example FAILED to compute the solution..."

I want to apply zero neumann boundary condition on Z. Really not understanding. Please help. I have attached the source file. Please help.






I have written up a couple thoughts on your problem and attached them here below.  It turns out that as written, your problem is not consistent, meaning your right hand side data and your boundary conditions cannot both be true at the same time.  I propose a modification and a family of weak solutions you can test against. 

If I may ask, what is the purpose of solving this specific problem?  Is there a scientific problem which reduces to this system or did you just make it up as a test ?

I don't know your background so I wrote out my thoughts assuming you have some exposure to applied analysis or functional analysis but I don't go deep into it at all.  I use some standard results which can be looked up if interest exists. 

I have also not modified your code problem but merely suggested what you could do to try fixing it.

This is outside the realm of what the Poisson solver was originally designed to do (it is a toy problem solver and has a few basic supported boundary conditions) but it should probably be able to solve a problem like this that actually has a solution.

I look forward to your thoughts.