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^{-x^2/\nu} e^{-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 <petscdmda.h>
30: #include <petscksp.h>
31: #include <petscpcmg.h>
33: extern PetscErrorCode ComputeMatrix(DM,Vec,Mat,Mat,MatStructure*);
34: extern PetscErrorCode ComputeRHS(DM,Vec,Vec);
35: extern PetscErrorCode VecView_VTK(Vec, const char [], const char []);
37: typedef enum {DIRICHLET, NEUMANN} BCType;
39: typedef struct {
40: PetscScalar rho;
41: PetscScalar nu;
42: BCType bcType;
43: } UserContext;
47: int main(int argc,char **argv)
48: {
49: KSP ksp;
50: DM da;
51: UserContext user;
52: const char *bcTypes[2] = {"dirichlet","neumann"};
54: PetscInt bc;
55: Vec b,x;
57: PetscInitialize(&argc,&argv,(char *)0,help);
59: KSPCreate(PETSC_COMM_WORLD,&ksp);
60: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
61: DMSetApplicationContext(da,&user);
63: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMqq");
64: user.rho = 1.0;
65: PetscOptionsScalar("-rho", "The conductivity", "ex29.c", user.rho, &user.rho, PETSC_NULL);
66: user.nu = 0.1;
67: PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", user.nu, &user.nu, PETSC_NULL);
68: bc = (PetscInt)DIRICHLET;
69: PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
70: user.bcType = (BCType)bc;
71: PetscOptionsEnd();
73: DMSetFunction(da,ComputeRHS);
74: DMSetJacobian(da,ComputeMatrix);
75: KSPSetDM(ksp,da);
76: KSPSetFromOptions(ksp);
77: KSPSetUp(ksp);
78: KSPSolve(ksp,PETSC_NULL,PETSC_NULL);
79: KSPGetSolution(ksp,&x);
80: KSPGetRhs(ksp,&b);
81: VecView_VTK(b, "rhs.vtk", bcTypes[user.bcType]);
82: VecView_VTK(x, "solution.vtk", bcTypes[user.bcType]);
84: DMDestroy(&da);
85: KSPDestroy(&ksp);
86: PetscFinalize();
88: return 0;
89: }
93: PetscErrorCode ComputeRHS(DM da,Vec x, Vec b)
94: {
95: UserContext *user;
97: PetscInt i,j,mx,my,xm,ym,xs,ys;
98: PetscScalar Hx,Hy;
99: PetscScalar **array;
102: DMGetApplicationContext(da,&user);
103: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);
104: Hx = 1.0 / (PetscReal)(mx-1);
105: Hy = 1.0 / (PetscReal)(my-1);
106: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
107: DMDAVecGetArray(da, b, &array);
108: for (j=ys; j<ys+ym; j++){
109: for(i=xs; i<xs+xm; i++){
110: array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
111: }
112: }
113: DMDAVecRestoreArray(da, b, &array);
114: VecAssemblyBegin(b);
115: VecAssemblyEnd(b);
117: /* force right hand side to be consistent for singular matrix */
118: /* note this is really a hack, normally the model would provide you with a consistent right handside */
119: if (user->bcType == NEUMANN) {
120: MatNullSpace nullspace;
122: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
123: MatNullSpaceRemove(nullspace,b,PETSC_NULL);
124: MatNullSpaceDestroy(&nullspace);
125: }
126: return(0);
127: }
129:
132: PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscScalar centerRho, PetscScalar *rho)
133: {
135: if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
136: *rho = centerRho;
137: } else {
138: *rho = 1.0;
139: }
140: return(0);
141: }
145: PetscErrorCode ComputeMatrix(DM da, Vec x,Mat J,Mat jac,MatStructure *str)
146: {
147: UserContext *user;
148: PetscScalar centerRho;
150: PetscInt i,j,mx,my,xm,ym,xs,ys,num;
151: PetscScalar v[5],Hx,Hy,HydHx,HxdHy,rho;
152: MatStencil row, col[5];
155: DMGetApplicationContext(da,&user);
156: centerRho = user->rho;
157: DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);
158: Hx = 1.0 / (PetscReal)(mx-1);
159: Hy = 1.0 / (PetscReal)(my-1);
160: HxdHy = Hx/Hy;
161: HydHx = Hy/Hx;
162: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
163: for (j=ys; j<ys+ym; j++){
164: for(i=xs; i<xs+xm; i++){
165: row.i = i; row.j = j;
166: ComputeRho(i, j, mx, my, centerRho, &rho);
167: if (i==0 || j==0 || i==mx-1 || j==my-1) {
168: if (user->bcType == DIRICHLET) {
169: v[0] = 2.0*rho*(HxdHy + HydHx);
170: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
171: } else if (user->bcType == NEUMANN) {
172: num = 0;
173: if (j!=0) {
174: v[num] = -rho*HxdHy; col[num].i = i; col[num].j = j-1;
175: num++;
176: }
177: if (i!=0) {
178: v[num] = -rho*HydHx; col[num].i = i-1; col[num].j = j;
179: num++;
180: }
181: if (i!=mx-1) {
182: v[num] = -rho*HydHx; col[num].i = i+1; col[num].j = j;
183: num++;
184: }
185: if (j!=my-1) {
186: v[num] = -rho*HxdHy; col[num].i = i; col[num].j = j+1;
187: num++;
188: }
189: v[num] = (num/2.0)*rho*(HxdHy + HydHx); col[num].i = i; col[num].j = j;
190: num++;
191: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
192: }
193: } else {
194: v[0] = -rho*HxdHy; col[0].i = i; col[0].j = j-1;
195: v[1] = -rho*HydHx; col[1].i = i-1; col[1].j = j;
196: v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
197: v[3] = -rho*HydHx; col[3].i = i+1; col[3].j = j;
198: v[4] = -rho*HxdHy; col[4].i = i; col[4].j = j+1;
199: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
200: }
201: }
202: }
203: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
204: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
205: if (user->bcType == NEUMANN) {
206: MatNullSpace nullspace;
208: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
209: MatSetNullSpace(jac,nullspace);
210: MatNullSpaceDestroy(&nullspace);
211: }
213: return(0);
214: }
218: PetscErrorCode VecView_VTK(Vec x, const char filename[], const char bcName[])
219: {
220: MPI_Comm comm;
221: DM da;
222: Vec coords;
223: PetscViewer viewer;
224: PetscScalar *array, *values;
225: PetscInt n, N, maxn, mx, my, dof;
226: PetscInt i, p;
227: MPI_Status status;
228: PetscMPIInt rank, size, tag, nn;
229: PetscErrorCode ierr;
230: PetscBool flg;
233: PetscObjectGetComm((PetscObject) x, &comm);
234: PetscViewerASCIIOpen(comm, filename, &viewer);
236: VecGetSize(x, &N);
237: VecGetLocalSize(x, &n);
238: PetscObjectQuery((PetscObject) x, "DM", (PetscObject *) &da);
239: if (!da) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
240: PetscTypeCompare((PetscObject)da,DMDA,&flg);
241: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
243: DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0, &dof,0,0,0,0,0);
245: PetscViewerASCIIPrintf(viewer, "# vtk DataFile Version 2.0\n");
246: PetscViewerASCIIPrintf(viewer, "Inhomogeneous Poisson Equation with %s boundary conditions\n", bcName);
247: PetscViewerASCIIPrintf(viewer, "ASCII\n");
248: /* get coordinates of nodes */
249: DMDAGetCoordinates(da, &coords);
250: if (!coords) {
251: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
252: DMDAGetCoordinates(da, &coords);
253: }
254: PetscViewerASCIIPrintf(viewer, "DATASET RECTILINEAR_GRID\n");
255: PetscViewerASCIIPrintf(viewer, "DIMENSIONS %d %d %d\n", mx, my, 1);
256: VecGetArray(coords, &array);
257: PetscViewerASCIIPrintf(viewer, "X_COORDINATES %d double\n", mx);
258: for(i = 0; i < mx; i++) {
259: PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*2]));
260: }
261: PetscViewerASCIIPrintf(viewer, "\n");
262: PetscViewerASCIIPrintf(viewer, "Y_COORDINATES %d double\n", my);
263: for(i = 0; i < my; i++) {
264: PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*mx*2+1]));
265: }
266: PetscViewerASCIIPrintf(viewer, "\n");
267: PetscViewerASCIIPrintf(viewer, "Z_COORDINATES %d double\n", 1);
268: PetscViewerASCIIPrintf(viewer, "%G\n", 0.0);
269: VecRestoreArray(coords, &array);
270: PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", N);
271: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", dof);
272: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
273: VecGetArray(x, &array);
274: /* Determine maximum message to arrive */
275: MPI_Comm_rank(comm, &rank);
276: MPI_Comm_size(comm, &size);
277: MPI_Reduce(&n, &maxn, 1, MPIU_INT, MPI_MAX, 0, comm);
278: tag = ((PetscObject) viewer)->tag;
279: if (!rank) {
280: PetscMalloc((maxn+1) * sizeof(PetscScalar), &values);
281: for(i = 0; i < n; i++) {
282: PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
283: }
284: for(p = 1; p < size; p++) {
285: MPI_Recv(values, (PetscMPIInt) n, MPIU_SCALAR, p, tag, comm, &status);
286: MPI_Get_count(&status, MPIU_SCALAR, &nn);
287: for(i = 0; i < nn; i++) {
288: PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
289: }
290: }
291: PetscFree(values);
292: } else {
293: MPI_Send(array, n, MPIU_SCALAR, 0, tag, comm);
294: }
295: VecRestoreArray(x, &array);
296: PetscViewerFlush(viewer);
297: PetscViewerDestroy(&viewer);
298: return(0);
299: }