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*/