MKL Cookbook Example 1 -- a Critique

[mecej4]

The MKL Cookbook at https://software.intel.com/en-us/mkl_cookbook has a set of interesting and useful examples of solving physics problems using MKL. After reading the first of the cookbook examples, which is nicely described at https://software.intel.com/en-us/node/507039, and trying to verify the results given by the cookbook code at https://software.intel.com/en-us/mkl_cookbook_samples, I have some critical comments. I hope someone will view these comments as constructive criticism.

1. The brief description given at https://software.intel.com/en-us/mkl_cookbook says "Finding an approximate solution to a nonlinear equation", which would suggest to most people that the example concerns solving a nonlinear equation in one variable. The example, in fact, is about solving a nonlinear PDE using finite-difference discretization, and solving the resulting system of nonlinear algebraic equations using Pardiso.

2. The nonlinearity is introduced through the temperature dependence of thermal conductivity, as described on https://software.intel.com/en-us/node/507039 under "Discussion" (two different but similar symbols are used for temperature -- v and the Greek letter 'nu' -- please fix this). To treat this nonlinearity, the discretized equations are solved iteratively, using the solution of each iteration to recompute the coefficients to be used in the next iteration. In principle, this is fine, and there are many applications where such iterations are necessary. However, in the specific problem chosen this iteration is unnecessary. Gustav Kirchhoff showed in 1894 that using w = \int \mu dv as the dependent variable reduces the variable conductivity problem to an equivalent constant conductivity problem.

3. I added 'printf' statements to the code in file pardiso_nonlinear.c to print out the final solution for temperature. The solution is not only incorrect in the interior of the unit cube domain, but shows values of 1 at the faces of the cube instead of the specified value of 0. I have not debugged the 300-line code to find out what is wrong. Instead, I give below two different ways of solving the problem and the results, so that you (Intel/MKL group) can debug the C code in the cookbook.

4. The example code employs a mesh with 6³ nodes, of which 4³ are interior nodes with unknown values of the dependent variable, and solves for these unknowns by calling Pardiso with n = 216. Exploiting the symmetry of the problem enables us to recognize that there are only four distinct temperature values! Let us name these four values as follows:

      x     y     z     w

   0.2  0.2  0.2     p

   0.2  0.4  0.2     q

   0.4  0.4  0.2     r

   0.4  0.4  0.4     s

Then, the heat equation is represented by A w = b, where A =

     6    -3     0     0
    -1     5    -2     0
     0    -2     4    -1
     0     0    -3     3

and b = h², with h = 0.2 being the mesh size, and w = [p q r s]^T. The solution is w = [   0.0196 ,   0.0260,    0.0351,     0.0484]^T. By reversing the Kirchhoff transformation, we recover the solution for temperature as T =   [  0.0180,    0.0233,    0.0305,    0.0403]^T.

5. Here is simple Fortran code to solve the problem using Gauss-Seidel iteration to solve the linear equations for w (denoted T in the code). To keep the code short, I have not exploited symmetry, using which would have reduced the number of variables. The output values agree with those given above in Item 4.

! Poisson equation in unit cube, with T = zero at faces. program lapl3d implicit none integer, parameter :: nx=5, ny=5, nz=5 ! number of segments real, dimension (0:nx,0:ny,0:nz) :: T = 0.0 ! Kirchhoff flow temperature real :: h=1.0/nx, dT, Tn integer :: i,j,k,iter=0 do dT=0 do i=1,nx-1 do j=1,ny-1 do k=1,nz-1 Tn = (T(i-1,j,k)+T(i+1,j,k)+ & T(i,j-1,k)+T(i,j+1,k)+ & T(i,j,k-1)+T(i,j,k+1) + h*h)/6.0 ! discrete Poisson eq. dT = max(dT,abs(T(i,j,k)-Tn)) T(i,j,k) = Tn ! Gauss-Seidel end do end do end do iter=iter+1 if(mod(iter,10).eq.0)write(*,'(1x,I4,2x,E10.3)')iter,dT if(iter.gt.500.or.dT.lt.1e-8)exit end do ! recover temperatures from Kirchhoff flow temperatures do i=1,nx-1 do j=1,ny-1 write(*,10)i,j,(sqrt(0.2*T(i,j,k)+0.01)-0.1,k=1,nz-1) end do end do 10 format(1x,2I3,<nx-1>F10.5) end program

[mecej4]

Historical point: Kirchhoff died in 1887 (see http://commons.wikimedia.org/wiki/File:Grabstein.Gustav.Kirchhoff.jpg) . The date reference in Item 2 is to his posthumously published lecture notes.

[Ying_H_Intel]

Hi Mecej4, 

Thank you much for the review of the cookbook and constructive comments.  I will escalate them to the corresponding authors to further check.

Thanks

Ying H.

Intel MKL Support  

[mecej4]

Intel, please respond!

[Ying_H_Intel]

Hi Mecej4,

Sorry for the delay. we are still investigating the problem internally. Some one in our team will cover it in next few days.

Thanks

Ying

[mecej4]

Ying,

Thanks, I'll wait. Since nearly a month had passed since I made the posting, I wrote a line as a reminder.

[MariaZh]

Hello mecej4,

Firstly I wand to thank you for all of your comments and apologize for the delay in answering.

I took into account your remarks and suppose that the result requires further consideration. I would like to take time to carefully investigate this problem.

Best regards,
Maria.

[Jeanette_F_Intel]

Hi, mecej4--

We posted an updated version of the recipe here: https://software.intel.com/en-us/node/507039, and updated sample code here: https://software.intel.com/en-us/mkl_cookbook_samples. I hope this fixes the issues that you mentioned.

Thank you so much for your valuable and constructive feedback.

--Jeanette