Actual source code: ex3.c

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

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

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

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

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

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

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

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

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

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

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

    127:   PetscCall(VecAssemblyBegin(u));
    128:   PetscCall(VecAssemblyEnd(u));
    129:   PetscCall(VecAssemblyBegin(b));
    130:   PetscCall(VecAssemblyEnd(b));

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

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

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

    145:   PetscCall(KSPSolve(ksp, b, u));

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

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

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

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

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

    218: /*TEST

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

    225: TEST*/