# Poisson Library Code Examples

## Cartesian case

The code below computes an approximate solution of a 2D Poisson problem
in the rectangle 0<x<1, 0 <y<1.

The following boundary conditions are imposed for the problem:

- Dirichlet boundary condition
  for 0≤y≤1.
- Neumann boundary condition
  for 0<x<1.

The exact solution is known to be
Example “C Example to Solve 2D Poisson Problem” implements the computations in C and Example “Fortran Example to Solve 2D Poisson Problem” provides Fortran 90 code. Mind that PL interface cannot be invoked from Fortran 77 due to restrictions imposed by the use of Intel MKL FFT interface.

The algorithm of computing the approximate solution of a Poisson Problem uses Poisson Library routines. Details of the Poisson problem being solved and the errors are printed out between the computed solution and the exact one.

Upon successful execution of Example “C Example to Solve 2D Poisson Problem” the following text is printed out (Example “Fortran Example to Solve 2D Poisson Problem” produces similar output):

Example of use of MKL Poisson Library

**********************************************

This example gives the solution of 2D Poisson problem

with the equation -u_xx-u_yy=f(x,y), 0<x<1, 0<y<1,

f(x,y)=(8*pi*pi)*sin(2*pi*x)*sin(2*pi*y),

and with the following boundary conditions:

u(0,y)=u(1,y)=1 (Dirichlet boundary conditions),

-u_y(x,0)=-2.0*pi*sin(2*pi*x) (Neumann boundary condition),

u_y(x,1)= 2.0*pi*sin(2*pi*x) (Neumann boundary condition).

-----------------------------------------------------------------------

In general, the error should be of order O(1.0/nx^2+1.0/ny^2)

For this example, the value of nx=ny is 6

The approximation error should be of order 1.0e-01, if everything is OK

-----------------------------------------------------------------------

DOUBLE PRECISION COMPUTATIONS

=======================================================================

The number of mesh intervals in x-direction is nx=6

The number of mesh intervals in y-direction is ny=6

In the mesh point (0.167,0.000) the error between the computed and the true solution

is equal to -7.505e-02

In the mesh point (0.167,0.167) the error betweenthe computed and the true solution

is equal to 4.432e-02

In the mesh point (0.167,0.333) the error betweenthe computed and the true solution

is equal to 6.309e-02

In the mesh point (0.167,0.500) the error between the computed and the true solution

is equal to -5.551e-16

In the mesh point (0.167,0.667) the error between the computed and the true solution

is equal to -6.309e-02

In the mesh point (0.167,0.833) the error between the computed and the true solution

is equal to -4.432e-02

In the mesh point (0.167,1.000) the error between the computed and the true solution

is equal to 7.505e-02

Double precision 2D Poisson example has successfully PASSED

through all steps of computation!

Note that actual figures in the error may slightly differ from the figures printed above depending on the architecture and operating system used to run the example.

C code of the computations is presented below:

## C Example to Solve 2D Poisson Problem

/*******************************************************************************

/* INTEL CONFIDENTIAL

/* Copyright(C) 2006-2008 Intel Corporation. All Rights Reserved.

/* The source code contained or described herein and all documents related to

/* the source code ("Material") are owned by Intel Corporation or its suppliers

/* or licensors. Title to the Material remains with Intel Corporation or its

/* suppliers and licensors. The Material contains trade secrets and proprietary

/* and confidential information of Intel or its suppliers and licensors. The

/* Material is protected by worldwide copyright and trade secret laws and

/* treaty provisions. No part of the Material may be used, copied, reproduced,

/* modified, published, uploaded, posted, transmitted, distributed or disclosed

/* in any way without Intel");s prior express written permission.

/* No license under any patent, copyright, trade secret or other intellectual

/* property right is granted to or conferred upon you by disclosure or delivery

/* of the Materials, either expressly, by implication, inducement, estoppel or

/* otherwise. Any license under such intellectual property rights must be

/* express and approved by Intel in writing.

/*

/*******************************************************************************

/* Content:

/* C double precision example of solving 2D Poisson problem in a

/* rectangular domain using MKL Poisson Library

/*

/*******************************************************************************/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

/* Include Poisson Library header files */

#include "mkl_poisson.h"

int main(void)

{

MKL_INT nx=6, ny=6;

double pi=3.14159265358979324;

MKL_INT ix, iy, i, stat;

MKL_INT ipar[128];

double ax, bx, ay, by, lx, ly, hx, hy, xi, yi, cx, cy, c1;

double *dpar, *f, *u, *bd_ax, *bd_bx, *bd_ay, *bd_by;

double q;

DFTI_DESCRIPTOR_HANDLE xhandle = 0;

char *BCtype;

/* Printing the header for the example */

printf("\n Example of use of MKL Poisson Library\n");

printf(" **********************************************\n\n");

printf(" This example gives the solution of 2D Poisson problem\n");

printf(" with the equation -u_xx-u_yy=f(x,y), 0<x<1, 0<y<1,\n");

printf(" f(x,y)=(8*pi*pi)*sin(2*pi*x)*sin(2*pi*y),\n");

printf(" and with the following boundary conditions:\n");

printf(" u(0,y)=u(1,y)=1 (Dirichlet boundary conditions),\n");

printf(" -u_y(x,0)=-2.0*pi*sin(2*pi*x) (Neumann boundary condition),\n");

printf(" u_y(x,1)= 2.0*pi*sin(2*pi*x) (Neumann boundary condition).\n");

printf(" -----------------------------------------------------------------------\n");

printf(" In general, the error should be of order O(1.0/nx^2+1.0/ny^2)\n");

printf(" For this example, the value of nx=ny is %d\n", nx);

printf(" The approximation error should be of order 1.0e-01, if everything is OK\n");

printf(" -----------------------------------------------------------------------\n");

printf(" DOUBLE PRECISION COMPUTATIONS \n");

printf(" =======================================================================\n\n");

dpar=(double*)malloc((13*nx/2+7)*sizeof(double));

f=(double*)malloc((nx+1)*(ny+1)*sizeof(double));

u=(double*)malloc((nx+1)*(ny+1)*sizeof(double));

bd_ax=(double*)malloc((ny+1)*sizeof(double));

bd_bx=(double*)malloc((ny+1)*sizeof(double));

bd_ay=(double*)malloc((nx+1)*sizeof(double));

bd_by=(double*)malloc((nx+1)*sizeof(double));

/* Defining the rectangular domain 0<x<1, 0<y<1 for 2D Poisson Solver */

ax=0.0E0;

bx=1.0E0;

ay=0.0E0;

by=1.0E0;

/*******************************************************************************

Setting the coefficient q to 0.

Note that this is the way to use Helmholtz Solver to solve Poisson problem!

*******************************************************************************/

q=0.0E0;

/* Computing the mesh size hx in x-direction */

lx=bx-ax;

hx=lx/nx;

/* Computing the mesh size hy in y-direction */

ly=by-ay;

hy=ly/ny;

/* Filling in the values of the TRUE solution u(x,y)=sin(2*pi*x)*sin(2*pi*y)+1

in the mesh points into the array u

Filling in the right-hand side f(x,y)=(8*pi*pi+q)*sin(2*pi*x)*sin(2*pi*y)+q

in the mesh points into the array f.

We choose the right-hand side to correspond to the TRUE solution of Poisson equation.

Here we are using the mesh sizes hx and hy computed before to compute

the coordinates (xi,yi) of the mesh points */

for(iy=0;iy<=ny;iy++)

{

for(ix=0;ix<=nx;ix++)

{

xi=hx*ix/lx;

yi=hy*iy/ly;

cx=sin(2*pi*xi);

cy=sin(2*pi*yi);

u[ix+iy*(nx+1)]=1.0E0*cx*cy;

f[ix+iy*(nx+1)]=(8.0E0*pi*pi)*u[ix+iy*(nx+1)];

u[ix+iy*(nx+1)]=u[ix+iy*(nx+1)]+1.0E0;

}

}

/* Setting the type of the boundary conditions on each side of the rectangular domain:

On the boundary laying on the line x=0(=ax) Dirichlet boundary condition will be used

On the boundary laying on the line x=1(=bx) Dirichlet boundary condition will be used

On the boundary laying on the line y=0(=ay) Neumann boundary condition will be used

On the boundary laying on the line y=1(=by) Neumann boundary condition will be used */

BCtype = "DDNN";

/* Setting the values of the boundary function G(x,y) that is equal to the TRUE solution

in the mesh points laying on Dirichlet boundaries */

for(iy=0;iy<=ny;iy++)

{

bd_ax[iy]=1.0E0;

bd_bx[iy]=1.0E0;

}

/* Setting the values of the boundary function g(x,y) that is equal to the normal derivative

of the TRUE solution in the mesh points laying on Neumann boundaries */

for(ix=0;ix<=nx;ix++)

{

bd_ay[ix]=-2.0*pi*sin(2*pi*ix/nx);

bd_by[ix]= 2.0*pi*sin(2*pi*ix/nx);

}

/* Initializing ipar array to make it free from garbage */

for(i=0;i<128;i++)

{

ipar[i]=0;

}

/* Initializing simple data structures of Poisson Library for 2D Poisson Solver */

d_init_Helmholtz_2D(&ax, &bx, &ay, &by, &nx, &ny, BCtype, &q, ipar, dpar, &stat);

if (stat!=0) goto FAILURE;

/* Initializing complex data structures of Poisson Library for 2D Poisson Solver

NOTE: Right-hand side f may be altered after the Commit step. If you want to keep it,

you should save it in another memory location! */

d_commit_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, &xhandle, ipar, dpar, &stat);

if (stat!=0) goto FAILURE;

/* Computing the approximate solution of 2D Poisson problem

NOTE: Boundary data stored in the arrays bd_ax, bd_bx, bd_ay, bd_by should not be changed

between the Commit step and the subsequent call to the Solver routine/*

Otherwise the results may be wrong. */

d_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, &xhandle, ipar, dpar, &stat);

if (stat!=0) goto FAILURE;

/* Cleaning the memory used by xhandle */

free_Helmholtz_2D(&xhandle, ipar, &stat);

if (stat!=0) goto FAILURE;

/* Now we can use xhandle to solve another 2D Poisson problem*/

/* Printing the results */

printf("The number of mesh intervals in x-direction is nx=%d\n", nx);

printf("The number of mesh intervals in y-direction is ny=%d\n\n", ny);

/* Watching the error along the line x=hx */

ix=1;

c1=0.0;

for(iy=0;iy<=ny;iy++)

{

printf("In the mesh point (%5.3f,%5.3f) the error between the computed and the true solution is equal to %10.3e\n", ix*hx, iy*hy, f[ix+iy*(nx+1)]-u[ix+iy*(nx+1)]);

if (c1<fabs(f[ix+iy*(nx+1)]-u[ix+iy*(nx+1)])) c1 = fabs(f[ix+iy*(nx+1)]-u[ix+iy*(nx+1)]);

}

if (c1>1.0e-01)

{

printf("The computed solution seems to be inaccurate. ");

goto FAILURE;

}

/* Success message to print if everything is OK */

printf("\n Double precision 2D Poisson example has successfully PASSED\n");

printf(" through all steps of computation!\n");

/* Jumping over failure message */

goto SUCCESS;

/* Failure message to print if something went wrong */

FAILURE: printf("\nDouble precision 2D Poisson example FAILED to compute the solution...\n");

SUCCESS: return 0;

/* End of the example code */

}

Fortran 90 code for the problem is presented below:

## Fortran Example to Solve 2D Poisson Problem

!*******************************************************************************

! INTEL CONFIDENTIAL

! Copyright(C) 2006-2008 Intel Corporation. All Rights Reserved.

! The source code contained or described herein and all documents related to

! the source code ("Material") are owned by Intel Corporation or its suppliers

! or licensors. Title to the Material remains with Intel Corporation or its

! suppliers and licensors. The Material contains trade secrets and proprietary

! and confidential information of Intel or its suppliers and licensors. The

! Material is protected by worldwide copyright and trade secret laws and

! treaty provisions. No part of the Material may be used, copied, reproduced,

! modified, published, uploaded, posted, transmitted, distributed or disclosed

! in any way without Intel's prior express written permission.

! No license under any patent, copyright, trade secret or other intellectual

! property right is granted to or conferred upon you by disclosure or delivery

! of the Materials, either expressly, by implication, inducement, estoppel or

! otherwise. Any license under such intellectual property rights must be

! express and approved by Intel in writing.

!

!*******************************************************************************

! Content:

! Fortran-90 double precision example of solving 2D Poisson problem in a

! rectangular domain using MKL Poisson Library

!

!*******************************************************************************

program Poisson_2D_double_precision

! Include modules defined by mkl_poisson.f90 and mkl_dfti.f90 header files

use mkl_poisson

implicit none

integer nx,ny

parameter(nx=6, ny=6)

double precision pi

parameter(pi=3.14159265358979324D0)

integer ix, iy, i, stat

integer ipar(128)

double precision ax, bx, ay, by, lx, ly, hx, hy, xi, yi, cx, cy, c1

double precision dpar(13*nx/2+7)

! Note that proper packing of data in right-hand side array f is

! automatically provided by the following declarations of the array

double precision f(nx+1,ny+1), u(nx+1,ny+1)

double precision bd_ax(ny+1), bd_bx(ny+1), bd_ay(nx+1), bd_by(nx+1)

double precision q

type(DFTI_DESCRIPTOR), pointer :: xhandle

character(4) BCtype

! Printing the header for the example

print *, ''

print *, ' Example of use of MKL Poisson Library'

print *, ' **********************************************'

print *, ''

print *, ' This example gives the solution of 2D Poisson problem'

print *, ' with the equation -u_xx-u_yy=f(x,y), 0<x<1, 0<y<1,'

print *, ' f(x,y)=(8*pi*pi)*sin(2*pi*x)*sin(2*pi*y),'

print *, ' and with the following boundary conditions:'

print *, ' u(0,y)=u(1,y)=1 (Dirichlet boundary conditions),'

print *, ' -u_y(x,0)=-2.0*pi*sin(2*pi*x) (Neumann boundary condition),'

print *, ' u_y(x,1)= 2.0*pi*sin(2*pi*x) (Neumann boundary condition).'

print *, ' -----------------------------------------------------------------------'

print *, ' In general, the error should be of order O(1.0/nx^2+1.0/ny^2)'

print '(1x,a,I1)', ' For this example, the value of nx=ny is ', nx

print *, ' The approximation error should be of order 0.1D+0, if everything is OK'

print *, ' -----------------------------------------------------------------------'

print *, ' DOUBLE PRECISION COMPUTATIONS '

print *, ' ======================================================================='

print *, ''

! Defining the rectangular domain 0<x<1, 0<y<1 for 2D Poisson Solver

ax=0.0D0

bx=1.0D0

ay=0.0D0

by=1.0D0

!*******************************************************************************

! Setting the coefficient q to 0.

! Note that this is the way to use Helmholtz Solver to solve Poisson problem!

!*******************************************************************************

q=0.0D0

! Computing the mesh size hx in x-direction

lx=bx-ax

hx=lx/nx

! Computing the mesh size hy in y-direction

ly=by-ay

hy=ly/ny

! Filling in the values of the TRUE solution u(x,y)=sin(2*pi*x)*sin(2*pi*y)+1

! in the mesh points into the array u

! Filling in the right-hand side f(x,y)=(8*pi*pi+q)*sin(2*pi*x)*sin(2*pi*y)+q

! in the mesh points into the array f.

! We choose the right-hand side to correspond to the TRUE solution of Poisson equation.

! Here we are using the mesh sizes hx and hy computed before to compute

! the coordinates (xi,yi) of the mesh points

do iy=1,ny+1

do ix=1,nx+1

xi=hx*(ix-1)/lx

yi=hy*(iy-1)/ly

cx=dsin(2*pi*xi)

cy=dsin(2*pi*yi)

u(ix,iy)=1.0D0*cx*cy

f(ix,iy)=(8.0D0*pi**2)*u(ix,iy)

u(ix,iy)=u(ix,iy)+1.0D0

enddo

enddo

! Setting the type of the boundary conditions on each side of the rectangular domain:

! On the boundary laying on the line x=0(=ax) Dirichlet boundary condition will be used

! On the boundary laying on the line x=1(=bx) Dirichlet boundary condition will be used

! On the boundary laying on the line y=0(=ay) Neumann boundary condition will be used

! On the boundary laying on the line y=1(=by) Neumann boundary condition will be used

BCtype = 'DDNN'

! Setting the values of the boundary function G(x,y) that is equal to the TRUE solution

! in the mesh points laying on Dirichlet boundaries

do iy = 1,ny+1

bd_ax(iy) = 1.0D0

bd_bx(iy) = 1.0D0

enddo

! Setting the values of the boundary function g(x,y) that is equal to the normal derivative

! of the TRUE solution in the mesh points laying on Neumann boundaries

do ix = 1,nx+1

bd_ay(ix) = -2.0*pi*dsin(2*pi*(ix-1)/nx)

bd_by(ix) = 2.0*pi*dsin(2*pi*(ix-1)/nx)

enddo

! Initializing ipar array to make it free from garbage

do i=1,128

ipar(i)=0

enddo

! Initializing simple data structures of Poisson Library for 2D Poisson Solver

call d_init_Helmholtz_2D(ax, bx, ay, by, nx, ny, BCtype, q, ipar, dpar, stat)

if (stat.ne.0) goto 999

! Initializing complex data structures of Poisson Library for 2D Poisson Solver

! NOTE: Right-hand side f may be altered after the Commit step. If you want to keep it,

! you should save it in another memory location!

call d_commit_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, xhandle, ipar, dpar, stat)

if (stat.ne.0) goto 999

! Computing the approximate solution of 2D Poisson problem

! NOTE: Boundary data stored in the arrays bd_ax, bd_bx, bd_ay, bd_by should not be changed

! between the Commit step and the subsequent call to the Solver routine!

! Otherwise the results may be wrong.

call d_Helmholtz_2D(f, bd_ax, bd_bx, bd_ay, bd_by, xhandle, ipar, dpar, stat)

if (stat.ne.0) goto 999

! Cleaning the memory used by xhandle

call free_Helmholtz_2D(xhandle, ipar, stat)

if (stat.ne.0) goto 999

! Now we can use xhandle to solve another 2D Poisson problem

! Printing the results

write(*,10) nx

write(*,11) ny

print *, ''

! Watching the error along the line x=hx

ix=2

c1 = 0.0

do iy=1,ny+1

write(*,12) (ix-1)*hx, (iy-1)*hy, f(ix,iy)-u(ix,iy)

if (dabs(f(ix,iy)-u(ix,iy)).ge.c1) c1 = dabs(f(ix,iy)-u(ix,iy))

enddo

print *, ''

if (c1.ge.0.1D+0) then

print *, 'The computed solution seems to be inaccurate.'

goto 999

endif

! Success message to print if everything is OK

print *, ' Double precision 2D Poisson example has successfully PASSED'

print *, ' through all steps of computation!'

! Jumping over failure message

go to 1

! Failure message to print if something went wrong

999 print *, 'Double precision 2D Poisson example FAILED to compute the solution...'

1 continue

10 format(1x,'The number of mesh intervals in x-direction is nx=',I1)

11 format(1x,'The number of mesh intervals in y-direction is ny=',I1)

12 format(1x,'In the mesh point (',F5.3,',',F5.3,') the error between the computed

1 and the true solution is equal to ', E10.3)

! End of the example code

end

## Spherical case

The code below computes an approximate solution of a Helmholtz problem on the entire sphere:
with aφ = 0, bφ = 2 π, aθ = 0, bθ = π, and a boundary condition
at the poles (periodic case).

The exact solution is known to be u(φ,θ) = cosθ .

Example “C Example to Solve Helmholtz Problem on a Sphere” implements the computations in C and Example “Fortran Example to Solve Helmholtz Problem on a Sphere” provides Fortran 90 code. Mind that PL interface cannot be invoked from Fortran 77 due to restrictions imposed by the use of Intel MKL FFT interface.

The algorithm of computing the approximate solution of a Poisson Problem uses Poisson Library routines. Details of the Poisson problem being solved and the errors are printed out between the computed solution and the exact one. Upon successful execution of Example “C Example to Solve Helmholtz Problem on a Sphere” the following text is printed out (Example “Fortran Example to Solve Helmholtz Problem on a Sphere” produces similar output):

Example of use of MKL Poisson Library

**********************************************

This example gives the solution of Helmholtz problem on a whole sphere

0<p<2*pi, 0<t<pi, with Helmholtz coefficient q=1 and right-hand side

f(p,t)=3*cos(t)

-----------------------------------------------------------------------

In general, the error should be of order O(1.0/np^2+1.0/nt^2)

For this example, the value of np=nt is 8

The approximation error should be of order 1.5e-01, if everything is OK

-----------------------------------------------------------------------

Note that np should be even to solve the PERIODIC problem!

-----------------------------------------------------------------------

DOUBLE PRECISION COMPUTATIONS

=======================================================================

The number of mesh intervals in phi-direction is np=8

The number of mesh intervals in theta-direction is nt=8

In the mesh point (0.785,0.000) the error between the computed and the true solution

is equal to -1.402e-001

In the mesh point (0.785,0.393) the error between the computed and the true solution

is equal to -3.097e-002

In the mesh point (0.785,0.785) the error between the computed and the true solution

is equal to -6.036e-003

In the mesh point (0.785,1.178) the error between the computed and the true solution

is equal to 2.964e-004

In the mesh point (0.785,1.571) the error between the computed and the true solution

is equal to 4.979e-017

In the mesh point (0.785,1.963) the error between the computed and the true solution

is equal to -2.964e-004

In the mesh point (0.785,2.356) the error between the computed and the true solution

is equal to 6.036e-003

In the mesh point (0.785,2.749) the error between the computed and the true solution

is equal to 3.097e-002

In the mesh point (0.785,3.142) the error between the computed and the true solution

is equal to 1.402e-001

Double precision Helmholtz example on a whole sphere has successfully PASSED

through all steps of computation!

Note that actual figures in the error may slightly differ from the figures printed above depending on the architecture and operating system used to run the example.

C code for the problem is presented below:

## C Example to Solve Helmholtz Problem on a Sphere

/*******************************************************************************

/* INTEL CONFIDENTIAL

/* Copyright(C) 2006-2008 Intel Corporation. All Rights Reserved.

/* The source code contained or described herein and all documents related to

/* the source code ("Material") are owned by Intel Corporation or its suppliers

/* or licensors. Title to the Material remains with Intel Corporation or its

/* suppliers and licensors. The Material contains trade secrets and proprietary

/* and confidential information of Intel or its suppliers and licensors. The

/* Material is protected by worldwide copyright and trade secret laws and

/* treaty provisions. No part of the Material may be used, copied, reproduced,

/* modified, published, uploaded, posted, transmitted, distributed or disclosed

/* in any way without Intel's prior express written permission.

/* No license under any patent, copyright, trade secret or other intellectual

/* property right is granted to or conferred upon you by disclosure or delivery

/* of the Materials, either expressly, by implication, inducement, estoppel or

/* otherwise. Any license under such intellectual property rights must be

/* express and approved by Intel in writing.

/*

/*******************************************************************************

/* Content:

/* C double precision example of solving Helmholtz problem on a whole sphere

/* using MKL Poisson Library

/*

/*******************************************************************************/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

/* Include Poisson Library header files */

#include "mkl_poisson.h"

int main(void)

{

/* Note that the size of the transform np must be even !!! */

MKL_INT np=8, nt=8;

double pi=3.14159265358979324;

MKL_INT ip, it, i, stat;

MKL_INT ipar[128];

double ap, bp, at, bt, lp, lt, hp, ht, theta_i, ct, c1;

double *dpar, *f, *u;

double q;

DFTI_DESCRIPTOR_HANDLE handle_s = 0;

DFTI_DESCRIPTOR_HANDLE handle_c = 0;

/* Printing the header for the example */

printf("\n Example of use of MKL Poisson Library\n");

printf(" **********************************************\n\n");

printf(" This example gives the solution of Helmholtz problem on a whole sphere\n");

printf(" 0<p<2*pi, 0<t<pi, with Helmholtz coefficient q=1 and right-hand side\n");

printf(" f(p,t)=3*cos(t)\n");

printf(" -----------------------------------------------------------------------\n");

printf(" In general, the error should be of order O(1.0/np^2+1.0/nt^2)\n");

printf(" For this example, the value of np=nt is %d\n", np);

printf(" The approximation error should be of order 1.5e-01, if everything is OK\n");

printf(" -----------------------------------------------------------------------\n");

printf(" Note that np should be even to solve the PERIODIC problem!\n");

printf(" -----------------------------------------------------------------------\n");

printf(" DOUBLE PRECISION COMPUTATIONS \n");

printf(" =======================================================================\n\n");

dpar=(double*)malloc((5*np/2+nt+10)*sizeof(double));

f=(double*)malloc((np+1)*(nt+1)*sizeof(double));

u=(double*)malloc((np+1)*(nt+1)*sizeof(double));

/* Defining the rectangular domain on a sphere 0<p<2*pi, 0<t<pi

for Helmholtz Solver on a sphere */

/* Poisson Library will automatically detect that this problem is on a whole sphere! */

ap=0.0E0;

bp=2*pi;

at=0.0E0;

bt=pi;

/* Setting the coefficient q to 1.0E0 for Helmholtz problem */

/* If you like to solve Poisson problem, please set q to 0.0E0 */

q=1.0E0;

/* Computing the mesh size hp in phi-direction */

lp=bp-ap;

hp=lp/np;

/* Computing the mesh size ht in theta-direction */

lt=bt-at;

ht=lt/nt;

/* Filling in the values of the TRUE solution u(p,t)=cos(t)

in the mesh points into the array u

Filling in the right-hand side f(p,t)=(2+q)*cos(t)

in the mesh points into the array f.

We choose the right-hand side to correspond to the TRUE solution

of Helmholtz equation on a sphere.

Here we are using the mesh sizes hp and ht computed before to compute

the coordinates (phi_i,theta_i) of the mesh points */

for(it=0;it<=nt;it++)

{

for(ip=0;ip<=np;ip++)

{

theta_i=ht*it;

ct=cos(theta_i);

u[ip+it*(np+1)]=ct;

f[ip+it*(np+1)]=ct*(2.+q);

}

}

/* Initializing ipar array to make it free from garbage */

for(i=0;i<128;i++)

{

ipar[i]=0;

}

/* Initializing simple data structures of Poisson Library

for Helmholtz Solver on a sphere */

/* As we are looking for the solution on a whole sphere, this is a PERIDOC problem */

/* Therefore, the routines ending with "_p" are used to find the solution */

d_init_sph_p(&ap,&bp,&at,&bt,&np,&nt,&q,ipar,dpar,&stat);

if (stat!=0) goto FAILURE;

/* Initializing complex data structures of Poisson Library for Helmholtz Solver on a sphere

NOTE: Right-hand side f may be altered after the Commit step. If you want to keep it,

you should save it in another memory location! */

d_commit_sph_p(f,&handle_s,&handle_c,ipar,dpar,&stat);

if (stat!=0) goto FAILURE;

/* Computing the approximate solution of Helmholtz problem on a whole sphere */

d_sph_p(f,&handle_s,&handle_c,ipar,dpar,&stat);

if (stat!=0) goto FAILURE;

/* Cleaning the memory used by handle_s and handle_c */

free_sph_p(&handle_s,&handle_c,ipar,&stat);

if (stat!=0) goto FAILURE;

/* Now we can use handle_s and handle_c to solve another Helmholtz problem */

/* after a proper initialization */

/* Printing the results */

printf("The number of mesh intervals in phi-direction is np=%d\n", np);

printf("The number of mesh intervals in theta-direction is nt=%d\n\n", nt);

/* Watching the error along the line phi=hp */

ip=1;

c1=0.0;

for(it=0;it<=nt;it++)

{

printf("In the mesh point (%5.3f,%5.3f) the error between the computed and the true solution is equal to %10.3e\n", ip*hp, it*ht, f[ip+it*(np+1)]-u[ip+it*(np+1)]);

if (c1<fabs(f[ip+it*(np+1)]-u[ip+it*(np+1)])) c1 = fabs(f[ip+it*(np+1)]-u[ip+it*(np+1)]);

}

if (c1>1.5e-01)

{

printf("The computed solution seems to be inaccurate. ");

goto FAILURE;

}

/* Success message to print if everything is OK */

printf("\n Double precision Helmholtz example on a whole sphere has successfully PASSED\n");

printf(" through all steps of computation!\n");

/* Jumping over failure message */

goto SUCCESS;

/* Failure message to print if something went wrong */

FAILURE: printf("\nDouble precision Helmholtz example on a whole sphere has FAILED to compute the solution...\n");

SUCCESS: return 0;

/* End of the example code */

}

Fortran 90 code for the problem is presented below:

## Fortran Example to Solve Helmholtz Problem on a Sphere

!*******************************************************************************

! INTEL CONFIDENTIAL

! Copyright(C) 2006-2008 Intel Corporation. All Rights Reserved.

! The source code contained or described herein and all documents related to

! the source code ("Material") are owned by Intel Corporation or its suppliers

! or licensors. Title to the Material remains with Intel Corporation or its

! suppliers and licensors. The Material contains trade secrets and proprietary

! and confidential information of Intel or its suppliers and licensors. The

! Material is protected by worldwide copyright and trade secret laws and

! treaty provisions. No part of the Material may be used, copied, reproduced,

! modified, published, uploaded, posted, transmitted, distributed or disclosed

! in any way without Intel's prior express written permission.

! No license under any patent, copyright, trade secret or other intellectual

! property right is granted to or conferred upon you by disclosure or delivery

! of the Materials, either expressly, by implication, inducement, estoppel or

! otherwise. Any license under such intellectual property rights must be

! express and approved by Intel in writing.

!

!*******************************************************************************

! Content:

! Fortran-90 precision example of solving Helmholtz problem on a whole sphere

! using MKL Poisson Library

!

!*******************************************************************************

!

program d_sph_with_poles

! Include modules defined by mkl_poisson.f90 and mkl_dfti.f90 header files

use mkl_poisson

implicit none

integer np,nt

! Note that the size of the transform np must be even !!!

parameter(np=8,nt=8)

double precision pi

parameter(pi=3.14159265358979324D0)

double precision ap,bp,hp,at,bt,ht,q,lp,lt,theta_i,ct,c1

double precision u(np+1,nt+1),f(np+1,nt+1)

type(DFTI_DESCRIPTOR), pointer :: handle_s, handle_c

integer stat

double precision dpar(5*np/2+nt+10)

integer ip,it,i

integer ipar(128)

! Printing the header for the example

print *, ''

print *, ' Example of use of MKL Poisson Library'

print *, ' **********************************************'

print *, ''

print *, ' This example gives the solution of Helmholtz problem on a whole sphere'

print *, ' 0<p<2*pi, 0<t<pi, with Helmholtz coefficient q=1 and right-hand side'

print *, ' f(p,t)=3*cos(t)'

print *, ' -----------------------------------------------------------------------'

print *, ' In general, the error should be of order O(1.0/np^2+1.0/nt^2)'

print '(1x,a,I1)', ' For this example, the value of np=nt is ', np

print *, ' The approximation error should be of order 0.15D+0, if everything is OK'

print *, ' -----------------------------------------------------------------------'

print *, ' Note that np should be even to solve the PERIODIC problem!'

print *, ' -----------------------------------------------------------------------'

print *, ' DOUBLE PRECISION COMPUTATIONS '

print *, ' ======================================================================='

print *, ''

! Defining the rectangular domain on a sphere 0<p<2*pi, 0<t<pi for Helmholtz Solver on a sphere

! Poisson Library will automatically detect that this problem is on a whole sphere!

ap=0.0D0

bp=2*pi

at=0.0D0

bt=pi

! Setting the coefficient q to 1.0D0 for Helmholtz problem

! If you like to solve Poisson problem, please set q to 0.0D0

q=1.0D0

! Computing the mesh size hp in phi-direction

lp=bp-ap

hp=lp/np

! Computing the mesh size ht in theta-direction

lt=bt-at

ht=lt/nt

! Filling in the values of the TRUE solution u(p,t)=cos(t)

! in the mesh points into the array u

! Filling in the right-hand side f(p,t)=3*cos(t)

! in the mesh points into the array f.

! We choose the right-hand side to correspond to the TRUE solution of Helmholtz equation.

! Here we are using the mesh sizes hp and ht computed before to compute

! the coordinates (phi_i,theta_i) of the mesh points

do it=1,nt+1

do ip=1,np+1

theta_i=ht*(it-1)

ct=dcos(theta_i)

u(ip,it)=ct

f(ip,it)=(2.0D0+q)*ct

enddo

enddo

! Initializing ipar array to make it free from garbage

do i=1,128

ipar(i)=0

enddo

! Initializing simple data structures of Poisson Library for Helmholtz Solver on a sphere

! As we are looking for the solution on a whole sphere, this is a PERIDOC problem

! Therefore, the routines ending with "_P" are used to find the solution

call D_INIT_SPH_P(ap,bp,at,bt,np,nt,q,ipar,dpar,stat)

if (stat.ne.0) goto 999

! Initializing complex data structures of Poisson Library for Helmholtz Solver on a sphere

! NOTE: Right-hand side f may be altered after the Commit step. If you want to keep it,

! you should save it in another memory location!

call D_COMMIT_SPH_P(f,handle_s,handle_c,ipar,dpar,stat)

if (stat.ne.0) goto 999

! Computing the approximate solution of Helmholtz problem on a whole sphere

call D_SPH_P(f,handle_s,handle_c,ipar,dpar,stat)

if (stat.ne.0) goto 999

! Cleaning the memory used by handle_s and handle_c

call FREE_SPH_P(handle_s,handle_c,ipar,stat)

if (stat.ne.0) goto 999

! Now we can use handle_s and handle_c to solve another Helmholtz problem

! after a proper initialization

! Printing the results

write(*,10) np

write(*,11) nt

print *, ''

! Watching the error along the line phi=hp

ip=2

c1 = 0.0

do it=1,nt+1

write(*,12) (ip-1)*hp, (it-1)*ht, f(ip,it)-u(ip,it)

if (dabs(f(ip,it)-u(ip,it)).ge.c1) c1 = dabs(f(ip,it)-u(ip,it))

enddo

print *, ''

if (c1.ge.0.15D+0) then

print *, 'The computed solution seems to be inaccurate.'

goto 999

endif

! Success message to print if everything is OK

print *, ' Double precision Helmholtz example on a whole sphere has successfully'

print *, ' PASSED through all steps of computation!'

! Jumping over failure message

go to 1

! Failure message to print if something went wrong

999 print *, 'Double precision Helmholtz example on a whole sphere FAILED

1 to compute the solution...'

1 continue

10 format(1x,'The number of mesh intervals in phi-direction is np=',I1)

11 format(1x,'The number of mesh intervals in theta-direction is nt=',I1)

12 format(1x,'In the mesh point (',F5.3,',',F5.3,') the error between

1 the computed and the true solution is equal to ', E10.3)

! End of the example code

end

Parent topic: PDE Support Code Examples