Actual source code: ex3.c
1: static char help[] = "Tests no-net-flux constraints combined with a supplied near-nullspace.\n\n";
3: #include <petsc/private/pcbddcimpl.h>
5: static PetscErrorCode CheckConstraint(Mat C, Vec v)
6: {
7: Vec c, projected;
8: PetscReal error, norm;
10: PetscFunctionBeginUser;
11: PetscCall(MatCreateVecs(C, &projected, &c));
12: PetscCall(MatMult(C, v, c));
13: PetscCall(MatMultTranspose(C, c, projected));
14: PetscCall(VecAXPY(projected, -1.0, v));
15: PetscCall(VecNorm(projected, NORM_2, &error));
16: PetscCall(VecNorm(v, NORM_2, &norm));
17: PetscCheck(error < PETSC_SMALL * norm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Constraint span does not contain the supplied vector: relative error %g", (double)(error / norm));
18: PetscCall(VecDestroy(&projected));
19: PetscCall(VecDestroy(&c));
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: int main(int argc, char **args)
24: {
25: Mat A, B, local;
26: MatNullSpace nsp = NULL, attached;
27: ISLocalToGlobalMapping map, pmap;
28: KSP ksp = NULL;
29: PC pc;
30: PC_BDDC *bddc;
31: Vec modes[2], x, b, exact, v;
32: PetscScalar *values;
33: PetscInt indices[] = {0, 1, 2, 3, 4, 5};
34: PetscInt nlocal, npressure, pressure, start, end, nmodes = 2, nconstraints, expected = 0, offset = 0;
35: PetscMPIInt rank, size, active;
36: PetscBool constant = PETSC_FALSE, empty_rank = PETSC_FALSE, dependent = PETSC_FALSE, has_const = PETSC_FALSE, reset = PETSC_FALSE;
37: PetscReal error;
39: PetscFunctionBeginUser;
40: PetscCall(PetscInitialize(&argc, &args, NULL, help));
41: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
42: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
43: PetscCall(PetscOptionsGetInt(NULL, NULL, "-user_modes", &nmodes, NULL));
44: PetscCall(PetscOptionsGetBool(NULL, NULL, "-constant", &constant, NULL));
45: PetscCall(PetscOptionsGetBool(NULL, NULL, "-empty_rank", &empty_rank, NULL));
46: PetscCall(PetscOptionsGetBool(NULL, NULL, "-reset", &reset, NULL));
47: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dependent", &dependent, NULL));
48: PetscCheck(nmodes >= 0 && nmodes <= 2, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Use zero, one, or two user modes");
49: if (dependent) {
50: nmodes = 1;
51: constant = PETSC_FALSE;
52: }
53: active = size - (empty_rank ? 1 : 0);
54: nlocal = rank < active ? 6 : 0;
55: npressure = rank < active ? 1 : 0;
56: pressure = rank;
57: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, nlocal, indices, PETSC_COPY_VALUES, &map));
58: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, npressure, &pressure, PETSC_COPY_VALUES, &pmap));
59: PetscCall(MatCreateIS(PETSC_COMM_WORLD, 1, PETSC_DECIDE, PETSC_DECIDE, 6, 6, map, map, &A));
60: PetscCall(MatISSetPreallocation(A, 6, NULL, 6, NULL));
61: PetscCall(MatISGetLocalMat(A, &local));
62: for (PetscInt i = 0; i < nlocal; i++)
63: for (PetscInt j = 0; j < nlocal; j++) PetscCall(MatSetValue(local, i, j, i == j ? 2.0 : 0.1, INSERT_VALUES));
64: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
65: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
66: PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));
67: PetscCall(MatCreateIS(PETSC_COMM_WORLD, 1, npressure, PETSC_DECIDE, active, 6, pmap, map, &B));
68: PetscCall(MatISSetPreallocation(B, 6, NULL, 6, NULL));
69: PetscCall(MatISGetLocalMat(B, &local));
70: for (PetscInt i = 0; i < nlocal; i++) PetscCall(MatSetValue(local, 0, i, (rank ? -1.0 : active - 1.0) * (i + 1), INSERT_VALUES));
71: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
74: PetscCall(MatCreateVecs(A, &modes[0], &b));
75: PetscCall(VecDuplicate(modes[0], &modes[1]));
76: PetscCall(VecDuplicate(modes[0], &x));
77: PetscCall(VecDuplicate(modes[0], &exact));
78: PetscCall(VecGetOwnershipRange(modes[0], &start, &end));
79: for (PetscInt k = 0; k < 2; k++) {
80: PetscCall(VecGetArray(modes[k], &values));
81: for (PetscInt i = start; i < end; i++) values[i - start] = (i == 2 * k ? 1.0 : 0.0) - (i == 2 * k + 1 ? 1.0 : 0.0);
82: if (dependent && !k)
83: for (PetscInt i = start; i < end; i++) values[i - start] = i + 1;
84: PetscCall(VecRestoreArray(modes[k], &values));
85: PetscCall(VecNormalize(modes[k], NULL));
86: }
87: PetscCall(VecSet(exact, 1.0));
88: PetscCall(VecAXPY(exact, 0.5, modes[1]));
89: for (PetscInt step = 0; step < 5; step++) {
90: // Rebuild the interface, replace the user modes, remove them, and reset or recreate the solver.
91: if (step == 4 && reset) {
92: PetscCall(PCReset(pc));
93: PetscCall(PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL));
94: }
95: if (step == 0 || (step == 4 && !reset)) {
96: PetscCall(KSPDestroy(&ksp));
97: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
98: PetscCall(KSPSetType(ksp, KSPCG));
99: PetscCall(KSPSetOperators(ksp, A, A));
100: PetscCall(KSPGetPC(ksp, &pc));
101: PetscCall(PCSetType(pc, PCBDDC));
102: PetscCall(PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL));
103: PetscCall(KSPSetFromOptions(ksp));
104: }
105: if (step != 1) {
106: PetscCall(MatNullSpaceDestroy(&nsp));
107: has_const = step == 2 || step == 3 ? PETSC_FALSE : constant;
108: offset = step == 2 ? 1 : 0;
109: expected = step == 2 ? 1 : step == 3 ? 0 : nmodes;
110: if (has_const || expected) PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, has_const, expected, modes + offset, &nsp));
111: PetscCall(MatSetNearNullSpace(A, nsp));
112: }
113: bddc = (PC_BDDC *)pc->data;
114: if (step == 1) bddc->recompute_topography = PETSC_TRUE;
115: PetscCall(MatScale(A, 1.01));
116: PetscCall(MatMult(A, exact, b));
117: PetscCall(KSPSetOperators(ksp, A, A));
118: PetscCall(KSPSolve(ksp, b, x));
119: PetscCall(MatGetNearNullSpace(A, &attached));
120: PetscCheck(attached == nsp, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The supplied near-nullspace was replaced");
121: PetscCall(VecAXPY(x, -1.0, exact));
122: PetscCall(VecNorm(x, NORM_INFINITY, &error));
123: PetscCheck(error < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Solution error %g", (double)error);
125: PetscCall(MatGetSize(bddc->ConstraintMatrix, &nconstraints, NULL));
126: PetscCheck(nconstraints == (active > 1 && nlocal ? expected + (has_const ? 1 : 0) + (dependent && step != 2 && step != 3 ? 0 : 1) : 0), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected number of constraints: %" PetscInt_FMT, nconstraints);
127: if (active > 1 && nlocal) {
128: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nlocal, &v));
129: PetscCall(VecGetArray(v, &values));
130: for (PetscInt i = 0; i < nlocal; i++) values[i] = i + 1;
131: PetscCall(VecRestoreArray(v, &values));
132: PetscCall(CheckConstraint(bddc->ConstraintMatrix, v));
133: if (has_const) {
134: PetscCall(VecSet(v, 1.0));
135: PetscCall(CheckConstraint(bddc->ConstraintMatrix, v));
136: }
137: for (PetscInt k = offset; k < offset + expected; k++) {
138: PetscCall(VecGetArray(v, &values));
139: for (PetscInt i = 0; i < nlocal; i++) values[i] = (i == 2 * k ? 1.0 : 0.0) - (i == 2 * k + 1 ? 1.0 : 0.0);
140: if (dependent && !k)
141: for (PetscInt i = 0; i < nlocal; i++) values[i] = i + 1;
142: PetscCall(VecRestoreArray(v, &values));
143: PetscCall(CheckConstraint(bddc->ConstraintMatrix, v));
144: }
145: PetscCall(VecDestroy(&v));
146: }
147: }
149: // Exercise the constant flag of the flux space independently of its computation flag.
150: PetscCall(MatNullSpaceDestroy(&bddc->nonetflux));
151: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &bddc->nonetflux));
152: for (PetscInt compute = 0; compute < 2; compute++) {
153: bddc->compute_nonetflux = (PetscBool)compute;
154: PetscCall(MatNullSpaceDestroy(&nsp));
155: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, constant, nmodes, modes, &nsp));
156: PetscCall(MatSetNearNullSpace(A, nsp));
157: PetscCall(MatScale(A, 1.01));
158: PetscCall(MatMult(A, exact, b));
159: PetscCall(KSPSetOperators(ksp, A, A));
160: PetscCall(KSPSolve(ksp, b, x));
161: PetscCall(MatGetSize(bddc->ConstraintMatrix, &nconstraints, NULL));
162: PetscCheck(nconstraints == (active > 1 && nlocal ? nmodes + 1 : 0), PETSC_COMM_SELF, PETSC_ERR_PLIB, "The constant from the flux space was lost: %" PetscInt_FMT " constraints", nconstraints);
163: if (active > 1 && nlocal) {
164: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nlocal, &v));
165: PetscCall(VecSet(v, 1.0));
166: PetscCall(CheckConstraint(bddc->ConstraintMatrix, v));
167: PetscCall(VecDestroy(&v));
168: }
169: }
171: PetscCall(KSPDestroy(&ksp));
172: PetscCall(MatNullSpaceDestroy(&nsp));
173: PetscCall(VecDestroy(&modes[0]));
174: PetscCall(VecDestroy(&modes[1]));
175: PetscCall(VecDestroy(&x));
176: PetscCall(VecDestroy(&b));
177: PetscCall(VecDestroy(&exact));
178: PetscCall(MatDestroy(&A));
179: PetscCall(MatDestroy(&B));
180: PetscCall(ISLocalToGlobalMappingDestroy(&map));
181: PetscCall(ISLocalToGlobalMappingDestroy(&pmap));
182: PetscCall(PetscFinalize());
183: return 0;
184: }
186: /*TEST
188: testset:
189: requires: double
190: output_file: output/empty.out
191: args: -pc_bddc_use_change_of_basis 0 -pc_bddc_coarsening_ratio 1 -ksp_error_if_not_converged -ksp_rtol 1e-12
192: test:
193: suffix: modes
194: nsize: {{1 2 3}}
195: args: -user_modes {{0 1 2}} -constant {{0 1}}
196: test:
197: suffix: dependent
198: nsize: 2
199: args: -dependent
200: test:
201: suffix: empty
202: nsize: 3
203: args: -empty_rank -constant {{0 1}}
204: test:
205: suffix: reset
206: nsize: {{1 2 3}}
207: args: -reset -constant
208: test:
209: suffix: coarse_resize
210: nsize: {{2 3}}
211: args: -pc_bddc_coarsening_ratio 8 -reset {{0 1}}
213: TEST*/