Actual source code: ex20.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: #include <petscksp.h>

  8: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
  9: {
 10:   Ke[0]  = H / 6.0;
 11:   Ke[1]  = -.125 * H;
 12:   Ke[2]  = H / 12.0;
 13:   Ke[3]  = -.125 * H;
 14:   Ke[4]  = -.125 * H;
 15:   Ke[5]  = H / 6.0;
 16:   Ke[6]  = -.125 * H;
 17:   Ke[7]  = H / 12.0;
 18:   Ke[8]  = H / 12.0;
 19:   Ke[9]  = -.125 * H;
 20:   Ke[10] = H / 6.0;
 21:   Ke[11] = -.125 * H;
 22:   Ke[12] = -.125 * H;
 23:   Ke[13] = H / 12.0;
 24:   Ke[14] = -.125 * H;
 25:   Ke[15] = H / 6.0;
 26:   return PETSC_SUCCESS;
 27: }

 29: int main(int argc, char **args)
 30: {
 31:   Mat          C;
 32:   PetscMPIInt  rank, size;
 33:   PetscInt     i, m = 5, N, start, end, M;
 34:   PetscInt     idx[4];
 35:   PetscScalar  Ke[16];
 36:   PetscReal    h;
 37:   Vec          u, b;
 38:   KSP          ksp;
 39:   MatNullSpace nullsp;

 41:   PetscFunctionBeginUser;
 42:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 43:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 44:   N = (m + 1) * (m + 1); /* dimension of matrix */
 45:   M = m * m;             /* number of elements */
 46:   h = 1.0 / m;           /* mesh width */
 47:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 48:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 50:   /* Create stiffness matrix */
 51:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 52:   PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
 53:   PetscCall(MatSetFromOptions(C));
 54:   PetscCall(MatSetUp(C));
 55:   start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
 56:   end   = start + M / size + ((M % size) > rank);

 58:   /* Assemble matrix */
 59:   PetscCall(FormElementStiffness(h * h, Ke)); /* element stiffness for Laplacian */
 60:   for (i = start; i < end; i++) {
 61:     /* location of lower left corner of element */
 62:     /* node numbers for the four corners of element */
 63:     idx[0] = (m + 1) * (i / m) + (i % m);
 64:     idx[1] = idx[0] + 1;
 65:     idx[2] = idx[1] + m + 1;
 66:     idx[3] = idx[2] - 1;
 67:     PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES));
 68:   }
 69:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 70:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 72:   /* Create right-hand side and solution vectors */
 73:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
 74:   PetscCall(VecSetSizes(u, PETSC_DECIDE, N));
 75:   PetscCall(VecSetFromOptions(u));
 76:   PetscCall(PetscObjectSetName((PetscObject)u, "Approx. Solution"));
 77:   PetscCall(VecDuplicate(u, &b));
 78:   PetscCall(PetscObjectSetName((PetscObject)b, "Right hand side"));

 80:   PetscCall(VecSet(b, 1.0));
 81:   PetscCall(VecSetValue(b, 0, 1.2, ADD_VALUES));
 82:   PetscCall(VecSet(u, 0.0));

 84:   /* Solve linear system */
 85:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 86:   PetscCall(KSPSetOperators(ksp, C, C));
 87:   PetscCall(KSPSetFromOptions(ksp));
 88:   PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));

 90:   PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nullsp));
 91:   /*
 92:      The KSP solver will remove this nullspace from the solution at each iteration
 93:   */
 94:   PetscCall(MatSetNullSpace(C, nullsp));
 95:   /*
 96:      The KSP solver will remove from the right-hand side any portion in this nullspace, thus making the linear system consistent.
 97:   */
 98:   PetscCall(MatSetTransposeNullSpace(C, nullsp));
 99:   PetscCall(MatNullSpaceDestroy(&nullsp));

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

103:   /* Free work space */
104:   PetscCall(KSPDestroy(&ksp));
105:   PetscCall(VecDestroy(&u));
106:   PetscCall(VecDestroy(&b));
107:   PetscCall(MatDestroy(&C));
108:   PetscCall(PetscFinalize());
109:   return 0;
110: }

112: /*TEST

114:     test:
115:       args: -ksp_monitor_short

117: TEST*/