3D Pollution Modeling with gridMathematica

Anton Antonov

Wolfram Research, Inc.

Talk topics

■ Mathematical formulation

■ Advection modeling

■ Chemistry modeling

■ Parallel implementation

■ Data feeding

■ Visualization

Mathematical Modeling

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.

The advection sub-model is with open boundary conditions.

Advection modeling -- 3D Grid Generation

It is often desirable to do some sort of numerical zooming over certain areas of interest, since increasing the resolution over the whole area can make the simulation prohibitively large.

We can use Mathematica''s graphics function DensityPlot to create 2D grids with local refinements. From the structured data of GraphicsComplex objects we can construct a layered 3D grid.

Consider the function:

n=12;nd=12;
data=Flatten[Table[{i,j,Piecewise[{{5+(Cos[i π]+Sin[j π]+Cos[(i+1/2)π]+Sin[(j+1/2)π]),1/4n<i<3/4n&&1/4n<j<3/4n}},0]},{i,0,n},{j,0,n}],1];
ifunc=Interpolation[data];
Manipulate[Plot3D[ifunc[i,j],{i,0,n},{j,0,n},PlotPoints->nd,AspectRatio->Automatic,PlotRange->All,Mesh->mesh,MaxRecursion->maxRecursion],{mesh,{Automatic,All}},{maxRecursion,{Automatic,1}},SaveDefinitions->True]

In order to plot that function as a density, DensityPlot needs to construct a grid with refinement

Manipulate[gr=DensityPlot[ifunc[i,j],{i,0,n},{j,0,n},PlotPoints->nd+1,AspectRatio->Automatic,Mesh->mesh,ClippingStyle->Automatic,MaxRecursion->maxRecursion,ColorFunction ->colorFunction],{colorFunction,{Automatic,(1&)}},{mesh,{False,All}},{maxRecursion,{Automatic,1}},SaveDefinitions->True]

Below is the code that constructs a 3D layered grid construction:

Clear[AddZLevels];
AddZLevels[points_,triangles_,zLevels_]:=
Block[{points3D,triangles3D},
points3D=Map[Function[{z},Append[#,z]&/@points],zLevels];
triangles3D=NestList[#+Length[points]&,triangles,Length[zLevels]-1];
{points3D,triangles3D}
];
Clear[MakeIndexPrisms];
MakeIndexPrisms[triangles3D_]:=
MapThread[Join[#1,#2]&,#]&/@Partition[triangles3D,2,1];
Clear[GetPrism];
GetPrism[points3D_,prisms_,i_]:=Part[points3D,prisms[[i]]];
Clear[PrismGraph];
PrismGraph[prism_]:={Line[Append[Take[prism,3],prism[[1]]]],Line[Append[Take[prism,-3],prism[[-3]]]],MapThread[Line[{#1,#2}]&,{Take[prism,3],Take[prism,-3]}]};
Clear[CoordsLessQ];
CoordsLessQ[c1_, c2_] :=
  Which[
    c1[[1]] < c2[[1]], True,
    c1[[1]] > c2[[1]], False,
    True, c1[[2]] <= c2[[2]]
    ];
gr=DensityPlot[ifunc[i,j],{i,0,n},{j,0,n},PlotPoints->nd+1,AspectRatio->Automatic,Mesh->All,ClippingStyle->Automatic,MaxRecursion->1,ColorFunction ->(1&)];
grC=gr[[1]];
points=grC[[1]];
triangles=Cases[grC[[2]],Polygon[s_]:>s,∞][[1]];
pointsSorted=MapIndexed[{#1,#2[[1]]}&,points];
pointsSorted=Sort[pointsSorted,CoordsLessQ[#1[[1]],#2[[1]]]&];
points=pointsSorted[[All,1]];
indexesSorted=MapIndexed[{#1,#2[[1]]}&,pointsSorted[[All,2]]];
indexesSorted=Sort[indexesSorted,#1[[1]]<#2[[2]]&];
triangles=triangles/.(Rule@@@indexesSorted);
zLevels=Range[1,3];
{points3D,triangles3D}=AddZLevels[points,triangles,zLevels];
prisms=Flatten[MakeIndexPrisms[triangles3D],1];
points3D=Flatten[points3D,1];
Graphics3D[{GraphicsComplex[points3D,PrismGraph/@prisms]}]

Finite element integration

◼ Partition the 2D grid into classes of triangles with the relation defined as: (Subscript[△, 1])is in the same class with (Subscript[△, 2]) iff there is a vector v and scalar α such that v+(Subscript[α△, 1])=(Subscript[△, 2]). (Note, no rotations.)

◼ For each class of triangles find the FEM integrals for the prism that uses a triangle of this class as base.

PrismFEM.m

Time step and operator assembly

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]).

That last formula is computed using SparseArray's and LinearSolve.

(How many chemical steps per advection step?)

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]))).

Chemistry modeling -- single cell

From chemical equations to ODE's

reactions=
{"(Subscript[NO, 2])"+hν->"NO"+"O",
"O"+"(Subscript[O, 2])"+"M"->"(Subscript[O, 3])"+"M",
"(Subscript[O, 3])"+"NO"->"(Subscript[NO, 2])"+"(Subscript[O, 2])",
"RH"+"OH"->"(Subscript[RO, 2])"+"(Subscript[H, 2])O",
"HCHO"+hν->2"(Subscript[HO, 2])"+"CO",
"HCHO"+hν->"(Subscript[H, 2])"+"CO",
"HCHO"+"OH"->"(Subscript[HO, 2])"+"CO"+"(Subscript[H, 2])O",
"RCHO"+"OH"->"RC(O)(Subscript[O, 2])"+"(Subscript[H, 2])O",
"(Subscript[RO, 2])"+"NO"->"(Subscript[NO, 2])"+"RO",
"RC(O)(Subscript[O, 2])"+"NO"->"(Subscript[NO, 2])"+"(Subscript[RO, 2])"+"(Subscript[CO, 2])",
"RO"+"(Subscript[O, 2])"->"RCHO"+"(Subscript[HO, 2])",
"(Subscript[HO, 2])"+"NO"->"(Subscript[NO, 2])"+"OH",
"(Subscript[O, 3])"+hν->"(Subscript[O, 2])"+"O((()^1)D)",
"O((()^1)D)"+"(Subscript[H, 2])O"->2 "OH",
"OH"+"(Subscript[NO, 2])"+"M"->"(Subscript[HNO, 3])"+"M",
"(Subscript[HO, 2])"+"(Subscript[HO, 2])"->"(Subscript[H, 2])(Subscript[O, 2])"+"(Subscript[O, 2])",
"(Subscript[RO, 2])"+"(Subscript[HO, 2])"->"ROOH"+"(Subscript[O, 2])",
"RC(O)(Subscript[O, 2])"+"(Subscript[NO, 2])"+"M"->"PAN"+"M",
"PAN"->"RC(O)(Subscript[O, 2])"+"(Subscript[NO, 2])",
"(Subscript[HO, 2])"+"(Subscript[O, 3])"->"OH"+2"(Subscript[O, 2])"
};
reactions//ColumnForm

Format[c[n_,t___]]:=Subscript[c,n][t];
Format[Derivative[0,1][c][n_,t_,r___]]:=D[Subscript[c,n][t,r],t];
Format[k[n_,t___]]:=Subscript[k,n][t];
MapAt[ColumnForm,FormODE["(Subscript[O, 3])",reactions,c,{t,T,z}],2]//TraditionalForm

MapAt[ColumnForm,FormODE["OH",reactions,c,{t,T,z}],2]//TraditionalForm

Simulation with NDSolve

Get["~/MathFiles/DemoDVD/GridDemos/PollutionSimulation3D/GeneralizedMechanism.m"];
ic={0,0,2.463`*^12,0,0,2.463`*^12,2.463`*^11,0,0,5.4186`*^18,0,0,0,2.463`*^12,0,2.463`*^12,0,0,0};
nmin=400;
Manipulate[
Module[{timing,res,timePoints,state,pllabel,scalepoints},
Δt = timeStep;
timing=Timing[res = FoldList[
(state = MakeChemistryState[50, 298, 1, #2];{#2,SolveChemistry[state, #1[[2]], 50, 298, 1, Δt]}) & ,{0, ic}, 
      Table[t, {t, Δt, nmin*60, Δt}]];]//First;
timePoints=First/@res;
res=Transpose[#[[2]]&/@res];res=res/2.643*^13 ;(*convert to ppb*)
res=Map[Transpose[{timePoints,#1}]&,res];pllabel="time with Δt="<>ToString[Δt]<>" is "<>ToString[timing];
scalepoints=Table[If[Mod[i,50*60]==0,{i,(i)/60},{i,""}],{i,0,nmin*60,20*60}];ListLinePlot[MapThread[Tooltip[#1,#2]&,{res,species1}],FrameLabel->{"time,min","ppb"},Frame->True,FrameTicks->{scalepoints,Automatic},PlotRange->{{0,nmin*60},Automatic},GridLines->{First/@scalepoints,Automatic},PlotLabel->pllabel]
],{timeStep,{100,600,900}}]

QSSA equations derivation

Because of the large number of ODE systems that need to be solved the method Quasi Steady-State Approximation (QSSA) is often used. With QSSA instead of using ODE's solvers we do arithmetic assignments.

Consider the general form of an ODE derived from a chemical kinetic mechanism:

((∂(Subscript[c, s])(t))/(∂t))==(Subscript[P, s])((Subscript[c, 1]),(Subscript[c, 2]),…,(Subscript[c, n]))-(Subscript[Q, s])((Subscript[c, 1]),(Subscript[c, 2]),…,(Subscript[c, n])) (Subscript[c, s])(t), i=1…q

These equations have exponential solutions provided the (Subscript[P, s]) and (Subscript[Q, s]) are constants:

DSolve[{D[c[t],t]==P+Q c[t],c[t0]==c0},c[t],t][[1,1]]/.{c0->c[t0],t->t0+h}/.{t0->(Subscript[t, 0]),c->(Subscript[c, s])}//FullSimplify//Apart//TraditionalForm

Replacing all the chemical kinetics ODE's with the solution above is the basis of QSSA. (Also called Plain QSSA). Further improvements of this method can be found in:

L.O.Jay, A.Sandu, F.A.Potra, and G.R.Carmichael. Improved QSSA methods for atmospheric chemistry integration.Technical report 67, The University of Iowa, Department of Mathematics, Iowa City, February 1995. http : // citeseer.ist.psu.edu/article/jay95improved.html

Here we compare Plain QSSA with NDSolve. Plain QSSA is fast if it uses the same time step sizes as NDSolve based integration, but in order to get satisfactory result the time step need to be decreased which makes Plain QSSA around 5-10 times slower than using NDSolve.

Manipulate[
Module[{timing,res,timePoints,state,pllabel,scalepoints,factor},
Δtqssa=timeStep;
factor=(Δt/Δtqssa);timing=Quiet[Timing[res = FoldList[{#2,SolveChemistryQSSA[QSSACoeff, 
                #1[[2]], 50, 298, 1, #2, Δtqssa]} & ,{0, ic}, 
            Table[t, {t, Δtqssa, nmin*60, Δtqssa}]];]//First];
timePoints=First/@res;
res=Transpose[#[[2]]&/@res];
res=res/2.643*^13 ;(*convert to ppb*)
res=Map[Transpose[{timePoints,#1}]&,res];pllabel="time with Δt="<>ToString[Δt]<>" is "<>ToString[timing];
scalepoints=Table[If[Mod[i,50*60]==0,{i,(i)/60},{i,""}],{i,0,nmin*60,20*60}];ListLinePlot[MapThread[Tooltip[#1,#2]&,{res,species1}],FrameLabel->{"time,min","ppb"},Frame->True,FrameTicks->{scalepoints,Automatic},PlotRange->{{0,nmin*60},Automatic},GridLines->{First/@scalepoints,Automatic},PlotLabel->pllabel]
]
,{{timeStep,Δt},{10,20,50,100,Δt}}]

Chemistry modeling

The chemical reaction rates depend upon sun linght, temperature and altitude.

Using NDSolves framework we pre-initialize the equations for each layer of the grid.

The chemistry simulation is embarassingly parallel.

Parallel implementation

Matrix algorithms and PDE solvers in general (as in CFD) can be programmed to be parallel over a Cartesian grid of processors.

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?

Distribution of the locally refined grid

Collecting the grids on the master processor:

Graphics[RemoteEvaluate[Graphics[Drop[totalGridGR,-1]],chanels[[2]]],Axes -> True]

Distribution over eight processors:

RemoteEvaluate[Graphics[grC,Axes->True]]

On each processor is then constructed the 3D grid as described above, and each processor initializes and runs its own Finite Element Method.

Approximation of the wind field

Approximation of wind field: given wind speeds on a non-uniformly distributed points approximate them on a regular Cartesian grid and locally refined grid. This can be done using ListPlot3D.

The red dots are represent the x,y-positions and the values of, say, the x-component of the wind field. The black dots are the wind field approximation over a regular grid.

Wind3DApproximation.nb

Emissions data handling

Given the average emissions over the Illinois counties approximate them on a locally refined grid.

We have imported the county shape data from US county boundaries site: www.nationalatlas.gov , 
http://edcftp.cr.usgs.gov/pub/data/nationalatlas/countyp020.tar.gz

IllinoisCounties.nb

Summary

For this project were used the following Mathematica functions and data structures:

◼ Piecewise and DensityPlot's adaptive mesh generation.

◼ Piecewise with NIntegrate/Integrate,

◼ NDSolve and its framework,

◼ SparseArray's,

◼ LinearSolve and its preconditioners,

◼ GridLines of ListPlot3D.

Implementation was done using:

◼ OOP with Downvalues/SubValues,

◼ gridMathematica.

Acknowledgements: Schoeller Porter, Manjula Jayasuriya.