Air pollution modeling
with gridMathematica

Anton Antonov

Wolfram Research, Inc.

introduction animation

Introduction

Import["/home/usr1/antonov/MathFiles/PhDStudies/WolframTechnology2006/Pollution.JPG"]

Mathematical model

Chemical schemes

Advection simulation

Two types of parallel algorithms

Data feeding (using OOP)

Development directions

Summary

Throughout the talk we will discuss programming issues

Mathematical model

An air pollution model is given by a system of PDE's:

((∂(Subscript[c, s]))/(∂t))=-((∂(u (Subscript[c, s])))/(∂x))-((∂(v (Subscript[c, s])))/(∂y))-((∂(w (Subscript[c, s])))/(∂z))+(∂/(∂x))((Subscript[K, x])((∂(Subscript[c, s]))/(∂x)))+(∂/(∂y))((Subscript[K, y])((∂(Subscript[c, s]))/(∂y)))+(∂/(∂z))((Subscript[K, z])((∂(Subscript[c, s]))/(∂z)))+(Subscript[E, s])+(Subscript[Q, s])((Subscript[c, 1]),(Subscript[c, 2]),…,(Subscript[c, q]))-((Subscript[k, 1s])+(Subscript[k, 2s]))(Subscript[c, s]),

s=1,2,…,q.

the concentrations are denoted by (Subscript[c, s])

u,v,w and are wind velocities

(Subscript[K, x]),(Subscript[K, y]), and (Subscript[K, z]) are diffusion coefficients

the emission sources are described by (Subscript[E, s])

(Subscript[k, 1s]) and (Subscript[k, 2s]) are deposition coefficients

the chemical reactions are described by the non-linear functions Q((Subscript[c, 1]),(Subscript[c, 2]),…,(Subscript[c, q]))

the number of equations q is equal to the number of species in the model

Europe covered with steer tanks

Numerical treatment of the mathematical model

Splitting the processes

It is hard to treat the mathematical model directly. Different kinds of splitting are used.

((∂(Subscript[c, s]))/(∂t))=-((∂(u (Subscript[c, s])))/(∂x))-((∂(v (Subscript[c, s])))/(∂y)),

((∂(Subscript[c, s]))/(∂t))=(∂/(∂x))((Subscript[K, x])((∂(Subscript[c, s]))/(∂x)))+(∂/(∂y))((Subscript[K, y])((∂(Subscript[c, s]))/(∂y))),

((∂(Subscript[c, s]))/(∂t))=(Subscript[E, s])+(Subscript[Q, s])((Subscript[c, 1]),(Subscript[c, 2]),…,(Subscript[c, q])),

((∂(Subscript[c, s]))/(∂t))=-((Subscript[k, 1s])+(Subscript[k, 2s]))(Subscript[c, s]),

((∂(Subscript[c, s]))/(∂t))=-((∂(w (Subscript[c, s])))/(∂z))+(∂/(∂z))((Subscript[K, z])((∂(Subscript[c, s]))/(∂z))),

s=1,2,…,q.

Reducing to 2D

In practice often the last system is skipped. The simulations discussed in this talk use this splitting:

((∂(Subscript[c, s]))/(∂t))=-((∂(u (Subscript[c, s])))/(∂x))-((∂(v (Subscript[c, s])))/(∂y))+(∂/(∂x))((Subscript[K, x])((∂(Subscript[c, s]))/(∂x)))+(∂/(∂y))((Subscript[K, y])((∂(Subscript[c, s]))/(∂y))),

((∂(Subscript[c, s]))/(∂t))=(Subscript[E, s])+(Subscript[Q, s])((Subscript[c, 1]),(Subscript[c, 2]),…,(Subscript[c, q])),

s=1,2,…,q.

Grid interpretation

Steer tanks

grMap=Graphics[euMapScaled];
grBoxes=GraphTriagLine[elementsBoxes];
Graphics[{List@@grMap,Thickness[0.003],Red,List@@grBoxes,PointSize[0.012],Point@(Mean/@elementsBoxes)}]

Simulation rectangular mesh

grMesh=GraphTriagLine[elementsRectangular];
Graphics[{List@@grMap,Red,List@@grBoxes,PointSize[0.012],Blue,Thickness[0.003],List@@grMesh,Red,Point@(Mean/@elementsBoxes)}]

code

XSize=16;YSize=16;
mesh=Table[N[{x,y}],{x,0,1,1/XSize},{y,0,1,1/YSize}];
flatmesh=Flatten[mesh,1];
nodes=Table[i,{i,Length[Flatten[mesh,1]]}];
elements=Partition[#1,2,1]&/@Transpose[Partition[#1,2,1]&/@mesh];
elements=Flatten[Transpose[Map[Join[#1[[1]],Reverse[#1[[2]] ] ]&,elements,{2}]],1];
elementsBoxes=elements;
XSize=16;YSize=16;
mesh=Table[N[{x,y}],{x,1/XSize/2,1-1/XSize/2,1/XSize},{y,1/YSize/2,1-1/YSize/2,1/YSize}];
flatmesh=Flatten[mesh,1];
nodes=Table[i,{i,Length[Flatten[mesh,1]]}];
elements=Partition[#1,2,1]&/@Transpose[Partition[#1,2,1]&/@mesh];
elements=Flatten[Transpose[Map[Join[#1[[1]],Reverse[#1[[2]] ] ]&,elements,{2}]],1];
elementsRectangular=elements;

euMapScaled=Get["/home/antonov/MathFiles/PhDStudies/fin.el/euMapScaled.m"];

Chemical schemes

Photo-chemical mechanisms

Using ODE's

Ethene/(Subscript[NO, x]) scheme

(Subscript[NO, x])(Subscript[O, 3]) scheme

Using QSSA

With QSSA we have algebraic equations only. (No solving. Technically speaking, an O(0) approximation.)

Advection simulation

1D pseudo spectral

1D finite elements

Locally 1D finite elements implementation

2D finite elements

Using triangle finite elements

Using square finite elements

2D control volume elements

2D×1D finite elements

Time step with Crank-Nicholson

From the finite element dicretization we have the following matrix equation:

P ((∂g)/(∂t))=-A g.

The following is the Crank-Nicholson approximation:

g((Subscript[t, n])+Δt)=-(P^(-1))A g((Subscript[t, n]))+(∫_(Subscript[t, n])^(Subscript[t, n+1]))(P^(-1))A g (s)ds,

g((Subscript[t, n])+Δt)≈-(P^(-1))A g((Subscript[t, n]))+(((P^(-1))A g ((Subscript[t, n+1]))-(P^(-1))A g ((Subscript[t, n])))/(2 Δt))

From the general Crank-Nicholson form:

(Subscript[g, n+1]) = (Subscript[g, n])-Δt(θ (P^(-1)) A (Subscript[g, n+1])+(1-θ)(P^(-1)) A (Subscript[g, n])),

we derive the computational form:

(I+Δt θ (P^(-1))A)(Subscript[g, n+1])=(I-Δt(1-θ)(P^(-1))A)(Subscript[g, n]).

(How many chemical steps per advection step?)

Interpretation of the equations

G=Table[g[t,i],{i,1,Length[Nodes/.meshObject]}];

Thread[KMatrix.D[G,t]==AMatrix.G][[12]]

Thread[KMatrix.D[G,t]==AMatrix.G][[64]]//Chop

CoefficientArrays[%[[2]],{u,v,d}]//Normal

Show[ShowTriag[elements],Graphics[{PointSize[0.02],Point /@ Flatten[mesh,1],MapThread[Text[#2,#1,{-2,0}]&, {Flatten[mesh,1],nodes}]}],PlotRange->All,Axes->False]

Anisotropy of the schemes

Consider the equation

((∂c)/(∂t))=-Wcosϕ((∂c)/(∂x))-Wsinϕ((∂c)/(∂y))

The phase velocity of a wave is the rate at which the phase of the wave propagates in space.

The group velocity of a wave is the velocity with which the variations in the shape of the wave' s amplitude (modulation, envelope) propagate through space.

As Lighthill explains "... Perhaps the most fundamental property of the group velocity can be explained at once ... The energy of sinusoidal waves is propagated not at the wave speed W, but at the group velocity V".

Phase and group velocities for a square FEM grid for various CFL numbers, ((W Δt)/h)

Phase velocity log linear plot for ϕ=π/2

Group velocity log linear plot for ϕ=π/2

Two types of parallel algorithms

Parallelizing the physical domain

The domain can be partitioned among 4 processors across the x-axis in the following way:

VizualizePartition[{4,1},{64,64}]

The domains overlap. So basically, if FEM is used, we will run 4 FEM advection simulations in parallel. We may say that the parallelism is at the FEM level.

What will be the partition among 16 processors across the x-axis?

VizualizePartition[{16,1},{64,64}]

Partitioning on both axes leads to a better conditioned space domains.

VizualizePartition[{4,4},{64,64}]

Parallel advection test

What partitioning of the physical domain should be used?

-- if the grid has refinements?

-- or the processors have different computational abilities?

Parallelizing in the pseudo-space domain

It is easier to make parallel distributions of vectors (in Hilbert spaces).

If FEM is used this would mean that only one FEM is used for the advection simulation. The parallelism is pushed at the linear algebra level.

Data feeding (and OOP)

Dynamic and Static polymorphism

Design Patterns in Mathematica

The WindField object

Wind field data objects, through static inheritance and decoration.

WindField

Clear[WindField]

WindField[fileName_String,partitionSize_:32]["ReadFile"[]]:=
Module[{stream,windData,windXData,windYData},
stream=OpenRead[fileName];
windData = ReadList[stream,Number];
Close[stream];
WindField[fileName]["PartitionSize"]=partitionSize;windXData=Take[windData,{1,Length[windData]/2}];
windYData=Take[windData,{Length[windData]/2+1,-1}];
windXData=Partition[windXData,partitionSize^2];
windYData=Partition[windYData,partitionSize^2];
(*windXData=Map[Flatten[Transpose[Partition[#,32]]]&,windXData];
windYData=Map[Flatten[Transpose[Partition[#,32]]]&,windYData];*)WindField[fileName]["WindXData"]=windXData;
WindField[fileName]["WindYData"]=windYData;WindField[fileName]["WindXData"[i_Integer]]:=WindField[fileName]["WindXData"][[i]];
WindField[fileName]["WindYData"[i_Integer]]:=WindField[fileName]["WindYData"][[i]];
];

WindField[fileName_String]["SetTimeStep"[timeStepSeconds_?IntegerQ]]:=WindField[fileName]["TimeStep"]=timeStepSeconds;
WindField[fileName_String]["StartTime"]:=0;
WindField[fileName_String]["EndTime"]:=WindField[fileName]["TimeStep"]*Length[WindField[fileName]["WindXData"]];WindField[fileName_String]["SetSpaceStep"[sps_?NumberQ]]:=WindField[fileName]["SpaceStep"]=sps;
WindField[fileName_String]["SetScale"[scale_]]:=WindField[fileName]["Scale"]=scale;WindField[fileName_String]["Scale"]=(10^2);

Clear[MakeFlatMesh];
MakeFlatMesh[XSize_,YSize_,deltaX_,deltaY_]:=
Module[{X0,Y0,X1,Y1},
X0=deltaX/2;Y0=deltaY/2;
X1=X0+(XSize-1)*deltaX;Y1=Y0+(YSize-1)*deltaY;Flatten[Table[N[{x,y}],{x,X0,X1,deltaX},{y,Y0,Y1,deltaY}],1]
];
WindField[fileName_String]["FlatMesh"[dummy_]]:=WindField[fileName]["FlatMesh"[]];
WindField[fileName_String]["FlatMesh"[]]:=
Module[{deltaX},
deltaX=WindField[fileName]["SpaceStep"];
WindField[fileName]["FlatMesh"[WindField[fileName]["SpaceStep"]]]=MakeFlatMesh[32,32,deltaX,deltaX]
];

WindField[fileName_String]["WindXFunction"[index_?IntegerQ]]:=
Module[{windXData,flatmesh,intfx},
flatmesh=WindField[fileName]["FlatMesh"[WindField[fileName]["SpaceStep"]]];
windXData=WindField[fileName]@"WindXData"[index];
intfx=Interpolation[MapThread[Append[#2,#1]&,{windXData,flatmesh}],InterpolationOrder->1];WindField[fileName]["WindXFunction"[index]]=intfx
];

WindField[fileName_String]["WindYFunction"[index_?IntegerQ]]:=
Module[{windYData,flatmesh,intfy},
flatmesh=WindField[fileName]["FlatMesh"[WindField[fileName]["SpaceStep"]]];
windYData=WindField[fileName]@"WindYData"[index];
intfy=Interpolation[MapThread[Append[#2,#1]&,{windYData,flatmesh}],InterpolationOrder->1];WindField[fileName]["WindYFunction"[index]]=intfy
];

WindField[fileName_String]["Interpolate"[index_?IntegerQ]]:=
{WindField[fileName]["WindXFunction"[index]],WindField[fileName]["WindYFunction"[index]]};

WindField[fileName_String]["InterpolateField"[time_,xsize_?IntegerQ,ysize_?IntegerQ]]:=
Module[{i,θ,ts,sps,intfx0,intfy0,intfx1,intfy1,flatmesh,partitionSize},
ts=WindField[fileName]["TimeStep"];
sps=WindField[fileName]["SpaceStep"];
i=Floor[time/ts];
θ=((time-ts*i)/ts)//N;
intfx0=WindField[fileName]["WindXFunction"[i+1]];
intfy0=WindField[fileName]["WindYFunction"[i+1]];
intfx1=WindField[fileName]["WindXFunction"[i+2]];
intfy1=WindField[fileName]["WindYFunction"[i+2]];
partitionSize=WindField[fileName]["PartitionSize"];
flatmesh=MakeFlatMesh[xsize,ysize,partitionSize*sps/xsize,partitionSize*sps/ysize];
{(θ#2+(1-θ)#1)&[intfx0@@@flatmesh,intfx1@@@flatmesh],
(θ#2+(1-θ)#1)&[intfy0@@@flatmesh,intfy1@@@flatmesh]}
]/;0≤time<=WindField[fileName]["EndTime"];

Clear[MakeElementWindField];
MakeElementWindField[windField_]:=
Table[Mean[{windField[[i,j]]+windField[[i+1,j]]+windField[[i,j+1]]+windField[[i+1,j+1]]}],{i,1,Length[windField]-1},{j,1,Length[windField[[1]]]-1}]/;ArrayQ[windField];

WindFieldEnhanced

Clear[WindFieldEnhanced]

(* static inheritance *)
SubValues[WindFieldEnhanced]=ReplaceRepeated[SubValues[WindField],WindField->WindFieldEnhanced];

Clear[ReplaceWithCopies];
ReplaceWithCopies[field_?MatrixQ,factor_?IntegerQ]:=
Module[{t},
t=Map[Flatten[Map[Table[#,{factor}]&,#],1]&,field];
Flatten[Map[Table[#,{factor}]&,t],1]
];

WindFieldEnhanced[fileName_String]["Enhance"[index_?IntegerQ,factor_?IntegerQ]]:=
{Flatten[ReplaceWithCopies[Partition[WindFieldEnhanced[fileName]["WindXData"[index]],WindFieldEnhanced[fileName]["PartitionSize"]],factor],1],Flatten[ReplaceWithCopies[Partition[WindFieldEnhanced[fileName]["WindYData"[index]],WindFieldEnhanced[fileName]["PartitionSize"]],factor],1]};

WindFieldEnhanced[fileName_String]["InterpolateEnhancedField"[time_,factor_?IntegerQ]]:=
Module[{i,θ,ts,sps,wfx0,wfy0,wfx1,wfy1},
ts=WindFieldEnhanced[fileName]["TimeStep"];
sps=WindFieldEnhanced[fileName]["SpaceStep"];
i=Floor[time/ts];
θ=((time-ts*i)/ts)//N;
{wfx0,wfy0}=WindFieldEnhanced[fileName]["Enhance"[i+1,factor]];
{wfx1,wfy1}=WindFieldEnhanced[fileName]["Enhance"[i+2,factor]];
{θ wfx1+(1-θ)wfx0,θ wfy1+(1-θ)wfy0}
]/;0≤time<=WindFieldEnhanced[fileName]["EndTime"];

(* position of the overridden method *)pos=Position[SubValues[WindFieldEnhanced],WindFieldEnhanced[___]["InterpolateField"[___]],∞][[1,1]];
SubValues[WindFieldEnhanced]=Drop[SubValues[WindFieldEnhanced],{pos}];
WindFieldEnhanced[fileName_String]["InterpolateField"[time_,xsize_?IntegerQ,ysize_?IntegerQ]]:=
With[{partitionSize=WindFieldEnhanced[fileName]["PartitionSize"]},
If[IntegerQ[xsize/partitionSize]&&IntegerQ[ysize/partitionSize]&&xsize==ysize,
WindFieldEnhanced[fileName]["InterpolateEnhancedField"[time,xsize/partitionSize]]
]
];

wf=WindFieldEnhanced["/home/antonov/OODEM/data/metdata/data1995/jul95winSpaced"];
wf@"ReadFile"[]
wf@"PartitionSize"
wf@"SetTimeStep"[21600]
wf@"SetSpaceStep"[200*10^3]

wf@"PartitionSize"

{windx,windy}=wf@"InterpolateEnhancedField"[1,2];

{windx,windy}=wf@"InterpolateField"[100,64,64];

WindFieldSubRegion

Clear[WindFieldSubRegion];

WindFieldSubRegion[wfObj:_Symbol[_String],sc:{{_?IntegerQ,_?IntegerQ},{_?IntegerQ,_?IntegerQ}}]["Create"[]]:=
Block[{},
(* delegation *)
WindFieldSubRegion[wfObj][pr_]:=wfObj[pr];WindFieldSubRegion[wfObj][mtd_[args___]]:=wfObj[mtd[args]];
(* subregion coodinates *)
WindFieldSubRegion[wfObj]["SubregionCoordinates"]=sc;
WindFieldSubRegion[wfObj]["Sizes"]={sc[[1,2]]-sc[[1,1]]+1,sc[[2,2]]-sc[[2,1]]+1};
WindFieldSubRegion[wfObj]
]/;(Head[wfObj]===WindFieldEnhanced||Head[wfObj]===WindField);

WindFieldSubRegion[wfObj_]["InterpolateField"[time_,xsize_?IntegerQ,ysize_?IntegerQ]]:=
Module[{windx,windy,x0,x1,y0,y1},
{windx,windy}=wfObj["InterpolateField"[time,xsize,ysize]];
{{x0,x1},{y0,y1}}=WindFieldSubRegion[wfObj]["SubregionCoordinates"];
windx=Take[#,{y0,y1}]&/@Take[Partition[windx,xsize],{x0,x1}];windy=Take[#,{y0,y1}]&/@Take[Partition[windy,xsize],{x0,x1}];
{Flatten[windx],Flatten[windy]}
];

WindField object Test

wf=WindFieldEnhanced["/home/antonov/OODEM/data/metdata/data1995/jul95winSpaced"];
wf@"ReadFile"[]
wf@"PartitionSize"
wf@"SetTimeStep"[21600]
wf@"SetSpaceStep"[200*10^3]

wfsr=WindFieldSubRegion[wf,{{31,64},{1,33}}]@"Create"[]

{{x0,x1},{y0,y1}}=wfsr@"SubregionCoordinates"

{windx64,windy64}=wf@"InterpolateField"[100,64,64];

{windx,windy}=wfsr@"InterpolateField"[100,64,64];

flatmesh=MakeFlatMesh[64,64,100*10^3,100*10^3];
wfield=wf@"InterpolateField"[100,64,64];
wfield*=10^2;
wfield=MapThread[Line[{#1,#1+{#2,#3}}]&,{flatmesh,wfield[[1]],wfield[[2]]}];Graphics[Partition[wfield,64],PlotRange->All,AspectRatio->Automatic,GridLines->{Table[i 100*10^3,{i,0,64,8}],Table[i 100*10^3,{i,0,64,8}]},Axes->True]

flatmeshS=Flatten[Take[#,{y0,y1}]&/@Take[Partition[flatmesh,64],{x0,x1}],1];
wfield=wfsr@"InterpolateField"[100,64,64];
wfield*=10^2;
wfieldArr=MapThread[Line[{#1,#1+{#2,#3}}]&,{flatmeshS,wfield[[1]],wfield[[2]]}];Graphics[Partition[wfieldArr,x1-x0+1],PlotRange->All,AspectRatio->Automatic,GridLines->{Table[i 100*10^3,{i,0,64,8}],Table[i 100*10^3,{i,0,64,8}]},Axes->True]

The AdvectionOperator object

Instead of assembling the advection matrix at each time step, we can derive as a linear combination of matrices assembled at the data time points.

A((Subscript[t, n]))==θ A((Subscript[t, (k+1) m]))+(1-θ)A((Subscript[t, k m])),

k m≤n≤(k+1)m,

θ=(((Subscript[t, n])-(Subscript[t, k m]))/((Subscript[t, (k+1)m])-(Subscript[t, k m]))).

AdvectionOperator

Clear[AdvectionOperator];
AdvectionOperator[wf_]["MakeOperator"[time_]]:=
Module[{windx,windy,windXOnElements,windYOnElements,rules,totalXSize,totalYSize,xsize,ysize,t},
(* approximate the wind field for the current time *)
{totalXSize,totalYSize}=AdvectionOperator[wf]["Sizes"];
{windx,windy}=wf@"InterpolateField"[time,totalXSize,totalYSize];
t=wf["Sizes"];
If[MatchQ[t,{_?NumberQ,_?NumberQ}],
{xsize,ysize}=t,
{xsize,ysize}={totalXSize,totalYSize}
];
windXOnElements=Flatten[MakeElementWindField[Partition[windx,ysize]],1];
windYOnElements=Flatten[MakeElementWindField[Partition[windy,ysize]],1];

(* assemble the advection matrix *)
rules=Join@@(AKMatrixRules[meshObject,#,windXOnElements[[#]]&,windYOnElements[[#]]&]&/@(MeshOrbits/.meshObject));
AMatrix=SparseArray[rules]
];
AdvectionOperator[wf_]["Initialize"[xsize_:XSize,ysize_:YSize]]:=
Module[{},
AdvectionOperator[wf]["Sizes"]={xsize,ysize};
AdvectionOperator[wf]["BeginTime"]=0;
AdvectionOperator[wf]["EndTime"]=wf@"TimeStep";
AdvectionOperator[wf]["A0"]=AdvectionOperator[wf]["MakeOperator"[0]];
AdvectionOperator[wf]["A1"]=AdvectionOperator[wf]["MakeOperator"[wf@"TimeStep"]];
];
AdvectionOperator[wf_]["Approximate"[time_]]:=
Module[{θ},
If[AdvectionOperator[wf]["EndTime"]<time,
AdvectionOperator[wf]["BeginTime"]=AdvectionOperator[wf]["EndTime"];AdvectionOperator[wf]["EndTime"]=AdvectionOperator[wf]["EndTime"]+wf["TimeStep"];
AdvectionOperator[wf]["A0"]=AdvectionOperator[wf]["A1"];
AdvectionOperator[wf]["A1"]=AdvectionOperator[wf]["MakeOperator"[AdvectionOperator[wf]["EndTime"]]];
];
θ=(time-AdvectionOperator[wf]["BeginTime"])/wf["TimeStep"];
Which[
θ==0,AdvectionOperator[wf]["A0"],
θ==1,AdvectionOperator[wf]["A1"],
True,
θ AdvectionOperator[wf]["A1"]+(1-θ)AdvectionOperator[wf]["A0"]
]
];

AdvectionOperator object Test

WindField["/home/antonov/OODEM/data/metdata/data1995/jul95winSpaced"]
wf@"ReadFile"[]
wf@"SetTimeStep"[21600]
wf@"SetSpaceStep"[200*10^3]
advop = AdvectionOperator[wf];
advop@"Initialize"[64,64]
AMatrix = advop@"Approximate"[100]
AAMatrix = advop@"MakeOperator"[100]; // Timing
AAMatrix == AMatrix

$MessagePrePrint=.

Parallel runs

Implementation studies

Parallel run

Profiles

3 pollutants for 100 steps  8 x1 

with 96 x96 : 
Advection : 1.5 seconds
Exchange : 22 seconds
Chemistry : 89 seconds
Assembly : 4 seconds

with 192 x192 : 
Advection : 7 seconds
Exchange : 40 seconds
Chemistry : 345 seconds
Assembly : 13 seconds

IDF demo profiles

PSC1 - 4 slaves, for a total of 32 cores.
Each machine was linked via Mellanox Infiniband.

8 cores : 1 - step : 29 seconds
10 - steps : 515 seconds

16 cores : 1 - step : 17 seconds
10 - steps : 184 seconds

32 cores : 1 - step : 10 seconds
10 - steps : 105 seconds

Development directions

Larger chemical schemes

3D advection

Grid refinement

Summary

For these project were used the following

Piecewise with NIntegrate/Integrate,
NDSolve and its framework,
SparseArray's,
LinearSolve and its preconditioners, 
OOP with Downvalues/Subvalues, 
Parallel Computational Toolkit.

Acknowledgements: Schoeller Porter, Manjula Jayasuriya.