Actual source code: ex3.c


  2: static char help[] = "Bilinear elements on the unit square for Laplacian.  To test the parallel\n\
  3: matrix assembly, the matrix is intentionally laid out across processors\n\
  4: differently from the way it is assembled.  Input arguments are:\n\
  5:   -m <size> : problem size\n\n";

  7: /*
  8:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  9:   automatically includes:
 10:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 11:      petscmat.h - matrices
 12:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 13:      petscviewer.h - viewers               petscpc.h  - preconditioners
 14: */
 15: #include <petscksp.h>

 17: /* Declare user-defined routines */
 18: extern PetscErrorCode FormElementStiffness(PetscReal, PetscScalar *);
 19: extern PetscErrorCode FormElementRhs(PetscScalar, PetscScalar, PetscReal, PetscScalar *);

 21: int main(int argc, char **args)
 22: {
 23:   Vec         u, b, ustar; /* approx solution, RHS, exact solution */
 24:   Mat         A;           /* linear system matrix */
 25:   KSP         ksp;         /* Krylov subspace method context */
 26:   PetscInt    N;           /* dimension of system (global) */
 27:   PetscInt    M;           /* number of elements (global) */
 28:   PetscMPIInt rank;        /* processor rank */
 29:   PetscMPIInt size;        /* size of communicator */
 30:   PetscScalar Ke[16];      /* element matrix */
 31:   PetscScalar r[4];        /* element vector */
 32:   PetscReal   h;           /* mesh width */
 33:   PetscReal   norm;        /* norm of solution error */
 34:   PetscScalar x, y;
 35:   PetscInt    idx[4], count, *rows, i, m = 5, start, end, its;

 38:   PetscInitialize(&argc, &args, (char *)0, help);
 39:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 40:   N = (m + 1) * (m + 1);
 41:   M = m * m;
 42:   h = 1.0 / m;
 43:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 44:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 47:          Compute the matrix and right-hand-side vector that define
 48:          the linear system, Au = b.
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 51:   /*
 52:      Create stiffness matrix
 53:   */
 54:   MatCreate(PETSC_COMM_WORLD, &A);
 55:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N);
 56:   MatSetFromOptions(A);
 57:   MatSeqAIJSetPreallocation(A, 9, NULL);
 58:   MatMPIAIJSetPreallocation(A, 9, NULL, 8, NULL);
 59:   start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
 60:   end   = start + M / size + ((M % size) > rank);

 62:   /*
 63:      Assemble matrix
 64:   */
 65:   FormElementStiffness(h * h, Ke);
 66:   for (i = start; i < end; i++) {
 67:     /* node numbers for the four corners of element */
 68:     idx[0] = (m + 1) * (i / m) + (i % m);
 69:     idx[1] = idx[0] + 1;
 70:     idx[2] = idx[1] + m + 1;
 71:     idx[3] = idx[2] - 1;
 72:     MatSetValues(A, 4, idx, 4, idx, Ke, ADD_VALUES);
 73:   }
 74:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 75:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 77:   /*
 78:      Create right-hand-side and solution vectors
 79:   */
 80:   VecCreate(PETSC_COMM_WORLD, &u);
 81:   VecSetSizes(u, PETSC_DECIDE, N);
 82:   VecSetFromOptions(u);
 83:   PetscObjectSetName((PetscObject)u, "Approx. Solution");
 84:   VecDuplicate(u, &b);
 85:   PetscObjectSetName((PetscObject)b, "Right hand side");
 86:   VecDuplicate(b, &ustar);
 87:   VecSet(u, 0.0);
 88:   VecSet(b, 0.0);

 90:   /*
 91:      Assemble right-hand-side vector
 92:   */
 93:   for (i = start; i < end; i++) {
 94:     /* location of lower left corner of element */
 95:     x = h * (i % m);
 96:     y = h * (i / m);
 97:     /* node numbers for the four corners of element */
 98:     idx[0] = (m + 1) * (i / m) + (i % m);
 99:     idx[1] = idx[0] + 1;
100:     idx[2] = idx[1] + m + 1;
101:     idx[3] = idx[2] - 1;
102:     FormElementRhs(x, y, h * h, r);
103:     VecSetValues(b, 4, idx, r, ADD_VALUES);
104:   }
105:   VecAssemblyBegin(b);
106:   VecAssemblyEnd(b);

108:   /*
109:      Modify matrix and right-hand-side for Dirichlet boundary conditions
110:   */
111:   PetscMalloc1(4 * m, &rows);
112:   for (i = 0; i < m + 1; i++) {
113:     rows[i]             = i;               /* bottom */
114:     rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
115:   }
116:   count = m + 1; /* left side */
117:   for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
118:   count = 2 * m; /* left side */
119:   for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
120:   for (i = 0; i < 4 * m; i++) {
121:     y = h * (rows[i] / (m + 1));
122:     VecSetValues(u, 1, &rows[i], &y, INSERT_VALUES);
123:     VecSetValues(b, 1, &rows[i], &y, INSERT_VALUES);
124:   }
125:   MatZeroRows(A, 4 * m, rows, 1.0, 0, 0);
126:   PetscFree(rows);

128:   VecAssemblyBegin(u);
129:   VecAssemblyEnd(u);
130:   VecAssemblyBegin(b);
131:   VecAssemblyEnd(b);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:                 Create the linear solver and set various options
135:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

137:   KSPCreate(PETSC_COMM_WORLD, &ksp);
138:   KSPSetOperators(ksp, A, A);
139:   KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
140:   KSPSetFromOptions(ksp);

142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:                       Solve the linear system
144:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

146:   KSPSolve(ksp, b, u);

148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149:                       Check solution and clean up
150:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

152:   /* Check error */
153:   VecGetOwnershipRange(ustar, &start, &end);
154:   for (i = start; i < end; i++) {
155:     y = h * (i / (m + 1));
156:     VecSetValues(ustar, 1, &i, &y, INSERT_VALUES);
157:   }
158:   VecAssemblyBegin(ustar);
159:   VecAssemblyEnd(ustar);
160:   VecAXPY(u, -1.0, ustar);
161:   VecNorm(u, NORM_2, &norm);
162:   KSPGetIterationNumber(ksp, &its);
163:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its);

165:   /*
166:      Free work space.  All PETSc objects should be destroyed when they
167:      are no longer needed.
168:   */
169:   KSPDestroy(&ksp);
170:   VecDestroy(&u);
171:   VecDestroy(&ustar);
172:   VecDestroy(&b);
173:   MatDestroy(&A);

175:   /*
176:      Always call PetscFinalize() before exiting a program.  This routine
177:        - finalizes the PETSc libraries as well as MPI
178:        - provides summary and diagnostic information if certain runtime
179:          options are chosen (e.g., -log_view).
180:   */
181:   PetscFinalize();
182:   return 0;
183: }

185: /* --------------------------------------------------------------------- */
186: /* element stiffness for Laplacian */
187: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
188: {
190:   Ke[0]  = H / 6.0;
191:   Ke[1]  = -.125 * H;
192:   Ke[2]  = H / 12.0;
193:   Ke[3]  = -.125 * H;
194:   Ke[4]  = -.125 * H;
195:   Ke[5]  = H / 6.0;
196:   Ke[6]  = -.125 * H;
197:   Ke[7]  = H / 12.0;
198:   Ke[8]  = H / 12.0;
199:   Ke[9]  = -.125 * H;
200:   Ke[10] = H / 6.0;
201:   Ke[11] = -.125 * H;
202:   Ke[12] = -.125 * H;
203:   Ke[13] = H / 12.0;
204:   Ke[14] = -.125 * H;
205:   Ke[15] = H / 6.0;
206:   return 0;
207: }
208: /* --------------------------------------------------------------------- */
209: PetscErrorCode FormElementRhs(PetscScalar x, PetscScalar y, PetscReal H, PetscScalar *r)
210: {
212:   r[0] = 0.;
213:   r[1] = 0.;
214:   r[2] = 0.;
215:   r[3] = 0.0;
216:   return 0;
217: }

219: /*TEST

221:    test:
222:       suffix: 1
223:       nsize: 2
224:       args: -ksp_monitor_short

226: TEST*/