Solving Elliptical PDE's with gridMathematica

Richard Gass

Department of Physics
University of Cincinnati

Introduction

•	Review of gridMathematica

•	Solving Elliptical PDE's with eigenfunction expansions
	
•	Some Examples

•	Scaling and Architecture

gridMathematica

You need
• A manager kernel
• A collection of worker kernels
• And either the Parallel Computing Toolkit or Cluster Management software

Running gridMathematica

With PCT

Needs["Parallel`"];
remotenodes={"freddo.gig1.loc","romano.gig1.loc","cortado.gig1.loc","aulait.gig1.loc","conleche.gig1.loc","robusta.gig1.loc","frappuccino.gig1.loc","mezzo.gig1.loc"};
$AvailableMachines=Map[RemoteMachine,remotenodes];
LaunchSlaves[];

With Cluster Management Software

Needs["ClusterIntegration`"];
ClusterSetEngine[PBS -> {
        "ClusterEnginePath" -> "/usr/local/torque-2.1.8", "MathKernelCommand" -> "env OMP_NUM_THREADS=1 /usr/local/Wolfram/Mathematica6.0.2/OSC/math",
        "ToQueue" -> True,
        "NativeSpecification" -> {"QueueName" -> None}}
    ];
CPUs = 4;
ClusterLaunchSlaves[CPUs];

Solving Elliptical PDE's by Eigenfunction Expansion

We want to solve A ψ=λ ψ.

Expand ψ as ψ=(∑_n)(Subscript[a, n]) (Subscript[ϕ, a])where the basis functions (Subscript[ϕ, n]) form a complete orthonormal set such that <(Subscript[ϕ, n])|(Subscript[ϕ, m])> = (Subscript[δ, nm]).

The eigenvalues (Subscript[λ, n]) are given by the eigenvalues of the matrix M defined by (Subscript[M, m,n])= <(Subscript[ϕ, m])|A |(Subscript[ϕ, n]) >and the eigenfunctions are (Subscript[ψ, n])=M (Subscript[ϕ, n]).

Method is very easy to program. 
Can be slow if the matrix elements are hard to compute.
Scales badly with the number of basis functions.
Easy to run in parallel.

Example

Particle in a box with a cut-off screened Coulomb potential.

V(x)=(Subscript[V, 0]) (e^(-r/ρ))/(r+(Subscript[r, 0]))

Plot3D[V[x,y],{x,-L/2,L/2},{y,-L/2,L/2},AxesLabel->{"x","y","V(x,y)"}]

((-(ℏ^2))/(2m))(∇^2)ψ(x,y)+V(x,y)ψ(x,y)=E ψ(x,y)

Use particle in a box basis functions (Subscript[ϕ, i,j])(x,y)=(2/L) sin(i π (x+L/2)/L) sin(j π (y+L/2)/L).

Precompute the kinetic part: (Subscript[K, ijkl])=<(Subscript[ϕ, ij])|((-(ℏ^2))/(2m))(∇^2)|(Subscript[ϕ, kl])>and use NIntegrate to compute (Subscript[V, ijkl])=<(Subscript[ϕ, ij])|V(x,y)|(Subscript[ϕ, kl])>.

Code

Clear[V,ℏ,L,m,a,ρ,r0,HV,V]
ϕ[i_,j_,x_,y_]=((Sqrt[2/L]) Sin[i π( x+L/2)/L])((Sqrt[2/L]) Sin[j π( y+L/2)/L]);
HK[i_,j_,i_,j_]=Integrate[(ϕ[i,j,x,y])(-(ℏ^2)/(2m)( D[ϕ[i,j,x,y],{x,2}]+D[ϕ[i,j,x,y],{y,2}])),{x,-L/2,L/2},{y,-L/2,L/2},Assumptions->{i∈Integers,j∈Integers,L>0}];
HK[i_,j_,k_,l_]=Integrate[(ϕ[i,j,x,y])(-(ℏ^2)/(2m)( D[ϕ[k,l,x,y],{x,2}]+D[ϕ[k,l,x,y],{y,2}])),{x,-L/2,L/2},{y,-L/2,L/2},Assumptions->{i∈Integers,j∈Integers,k∈Integers,l∈Integers}];


HV[i_,j_,x_,y_]=V[x,y]ϕ[i,j,x,y];
V[x_,y_]=V0 Exp[-r/ρ]/(r+r0);
ρ=50;
r0=.5;
V0=-5;
r=(Sqrt[(x^2)+(y^2)]);
L=10;
ℏ=1;
m=1;

numberOfBasisFunctions=12;
a[i_,j_,k_,n_]:=HK[i,j,k,n]+NIntegrate[ϕ[i,j,x,y] HV[k,n,x,y],{x,-L/2,L/2},{y,-L/2,L/2},AccuracyGoal->3];
M=Table[a[i,j,l,n],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions},{l,1,numberOfBasisFunctions},{n,1,numberOfBasisFunctions}];
M=Partition[Flatten[M],(numberOfBasisFunctions^2)];
results=Chop[Eigenvalues[M], (10^(-9))];
basisFunctions=Flatten[Table[ϕ[i,j,x,y],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions}]];
wavefunctions=Eigenvectors[M].basisFunctions/Norm[Eigenvectors[M].basisFunctions];

Results and scaling

The n=16 case requires (16^4)=65536 integrals.

A more interesting case

Particle in a harmonic potential confined in an irregular region.

V[x_,y_]=(1/2)((x^2)+(y^2))+100Boole[(x^2)+Sin[x y]+ (y^2)/4>1]

Think of this as a single particle in an irregular nano-structure. Changing the factor in front of Boole allows us to control how "hard" the confinement is.

Tuning the integration

Computation of the matrix elements <(Subscript[ϕ, ij])|V(x,y)|(Subscript[ϕ, kl])>is very slow and sometimes wrong with the default method in NIntegrate.

Plot3D[ϕ[1,2,x,y]V[x,y]ϕ[3,4,x,y],{x,-L/2,L/2},{y,-L/2,L/2},PlotLabel->"Typical integrand"]

Discontinuity is over-sampled.

points=Quiet[Reap[NIntegrate[ϕ[1,2,x,y]V[x,y]ϕ[3,4,x,y],{x,-L/2,L/2},{y,-L/2,L/2},PrecisionGoal->3,EvaluationMonitor:>Sow[{x,y}]]][[2,1]]];
Graphics[{Red,Point[points]},Axes->True,AxesLabel->{"x","y"},PlotLabel->"Points sampled by NIntegrate"]

Using a Monte-Carlo method works better, but there is a better way still.

Clever Method due to Rob Knapp and Andrew Moylan

•	Convert the integral into an integral over a smooth integrand

•	Use NMaximize to find the right most edge of the region.

•	Then use FindRoot to compute the top and bottom boundaries

•	Finally use FunctionInterpolation to precompute the region 
	boundaries
	
•	We can now write the original integral as 

(∫_box)(Subscript[ϕ, ij])(Subscript[V, outside])(Subscript[ϕ, kl])dx dy+(∫_region)(Subscript[ϕ, ij])(Subscript[V, diff])(Subscript[ϕ, kl])dx dy 

where (Subscript[V, outside]) is the potential out the confining region and (Subscript[V, diff]) is the difference between the potential outside the confining region and inside the confining region.

Code

{xedge,yedge}={x,y}/.NMaximize[{y,(x^2)+Sin[x y]+ (y^2)/4==1},{x,y}][[2]];
topboundary[y_/;0≤y≤yedge]:=x/.FindRoot[(x^2)+Sin[x y]+ (y^2)/4==1,{x,0,xedge,1.5}]
topboundary[y_/;y<0]:=-bottomboundary[-y];
bottomboundary[y_/;0≤y≤yedge]:=x/.FindRoot[(x^2)+Sin[x y]+ (y^2)/4==1,{x,-1.,-1.5,xedge}]
bottomboundary[y_/;y<0]:=-topboundary[-y]
Quiet[
top=FunctionInterpolation[topboundary[y],{y,-yedge,yedge}];
bottom=FunctionInterpolation[bottomboundary[y],{y,-yedge,yedge}];
];
integral[ϕ1_,ϕ2_,options___]:=NIntegrate[ϕ1 Voutside[x,y]ϕ2,{y,-L/2,L/2},{x,-L/2,L/2},options]+NIntegrate[ϕ1 Vdiff[x,y]ϕ2,{y,-yedge,yedge},{x,bottom[y],top[y]},options]

Integration Results

•	Knapp-Moylan method is much faster. Up to 100  	times faster.

Potential Problems

The higher order terms oscillate rapidly and convergence is slow.

Many of the off-diagonal terms are zero. We will need to use the AccuracyGoal option in NIntegrate for this. The standard is AccuracyGoal->6. This is too slow and is 
overkill. Use AccuracyGoal->3.

Parallel Computing

Run on 18 processors on OSC's new supercomputer Glenn.

Set Up

Create a file called my_request.job with the following code

# PBS - l walltime = 50 : 00 : 00
# PBS - l nodes = 5 : ppn = 4
# PBS - N nef
# PBS - j oe
# PBS - S/bin/ksh
cd $Home
/usr/local/Wolfram/Mathematica6 .0 .2/OSC/math - noprompt - run "<<NEF.m"

Submission

qsub my_request.job

The file NEF.m

Needs["Parallel`Debug`"];
Needs["ClusterIntegration`"];
ClusterSetEngine[PBS -> {
        "ClusterEnginePath" -> "/usr/local/torque-2.1.8", "MathKernelCommand" -> "env OMP_NUM_THREADS=1 /usr/local/Wolfram/Mathematica6.0.2/OSC/math",
        "ToQueue" -> True,
        "NativeSpecification" -> {"QueueName" -> None}}
    ];
CPUs = 18;
ClusterLaunchSlaves[CPUs];
Clear[L];{xedge,yedge}={x,y}/.NMaximize[{y,x^2+Sin[x y]+ y^2/4==1},{x,y}][[2]];
topboundary[y_/;0<=y<=yedge]:=x/.FindRoot[x^2+Sin[x y]+ y^2/4==1,{x,0,xedge,1.5}];
topboundary[y_/;y<0]:=-bottomboundary[-y];
bottomboundary[y_/;0<=y<=yedge]:=x/.FindRoot[x^2+Sin[x y]+ y^2/4==1,{x,-1.,-1.5,xedge}];
bottomboundary[y_/;y<0]:=-topboundary[-y];
Quiet[
top=FunctionInterpolation[topboundary[y],{y,-yedge,yedge}];
bottom=FunctionInterpolation[bottomboundary[y],{y,-yedge,yedge}]
];
Voutside[x_,y_]=(1/2)(x^2+y^2)+100;
Vdiff[x_,y_]=-100;
integral[phi1_,phi2_,options___]:=NIntegrate[phi1 Voutside[x,y]phi2,{y,-L/2,L/2},{x,-L/2,L/2},options]+NIntegrate[phi1 Vdiff[x,y]phi2,{y,-yedge,yedge},{x,bottom[y],top[y]},options];
timeone=AbsoluteTime[];
phi[i_,j_,x_,y_]=((Sqrt[2/L] Sin[i Pi ( x+L/2)/L])(Sqrt[2/L] Sin[j Pi ( y+L/2)/L]));
HK[i_,j_,k_,l_]=Integrate[(phi[i,j,x,y])(-hbar^2/(2m)( D[phi[k,l,x,y],{x,2}]+D[phi[k,l,x,y],{y,2}])),{x,-L/2,L/2},{y,-L/2,L/2},Assumptions->{Element[i,Integers],Element[j,Integers],Element[k,Integers],Element[l,Integers]}];
HK[i_,j_,i_,j_]=Integrate[(phi[i,j,x,y])(-hbar^2/(2m)( D[phi[i,j,x,y],{x,2}]+D[phi[i,j,x,y],{y,2}])),{x,-L/2,L/2},{y,-L/2,L/2},Assumptions->{Element[i,Integers],Element[j,Integers],L>0}];
L=10;
hbar=1;
m=1;
numberOfBasisFunctions=CPUs;
ExportEnvironment[L,hbar,HK,phi,integral,m,Voutside,Vdiff,top,bottom,xedge,yedge,topboundary,bottomboundary,a,numberOfBasisFunctions];
a[i_,j_,k_,n_]:=HK[i,j,k,n]+integral[phi[i,j,x,y],phi[k,n,x,y],AccuracyGoal->3,PrecisionGoal->3];
M=ParallelTable[a[i,j,l,n],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions},{l,1,numberOfBasisFunctions},{n,1,numberOfBasisFunctions}];
M1=Partition[Flatten[M],numberOfBasisFunctions^2];
results=Chop[Eigenvalues[M1], 10^-9];
basisFunctions=Flatten[Table[phi[i,j,x,y],{i,1,numberOfBasisFunctions},{j,1,numberOfBasisFunctions}]];
wavefunctions=Eigenvectors[M1].basisFunctions/Norm[Eigenvectors[M1].basisFunctions];
Put[results,StringJoin["energies",ToString[Date[][[3]]],"_",ToString[Date[][[4]]],"_",ToString[Date[][[5]]]]];
wavefunctionData=Take[wavefunctions,{Length[wavefunctions]-9,Length[wavefunctions]}];
Put[wavefunctionData,StringJoin["wavefunctions",ToString[Date[][[3]]],"_",ToString[Date[][[4]]],"_",ToString[Date[][[5]]]]];
Quit[];

Results

Ground State energy is accurate to better than one percent.

Scaling and Architecture

More is not necessarily better

How much information needs to be passed between the processors? Best case is with computationally intensive tasks for each processor and little communication.

Architecture matters

• Multiple machines on a standard network
	• Dedicated cluster with Infiniband or similar 		connections
	• Cores on the same chip or motherboard

Example - Run on a quad core MacPro

Four is better than one

t1=AbsoluteTime[];
data=ParallelTable[RandomReal[],{i,1,64},{j,1,25},{k,1,25},{l,1,25}];
t2=AbsoluteTime[];
t2-t1

0.096445

t3=AbsoluteTime[];
data=Table[RandomReal[],{i,1,64},{j,1,25},{k,1,25},{l,1,25}];
t4=AbsoluteTime[];
t4-t3

0.115544

One is better than four

t1=AbsoluteTime[];
data=ParallelTable[RandomReal[],{i,1,64},{j,1,50},{k,1,50},{l,1,50}];
t2=AbsoluteTime[];
t2-t1

1.035505

t3=AbsoluteTime[];
data=Table[RandomReal[],{i,1,64},{j,1,50},{k,1,50},{l,1,50}];
t4=AbsoluteTime[];
t4-t3

0.880863

Acknowledgments

SC08 and the University of Northern Iowa
Rob Knapp and Andrew Moylan
This work was supported in part by an allocation of computing time from the Ohio Supercomputer Center