Actual source code: ex29.c
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 2d
4: Processors: n
5: T*/
7: /*
8: Added at the request of Marc Garbey.
10: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
12: div \rho grad u = f, 0 < x,y < 1,
14: with forcing function
16: f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}
18: with Dirichlet boundary conditions
20: u = f(x,y) for x = 0, x = 1, y = 0, y = 1
22: or pure Neumman boundary conditions
24: This uses multigrid to solve the linear system
25: */
27: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
29: #include petscda.h
30: #include petscksp.h
31: #include petscmg.h
32: #include petscdmmg.h
38: typedef enum {DIRICHLET, NEUMANN} BCType;
40: typedef struct {
41: PetscScalar rho;
42: PetscScalar nu;
43: BCType bcType;
44: } UserContext;
48: int main(int argc,char **argv)
49: {
50: DMMG *dmmg;
51: DA da;
52: UserContext user;
53: PetscReal norm;
54: const char *bcTypes[2] = {"dirichlet","neumann"};
56: PetscInt l,bc;
58: PetscInitialize(&argc,&argv,(char *)0,help);
60: DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
61: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
62: DMMGSetDM(dmmg,(DM)da);
63: DADestroy(da);
64: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
65: DMMGSetUser(dmmg,l,&user);
66: }
68: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMMG");
69: user.rho = 1.0;
70: PetscOptionsScalar("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, PETSC_NULL);
71: user.nu = 0.1;
72: PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, PETSC_NULL);
73: bc = (PetscInt)DIRICHLET;
74: PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
75: user.bcType = (BCType)bc;
76: PetscOptionsEnd();
78: DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
79: if (user.bcType == NEUMANN) {
80: DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
81: }
83: DMMGSolve(dmmg);
85: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
86: VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
87: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
88: /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm); */
89: VecAssemblyBegin(DMMGGetx(dmmg));
90: VecAssemblyEnd(DMMGGetx(dmmg));
91: VecView_VTK(DMMGGetRHS(dmmg), "rhs.vtk", bcTypes[user.bcType]);
92: VecView_VTK(DMMGGetx(dmmg), "solution.vtk", bcTypes[user.bcType]);
94: DMMGDestroy(dmmg);
95: PetscFinalize();
97: return 0;
98: }
102: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
103: {
104: DA da = (DA)dmmg->dm;
105: UserContext *user = (UserContext *) dmmg->user;
107: PetscInt i,j,mx,my,xm,ym,xs,ys;
108: PetscScalar Hx,Hy;
109: PetscScalar **array;
112: DAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0);
113: Hx = 1.0 / (PetscReal)(mx-1);
114: Hy = 1.0 / (PetscReal)(my-1);
115: DAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
116: DAVecGetArray(da, b, &array);
117: for (j=ys; j<ys+ym; j++){
118: for(i=xs; i<xs+xm; i++){
119: array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
120: }
121: }
122: DAVecRestoreArray(da, b, &array);
123: VecAssemblyBegin(b);
124: VecAssemblyEnd(b);
126: /* force right hand side to be consistent for singular matrix */
127: /* note this is really a hack, normally the model would provide you with a consistent right handside */
128: if (user->bcType == NEUMANN) {
129: MatNullSpace nullspace;
131: KSPGetNullSpace(dmmg->ksp,&nullspace);
132: MatNullSpaceRemove(nullspace,b,PETSC_NULL);
133: }
134: return(0);
135: }
137:
140: PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscScalar centerRho, PetscScalar *rho)
141: {
143: if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
144: *rho = centerRho;
145: } else {
146: *rho = 1.0;
147: }
148: return(0);
149: }
153: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J,Mat jac)
154: {
155: DA da = (DA) dmmg->dm;
156: UserContext *user = (UserContext *) dmmg->user;
157: PetscScalar centerRho = user->rho;
159: PetscInt i,j,mx,my,xm,ym,xs,ys,num;
160: PetscScalar v[5],Hx,Hy,HydHx,HxdHy,rho;
161: MatStencil row, col[5];
164: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
165: Hx = 1.0 / (PetscReal)(mx-1);
166: Hy = 1.0 / (PetscReal)(my-1);
167: HxdHy = Hx/Hy;
168: HydHx = Hy/Hx;
169: DAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
170: for (j=ys; j<ys+ym; j++){
171: for(i=xs; i<xs+xm; i++){
172: row.i = i; row.j = j;
173: ComputeRho(i, j, mx, my, centerRho, &rho);
174: if (i==0 || j==0 || i==mx-1 || j==my-1) {
175: if (user->bcType == DIRICHLET) {
176: v[0] = 2.0*rho*(HxdHy + HydHx);
177: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
178: } else if (user->bcType == NEUMANN) {
179: num = 0;
180: if (j!=0) {
181: v[num] = -rho*HxdHy; col[num].i = i; col[num].j = j-1;
182: num++;
183: }
184: if (i!=0) {
185: v[num] = -rho*HydHx; col[num].i = i-1; col[num].j = j;
186: num++;
187: }
188: if (i!=mx-1) {
189: v[num] = -rho*HydHx; col[num].i = i+1; col[num].j = j;
190: num++;
191: }
192: if (j!=my-1) {
193: v[num] = -rho*HxdHy; col[num].i = i; col[num].j = j+1;
194: num++;
195: }
196: v[num] = (num/2.0)*rho*(HxdHy + HydHx); col[num].i = i; col[num].j = j;
197: num++;
198: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
199: }
200: } else {
201: v[0] = -rho*HxdHy; col[0].i = i; col[0].j = j-1;
202: v[1] = -rho*HydHx; col[1].i = i-1; col[1].j = j;
203: v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
204: v[3] = -rho*HydHx; col[3].i = i+1; col[3].j = j;
205: v[4] = -rho*HxdHy; col[4].i = i; col[4].j = j+1;
206: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
207: }
208: }
209: }
210: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
211: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
212: return(0);
213: }
217: PetscErrorCode VecView_VTK(Vec x, const char filename[], const char bcName[])
218: {
219: MPI_Comm comm;
220: DA da;
221: Vec coords;
222: PetscViewer viewer;
223: PetscScalar *array, *values;
224: PetscInt n, N, maxn, mx, my, dof;
225: PetscInt i, p;
226: MPI_Status status;
227: PetscMPIInt rank, size, tag, nn;
228: PetscErrorCode ierr;
231: PetscObjectGetComm((PetscObject) x, &comm);
232: PetscViewerASCIIOpen(comm, filename, &viewer);
234: VecGetSize(x, &N);
235: VecGetLocalSize(x, &n);
236: PetscObjectQuery((PetscObject) x, "DA", (PetscObject *) &da);
237: if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DA");
239: DAGetInfo(da, 0, &mx, &my, 0,0,0,0, &dof,0,0,0);
241: PetscViewerASCIIPrintf(viewer, "# vtk DataFile Version 2.0\n");
242: PetscViewerASCIIPrintf(viewer, "Inhomogeneous Poisson Equation with %s boundary conditions\n", bcName);
243: PetscViewerASCIIPrintf(viewer, "ASCII\n");
244: /* get coordinates of nodes */
245: DAGetCoordinates(da, &coords);
246: if (!coords) {
247: DASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
248: DAGetCoordinates(da, &coords);
249: }
250: PetscViewerASCIIPrintf(viewer, "DATASET RECTILINEAR_GRID\n");
251: PetscViewerASCIIPrintf(viewer, "DIMENSIONS %d %d %d\n", mx, my, 1);
252: VecGetArray(coords, &array);
253: PetscViewerASCIIPrintf(viewer, "X_COORDINATES %d double\n", mx);
254: for(i = 0; i < mx; i++) {
255: PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*2]));
256: }
257: PetscViewerASCIIPrintf(viewer, "\n");
258: PetscViewerASCIIPrintf(viewer, "Y_COORDINATES %d double\n", my);
259: for(i = 0; i < my; i++) {
260: PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*mx*2+1]));
261: }
262: PetscViewerASCIIPrintf(viewer, "\n");
263: PetscViewerASCIIPrintf(viewer, "Z_COORDINATES %d double\n", 1);
264: PetscViewerASCIIPrintf(viewer, "%G\n", 0.0);
265: VecRestoreArray(coords, &array);
266: VecDestroy(coords);
267: PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", N);
268: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", dof);
269: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
270: VecGetArray(x, &array);
271: /* Determine maximum message to arrive */
272: MPI_Comm_rank(comm, &rank);
273: MPI_Comm_size(comm, &size);
274: MPI_Reduce(&n, &maxn, 1, MPIU_INT, MPI_MAX, 0, comm);
275: tag = ((PetscObject) viewer)->tag;
276: if (!rank) {
277: PetscMalloc((maxn+1) * sizeof(PetscScalar), &values);
278: for(i = 0; i < n; i++) {
279: PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
280: }
281: for(p = 1; p < size; p++) {
282: MPI_Recv(values, (PetscMPIInt) n, MPIU_SCALAR, p, tag, comm, &status);
283: MPI_Get_count(&status, MPIU_SCALAR, &nn);
284: for(i = 0; i < nn; i++) {
285: PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
286: }
287: }
288: PetscFree(values);
289: } else {
290: MPI_Send(array, n, MPIU_SCALAR, 0, tag, comm);
291: }
292: VecRestoreArray(x, &array);
293: PetscViewerFlush(viewer);
294: PetscViewerDestroy(viewer);
295: return(0);
296: }