Skip to content

Commit b9a0331

Browse files
ahmadiaahmadia
committed
Merge pull request #2 from ahmadia/sc_2013
[WIP] Distributed Computing, Sparse Linear Algebra, and Anatomy of PyClaw
2 parents 3ba2473 + 21b8f9a commit b9a0331

17 files changed

+2641
-0
lines changed

‎3.1_Distributed_Computing.ipynb‎

Lines changed: 1359 additions & 0 deletions
Large diffs are not rendered by default.

‎petsc4py-tutorial/Bratu3D.pyx‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
from petsc4py.PETSc cimport Vec, PetscVec
2+
from petsc4py.PETSc cimport Mat, PetscMat
3+
from petsc4py.PETSc cimport DA, PetscDM
4+
from petsc4py.PETSc cimport SNES, PetscSNES
5+
6+
from petsc4py.PETSc import Error
7+
8+
cdef extern from "Bratu3Dimpl.h":
9+
ctypedef struct Params:
10+
double lambda_
11+
int FormInitGuess(PetscDM da, PetscVec x, Params *p)
12+
int FormFunction (PetscDM da, PetscVec x, PetscVec F, Params *p)
13+
int FormJacobian (PetscDM da, PetscVec x, PetscMat J, Params *p)
14+
15+
def formInitGuess(Vec x, DA da, double lambda_):
16+
cdef int ierr
17+
cdef Params p = {"lambda_" : lambda_}
18+
ierr = FormInitGuess(da.dm, x.vec, &p)
19+
if ierr != 0: raise Error(ierr)
20+
21+
def formFunction(SNES snes, Vec x, Vec f, DA da, double lambda_):
22+
cdef int ierr
23+
cdef Params p = {"lambda_" : lambda_}
24+
ierr = FormFunction(da.dm, x.vec, f.vec, &p)
25+
if ierr != 0: raise Error(ierr)
26+
27+
def formJacobian(SNES snes, Vec x, Mat J, Mat P, DA da, double lambda_):
28+
cdef int ierr
29+
cdef Params p = {"lambda_" : lambda_}
30+
ierr = FormJacobian(da.dm, x.vec, P.mat, &p)
31+
if ierr != 0: raise Error(ierr)
32+
if J != P: J.assemble() # for matrix-free operator
33+
return Mat.Structure.SAME_NONZERO_PATTERN

‎petsc4py-tutorial/Bratu3Dimpl.c‎

Lines changed: 301 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,301 @@
1+
/* ------------------------------------------------------------------------
2+
3+
Solid Fuel Ignition (SFI) problem. This problem is modeled by the
4+
partial differential equation
5+
6+
-Laplacian(u) - lambda * exp(u) = 0, 0 < x,y,z < 1,
7+
8+
with boundary conditions
9+
10+
u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
11+
12+
A finite difference approximation with the usual 7-point stencil
13+
is used to discretize the boundary value problem to obtain a
14+
nonlinear system of equations. The problem is solved in a 3D
15+
rectangular domain, using distributed arrays (DAs) to partition
16+
the parallel grid.
17+
18+
------------------------------------------------------------------------- */
19+
20+
#include "Bratu3Dimpl.h"
21+
22+
#if PETSC_VERSION_(3,1,0)
23+
#define DMDAGetInfo(da,dim,M,N,P,m,n,p,dof,s,bx,by,bz,st) \
24+
DAGetInfo((DA)da,dim,M,N,P,m,n,p,dof,s,bx,st)
25+
#define DMDAGetCorners(da,x,y,z,m,n,p) DAGetCorners((DA)da,x,y,z,m,n,p)
26+
#define DMDAVecGetArray(da,v,a) DAVecGetArray((DA)da,v,a)
27+
#define DMDAVecRestoreArray(da,v,a) DAVecRestoreArray((DA)da,v,a)
28+
#endif
29+
30+
#undef __FUNCT__
31+
#define __FUNCT__ "FormInitGuess"
32+
PetscErrorCode FormInitGuess(DM da, Vec X, Params *p)
33+
{
34+
PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
35+
PetscReal lambda,temp1,hx,hy,hz,tempk,tempj;
36+
PetscScalar ***x;
37+
PetscErrorCode ierr;
38+
39+
PetscFunctionBegin;
40+
41+
ierr = DMDAGetInfo(da,PETSC_IGNORE,
42+
&Mx,&My,&Mz,
43+
PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
44+
PETSC_IGNORE,PETSC_IGNORE,
45+
PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
46+
PETSC_IGNORE);
47+
48+
lambda = p->lambda_;
49+
hx = 1.0/(PetscReal)(Mx-1);
50+
hy = 1.0/(PetscReal)(My-1);
51+
hz = 1.0/(PetscReal)(Mz-1);
52+
temp1 = lambda/(lambda + 1.0);
53+
54+
/*
55+
Get a pointer to vector data.
56+
57+
- For default PETSc vectors, VecGetArray() returns a pointer to
58+
the data array. Otherwise, the routine is implementation
59+
dependent.
60+
61+
- You MUST call VecRestoreArray() when you no longer need access
62+
to the array.
63+
*/
64+
ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr);
65+
66+
/*
67+
Get local grid boundaries (for 3-dimensional DA):
68+
69+
- xs, ys, zs: starting grid indices (no ghost points)
70+
71+
- xm, ym, zm: widths of local grid (no ghost points)
72+
*/
73+
ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
74+
75+
/*
76+
Compute initial guess over the locally owned part of the grid
77+
*/
78+
for (k=zs; k<zs+zm; k++) {
79+
tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
80+
for (j=ys; j<ys+ym; j++) {
81+
tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
82+
for (i=xs; i<xs+xm; i++) {
83+
if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
84+
/* boundary conditions are all zero Dirichlet */
85+
x[k][j][i] = 0.0;
86+
} else {
87+
x[k][j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
88+
}
89+
}
90+
}
91+
}
92+
93+
/*
94+
Restore vector
95+
*/
96+
ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr);
97+
98+
PetscFunctionReturn(0);
99+
}
100+
101+
102+
#undef __FUNCT__
103+
#define __FUNCT__ "FormFunction"
104+
PetscErrorCode FormFunction(DM da, Vec X, Vec F, Params *p)
105+
{
106+
PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
107+
PetscReal two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
108+
PetscScalar u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
109+
Vec localX;
110+
PetscErrorCode ierr;
111+
112+
PetscFunctionBegin;
113+
114+
ierr = DMDAGetInfo(da,PETSC_IGNORE,
115+
&Mx,&My,&Mz,
116+
PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
117+
PETSC_IGNORE,PETSC_IGNORE,
118+
PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
119+
PETSC_IGNORE);
120+
121+
lambda = p->lambda_;
122+
hx = 1.0/(PetscReal)(Mx-1);
123+
hy = 1.0/(PetscReal)(My-1);
124+
hz = 1.0/(PetscReal)(Mz-1);
125+
sc = hx*hy*hz*lambda;
126+
hxhzdhy = hx*hz/hy;
127+
hyhzdhx = hy*hz/hx;
128+
hxhydhz = hx*hy/hz;
129+
130+
/*
131+
132+
*/
133+
ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);
134+
135+
/*
136+
Scatter ghost points to local vector,using the 2-step process
137+
DAGlobalToLocalBegin(),DAGlobalToLocalEnd(). By placing code
138+
between these two statements, computations can be done while
139+
messages are in transition.
140+
*/
141+
ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
142+
ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
143+
144+
/*
145+
Get pointers to vector data.
146+
*/
147+
ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr);
148+
ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
149+
150+
/*
151+
Get local grid boundaries.
152+
*/
153+
ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
154+
155+
/*
156+
Compute function over the locally owned part of the grid.
157+
*/
158+
for (k=zs; k<zs+zm; k++) {
159+
for (j=ys; j<ys+ym; j++) {
160+
for (i=xs; i<xs+xm; i++) {
161+
if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
162+
/* boundary points */
163+
f[k][j][i] = x[k][j][i] - 0.0;
164+
} else {
165+
/* interior grid points */
166+
u = x[k][j][i];
167+
u_east = x[k][j][i+1];
168+
u_west = x[k][j][i-1];
169+
u_north = x[k][j+1][i];
170+
u_south = x[k][j-1][i];
171+
u_up = x[k+1][j][i];
172+
u_down = x[k-1][j][i];
173+
u_xx = (-u_east + two*u - u_west)*hyhzdhx;
174+
u_yy = (-u_north + two*u - u_south)*hxhzdhy;
175+
u_zz = (-u_up + two*u - u_down)*hxhydhz;
176+
f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
177+
}
178+
}
179+
}
180+
}
181+
182+
/*
183+
Restore vectors.
184+
*/
185+
ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
186+
ierr = DMDAVecRestoreArray(da,localX,&x);CHKERRQ(ierr);
187+
ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr);
188+
ierr = PetscLogFlops(11.0*ym*xm);CHKERRQ(ierr);
189+
190+
PetscFunctionReturn(0);
191+
}
192+
193+
194+
#undef __FUNCT__
195+
#define __FUNCT__ "FormJacobian"
196+
PetscErrorCode FormJacobian(DM da, Vec X, Mat J, Params *p)
197+
{
198+
PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
199+
PetscReal lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
200+
PetscScalar v[7],***x;
201+
MatStencil col[7],row;
202+
Vec localX;
203+
PetscErrorCode ierr;
204+
205+
PetscFunctionBegin;
206+
207+
ierr = DMDAGetInfo(da,PETSC_IGNORE,
208+
&Mx,&My,&Mz,
209+
PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
210+
PETSC_IGNORE,PETSC_IGNORE,
211+
PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
212+
PETSC_IGNORE);
213+
214+
lambda = p->lambda_;
215+
hx = 1.0/(PetscReal)(Mx-1);
216+
hy = 1.0/(PetscReal)(My-1);
217+
hz = 1.0/(PetscReal)(Mz-1);
218+
sc = hx*hy*hz*lambda;
219+
hxhzdhy = hx*hz/hy;
220+
hyhzdhx = hy*hz/hx;
221+
hxhydhz = hx*hy/hz;
222+
223+
/*
224+
225+
*/
226+
ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr);
227+
228+
/*
229+
Scatter ghost points to local vector, using the 2-step process
230+
DAGlobalToLocalBegin(), DAGlobalToLocalEnd(). By placing code
231+
between these two statements, computations can be done while
232+
messages are in transition.
233+
*/
234+
ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
235+
ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
236+
237+
/*
238+
Get pointer to vector data.
239+
*/
240+
ierr = DMDAVecGetArray(da,localX,&x);CHKERRQ(ierr);
241+
242+
/*
243+
Get local grid boundaries.
244+
*/
245+
ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
246+
247+
/*
248+
Compute entries for the locally owned part of the Jacobian.
249+
250+
- Currently, all PETSc parallel matrix formats are partitioned by
251+
contiguous chunks of rows across the processors.
252+
253+
- Each processor needs to insert only elements that it owns
254+
locally (but any non-local elements will be sent to the
255+
appropriate processor during matrix assembly).
256+
257+
- Here, we set all entries for a particular row at once.
258+
259+
- We can set matrix entries either using either
260+
MatSetValuesLocal() or MatSetValues(), as discussed above.
261+
*/
262+
for (k=zs; k<zs+zm; k++) {
263+
for (j=ys; j<ys+ym; j++) {
264+
for (i=xs; i<xs+xm; i++) {
265+
row.k = k; row.j = j; row.i = i;
266+
/* boundary points */
267+
if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
268+
v[0] = 1.0;
269+
ierr = MatSetValuesStencil(J,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
270+
} else {
271+
/* interior grid points */
272+
v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j; col[0].i = i;
273+
v[1] = -hxhzdhy; col[1].k=k; col[1].j=j-1;col[1].i = i;
274+
v[2] = -hyhzdhx; col[2].k=k; col[2].j=j; col[2].i = i-1;
275+
v[3] = 2.0*(hyhzdhx+hxhzdhy+hxhydhz)-sc*PetscExpScalar(x[k][j][i]);col[3].k=row.k;col[3].j=row.j;col[3].i = row.i;
276+
v[4] = -hyhzdhx; col[4].k=k; col[4].j=j; col[4].i = i+1;
277+
v[5] = -hxhzdhy; col[5].k=k; col[5].j=j+1;col[5].i = i;
278+
v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j; col[6].i = i;
279+
ierr = MatSetValuesStencil(J,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
280+
}
281+
}
282+
}
283+
}
284+
ierr = DMDAVecRestoreArray(da,localX,&x);CHKERRQ(ierr);
285+
ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr);
286+
287+
/*
288+
Assemble matrix, using the 2-step process: MatAssemblyBegin(),
289+
MatAssemblyEnd().
290+
*/
291+
ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
292+
ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
293+
294+
/*
295+
Tell the matrix we will never add a new nonzero location to the
296+
matrix. If we do, it will generate an error.
297+
*/
298+
ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
299+
300+
PetscFunctionReturn(0);
301+
}

‎petsc4py-tutorial/Bratu3Dimpl.h‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
#ifndef BRATU3D_H
2+
#define BRATU3D_H
3+
4+
#include <petsc.h>
5+
6+
#if PETSC_VERSION_(3,1,0)
7+
#include <petscvec.h>
8+
#include <petscmat.h>
9+
#include <petscda.h>
10+
#endif
11+
12+
typedef struct Params {
13+
double lambda_;
14+
} Params;
15+
16+
PetscErrorCode FormInitGuess(DM da, Vec x, Params *p);
17+
PetscErrorCode FormFunction(DM da, Vec x, Vec F, Params *p);
18+
PetscErrorCode FormJacobian(DM da, Vec x, Mat J, Params *p);
19+
20+
#endif /* !BRATU3D_H */

‎petsc4py-tutorial/CavityFlow2D.pyx‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
from petsc4py.PETSc cimport Vec, PetscVec
2+
from petsc4py.PETSc cimport DA, PetscDM
3+
from petsc4py.PETSc cimport SNES, PetscSNES
4+
5+
from petsc4py.PETSc import Error
6+
7+
cdef extern from "CavityFlow2Dimpl.h":
8+
ctypedef struct Params:
9+
double alpha_, lidvelocity_, prandtl_, grashof_
10+
int FormInitGuess(PetscDM da, PetscVec x, Params *p)
11+
int FormFunction (PetscDM da, PetscVec x, PetscVec F, Params *p)
12+
13+
def formInitGuess(Vec x, DA da, double alpha_, lidvelocity_, prandtl_, grashof_):
14+
cdef int ierr
15+
cdef Params p = {"alpha_" : alpha_, "lidvelocity_" : lidvelocity_,
16+
"prandtl_" : prandtl_, "grashof_" : grashof_}
17+
ierr = FormInitGuess(da.dm, x.vec, &p)
18+
if ierr != 0: raise Error(ierr)
19+
20+
def formFunction(SNES snes, Vec x, Vec f, DA da, double alpha_, lidvelocity_, prandtl_, grashof_):
21+
cdef int ierr
22+
cdef Params p = {"alpha_" : alpha_, "lidvelocity_" : lidvelocity_,
23+
"prandtl_" : prandtl_, "grashof_" : grashof_}
24+
25+
ierr = FormFunction(da.dm, x.vec, f.vec, &p)
26+
if ierr != 0: raise Error(ierr)

0 commit comments

Comments
 (0)