Actual source code: ex2.c
1: static char help[] = "Tests BDDC Nedelec support and user-defined primal vertices.\n\n";
3: #include <petsc/private/pcbddcimpl.h>
5: int main(int argc, char **args)
6: {
7: Mat A, local, G;
8: ISLocalToGlobalMapping map;
9: IS primal = NULL, stored;
10: KSP ksp;
11: PC pc;
12: PC_BDDC *bddc;
13: Vec x, b, exact, residual;
14: PetscInt *indices, *eoffset, *voffset, *orders, *selected, *requested, *input;
15: PetscInt start, end, n, nedge, nlocal, nbase, nmesh, nv, nselected, nrequested = 0, ninput = 0, order = 1, copies = 1, field = PETSC_DECIDE, unselected = -1;
16: PetscMPIInt rank, size, active;
17: PetscBool *expected;
18: PetscBool local_primal = PETSC_FALSE, closed_loop = PETSC_FALSE, periodic = PETSC_FALSE, branch = PETSC_FALSE, mixed_order = PETSC_FALSE;
19: PetscBool faces = PETSC_FALSE, permuted = PETSC_FALSE, empty_rank = PETSC_FALSE, distributed_primal = PETSC_FALSE, no_primal = PETSC_FALSE, other_only = PETSC_FALSE;
20: PetscBool explicit_fields = PETSC_FALSE, gradient_global = PETSC_TRUE, conforming = PETSC_TRUE;
21: PetscBool dirichlet = PETSC_FALSE, neumann = PETSC_FALSE, nullspace = PETSC_FALSE, nullspace_explicit = PETSC_FALSE;
22: PetscBool near_nullspace = PETSC_FALSE, near_nullspace_explicit = PETSC_FALSE, transpose = PETSC_TRUE, setprimal = PETSC_FALSE, check_multilevel = PETSC_FALSE, flg;
24: PetscFunctionBeginUser;
25: PetscCall(PetscInitialize(&argc, &args, NULL, help));
26: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
27: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
28: PetscCall(PetscOptionsGetBool(NULL, NULL, "-local_primal", &local_primal, NULL));
29: PetscCall(PetscOptionsGetBool(NULL, NULL, "-closed_loop", &closed_loop, NULL));
30: PetscCall(PetscOptionsGetBool(NULL, NULL, "-periodic", &periodic, NULL));
31: PetscCall(PetscOptionsGetBool(NULL, NULL, "-branch", &branch, NULL));
32: PetscCall(PetscOptionsGetBool(NULL, NULL, "-mixed_order", &mixed_order, NULL));
33: PetscCall(PetscOptionsGetBool(NULL, NULL, "-faces", &faces, NULL));
34: PetscCall(PetscOptionsGetBool(NULL, NULL, "-permuted", &permuted, NULL));
35: PetscCall(PetscOptionsGetBool(NULL, NULL, "-empty_rank", &empty_rank, NULL));
36: PetscCall(PetscOptionsGetBool(NULL, NULL, "-distributed_primal", &distributed_primal, NULL));
37: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_primal", &no_primal, NULL));
38: PetscCall(PetscOptionsGetBool(NULL, NULL, "-other_only", &other_only, NULL));
39: PetscCall(PetscOptionsGetBool(NULL, NULL, "-explicit_fields", &explicit_fields, NULL));
40: PetscCall(PetscOptionsGetBool(NULL, NULL, "-gradient_global", &gradient_global, NULL));
41: PetscCall(PetscOptionsGetBool(NULL, NULL, "-conforming", &conforming, NULL));
42: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dirichlet", &dirichlet, NULL));
43: PetscCall(PetscOptionsGetBool(NULL, NULL, "-neumann", &neumann, NULL));
44: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nullspace", &nullspace, NULL));
45: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nullspace_explicit", &nullspace_explicit, NULL));
46: PetscCall(PetscOptionsGetBool(NULL, NULL, "-near_nullspace", &near_nullspace, NULL));
47: PetscCall(PetscOptionsGetBool(NULL, NULL, "-near_nullspace_explicit", &near_nullspace_explicit, NULL));
48: PetscCall(PetscOptionsGetBool(NULL, NULL, "-transpose", &transpose, NULL));
49: PetscCall(PetscOptionsGetBool(NULL, NULL, "-check_multilevel", &check_multilevel, NULL));
50: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc_bddc_nedelec_field_primal", &setprimal, NULL));
51: PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &order, NULL));
52: PetscCall(PetscOptionsGetInt(NULL, NULL, "-copies", &copies, NULL));
53: PetscCall(PetscOptionsGetInt(NULL, NULL, "-field", &field, NULL));
54: PetscCall(PetscOptionsGetInt(NULL, NULL, "-unselected_edge", &unselected, NULL));
55: PetscCheck(order >= 1 && order <= 3, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Use orders 1 through 3");
56: PetscCheck(copies > 0, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "The number of local subdomains must be positive");
57: active = size - (empty_rank ? 1 : 0);
58: PetscCheck(active >= 2, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "At least two nonempty subdomains are required");
59: if (!gradient_global) explicit_fields = PETSC_TRUE;
61: /* Seven fine edges form a path, ring, branched tree, or path plus a separate triangle.
62: Optional face edges are shared only by ranks 0 and 1 and meet the path at node 2. */
63: nmesh = faces ? 11 : 7;
64: nv = faces ? 12 : 8;
65: PetscCall(PetscMalloc3(nmesh + 1, &eoffset, nmesh + 1, &voffset, nmesh, &orders));
66: eoffset[0] = 0;
67: for (PetscInt e = 0; e < nmesh; e++) {
68: orders[e] = mixed_order ? 1 + e % 3 : order;
69: eoffset[e + 1] = eoffset[e] + orders[e];
70: voffset[e] = nv;
71: nv += orders[e] - 1;
72: }
73: nedge = eoffset[nmesh];
74: PetscCheck(unselected >= -1 && unselected < nmesh, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Invalid unselected mesh edge");
75: if (unselected >= 0) unselected = eoffset[unselected] + orders[unselected] / 2;
77: /* Local numbering interleaves fields and may differ on each rank. */
78: nbase = rank < active ? (rank < 2 ? nedge : eoffset[7]) + 4 : 0;
79: nlocal = copies * nbase;
80: PetscCall(PetscMalloc1(nlocal, &indices));
81: if (nlocal) {
82: PetscInt ne = nbase - 4;
84: indices[0] = nedge + 2;
85: for (PetscInt i = 0; i < ne; i++) indices[i + 1] = ne - i - 1;
86: indices[ne + 1] = nedge;
87: indices[ne + 2] = nedge + 1;
88: indices[ne + 3] = nedge + 3 + rank;
89: if (permuted) {
90: for (PetscInt i = 0; i < nbase; i++) {
91: PetscInt j = (3 * i + rank + 1) % nbase, tmp = indices[i];
93: indices[i] = indices[j];
94: indices[j] = tmp;
95: }
96: }
97: for (PetscInt c = 1; c < copies; c++) PetscCall(PetscArraycpy(indices + c * nbase, indices, nbase));
98: }
99: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, nlocal, indices, PETSC_COPY_VALUES, &map));
100: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
101: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, nedge + 3 + active, nedge + 3 + active));
102: PetscCall(MatSetType(A, MATIS));
103: PetscCall(MatISSetAllowRepeated(A, (PetscBool)(copies > 1)));
104: PetscCall(MatSetLocalToGlobalMapping(A, map, map));
105: PetscCall(MatISGetLocalMat(A, &local));
106: PetscCall(MatSeqAIJSetPreallocation(local, nbase, NULL));
107: for (PetscInt i = 0; i < nlocal; i++) {
108: for (PetscInt j = i / nbase * nbase; j < (i / nbase + 1) * nbase; j++) PetscCall(MatSetValue(local, i, j, i == j ? 2.0 + 0.1 * i : -1.0 / nbase, INSERT_VALUES));
109: }
110: if (copies > 1) {
111: PetscInt *sizes;
113: PetscCall(PetscMalloc1(copies, &sizes));
114: for (PetscInt c = 0; c < copies; c++) sizes[c] = nbase;
115: PetscCall(MatSetVariableBlockSizes(local, copies, sizes));
116: PetscCall(PetscFree(sizes));
117: }
118: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
119: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
120: if (dirichlet) {
121: PetscInt *rows;
123: PetscCall(PetscMalloc1(orders[0], &rows));
124: for (PetscInt i = 0; i < orders[0]; i++) rows[i] = i;
125: PetscCall(MatZeroRowsColumns(A, orders[0], rows, 1.0, NULL, NULL));
126: PetscCall(PetscFree(rows));
127: }
129: /* Differentiate degree-p Lagrange polynomials at p points on each mesh edge.
130: All rows on an edge have the same p+1 nodal columns, without stored zeros. */
131: PetscCall(MatGetLocalSize(A, &n, NULL));
132: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, gradient_global ? n : PETSC_DECIDE, PETSC_DECIDE, gradient_global ? nedge + 3 + active : nedge, nv, 4, NULL, 4, NULL, &G));
133: PetscCall(MatGetOwnershipRange(G, &start, &end));
134: for (PetscInt e = 0; e < nmesh; e++) {
135: PetscInt p = orders[e], v0 = e, v1 = e + 1;
137: if (closed_loop && e >= 4 && e < 7) {
138: v0 = e + 1;
139: v1 = e == 6 ? 5 : e + 2;
140: }
141: if (periodic && e == 6) v1 = 0;
142: if (branch && e >= 4 && e < 7) v0 = e == 4 ? 2 : e;
143: if (e >= 7) {
144: v0 = e == 7 ? 2 : e;
145: v1 = e + 1;
146: }
147: for (PetscInt r = 0; r < p; r++) {
148: PetscInt row = eoffset[e] + r;
149: PetscReal t = p == 1 ? 0.5 : (PetscReal)r / (p - 1);
151: if (row < start || row >= end) continue;
152: for (PetscInt j = 0; j <= p; j++) {
153: PetscInt col = j == 0 ? v0 : (j == p ? v1 : voffset[e] + j - 1);
154: PetscReal value = 0.0, denominator = 1.0;
156: for (PetscInt k = 0; k <= p; k++)
157: if (k != j) denominator *= (PetscReal)(j - k) / p;
158: for (PetscInt k = 0; k <= p; k++) {
159: PetscReal term = 1.0;
161: if (k == j) continue;
162: for (PetscInt l = 0; l <= p; l++)
163: if (l != j && l != k) term *= t - (PetscReal)l / p;
164: value += term / denominator;
165: }
166: PetscCheck(value != 0.0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Zero entry in test gradient");
167: PetscCall(MatSetValue(G, row, col, value, INSERT_VALUES));
168: }
169: }
170: }
171: PetscCall(MatAssemblyBegin(G, MAT_FINAL_ASSEMBLY));
172: PetscCall(MatAssemblyEnd(G, MAT_FINAL_ASSEMBLY));
173: if (nullspace || nullspace_explicit) {
174: MatNullSpace nsp;
175: Vec nodes, edges;
176: PetscReal norm;
178: PetscCall(MatCreateVecs(G, &nodes, &edges));
179: PetscCall(VecSet(nodes, 1.0));
180: PetscCall(MatMult(G, nodes, edges));
181: PetscCall(VecNorm(edges, NORM_INFINITY, &norm));
182: PetscCheck(norm < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Test gradient does not annihilate constants");
183: if (nullspace_explicit) {
184: PetscCall(VecNormalize(nodes, NULL));
185: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, 1, &nodes, &nsp));
186: } else PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp));
187: PetscCall(VecDestroy(&nodes));
188: PetscCall(VecDestroy(&edges));
189: PetscCall(MatSetNullSpace(G, nsp));
190: PetscCall(MatNullSpaceDestroy(&nsp));
191: }
192: if (near_nullspace || near_nullspace_explicit) {
193: MatNullSpace nsp;
195: if (near_nullspace_explicit) {
196: Vec mode;
198: PetscCall(MatCreateVecs(A, &mode, NULL));
199: PetscCall(VecSet(mode, 1.0));
200: PetscCall(VecNormalize(mode, NULL));
201: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, 1, &mode, &nsp));
202: PetscCall(VecDestroy(&mode));
203: } else PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp));
204: PetscCall(MatSetNearNullSpace(A, nsp));
205: PetscCall(MatNullSpaceDestroy(&nsp));
206: }
208: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
209: PetscCall(KSPSetOperators(ksp, A, A));
210: PetscCall(KSPGetPC(ksp, &pc));
211: PetscCall(PCSetType(pc, PCBDDC));
212: if (explicit_fields) {
213: IS fields[2];
214: PetscInt *rows;
215: PetscInt counts[2] = {0, 0};
217: PetscCall(PetscMalloc1(nlocal, &rows));
218: for (PetscInt i = 0; i < nlocal; i++)
219: if (indices[i] >= nedge) rows[counts[0]++] = i;
220: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, counts[0], rows, PETSC_COPY_VALUES, &fields[0]));
221: for (PetscInt i = 0; i < nlocal; i++)
222: if (indices[i] < nedge) rows[counts[1]++] = i;
223: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, counts[1], rows, PETSC_COPY_VALUES, &fields[1]));
224: PetscCall(PCBDDCSetDofsSplittingLocal(pc, 2, fields));
225: PetscCall(ISDestroy(&fields[0]));
226: PetscCall(ISDestroy(&fields[1]));
227: PetscCall(PetscFree(rows));
228: field = 1;
229: }
230: PetscCall(PCBDDCSetDiscreteGradient(pc, G, mixed_order ? 0 : order, field, gradient_global, conforming));
231: if (dirichlet || neumann) {
232: IS boundary;
233: PetscInt *rows;
234: PetscInt nr = 0;
236: PetscCall(PetscMalloc1(nlocal, &rows));
237: if (dirichlet) {
238: for (PetscInt i = 0; i < nlocal; i++)
239: if (indices[i] < eoffset[1]) rows[nr++] = i;
240: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, rows, PETSC_COPY_VALUES, &boundary));
241: PetscCall(PCBDDCSetDirichletBoundariesLocal(pc, boundary));
242: PetscCall(ISDestroy(&boundary));
243: }
244: if (neumann) {
245: nr = 0;
246: for (PetscInt i = 0; i < nlocal; i++)
247: if (indices[i] >= eoffset[7] && indices[i] < nedge) rows[nr++] = i;
248: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, nr, rows, PETSC_COPY_VALUES, &boundary));
249: PetscCall(PCBDDCSetNeumannBoundariesLocal(pc, boundary));
250: PetscCall(ISDestroy(&boundary));
251: }
252: PetscCall(PetscFree(rows));
253: }
255: /* Include another field, select a non-first edge dof, and allow unsorted duplicates.
256: Global input may come from a process that does not own any of the selected dofs. */
257: PetscCall(PetscMalloc1(nedge + 3, &selected));
258: nselected = nedge + 3;
259: PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-primal_edges", selected, &nselected, &flg));
260: if (!flg) {
261: selected[0] = 3;
262: selected[1] = 8;
263: nselected = faces ? 2 : 1;
264: }
265: if (other_only) nselected = 0;
266: PetscCall(PetscMalloc2(nselected + 1, &requested, nselected + 1, &input));
267: PetscCall(PetscCalloc1(nedge + 3, &expected));
268: if (!no_primal) {
269: requested[nrequested++] = nedge;
270: expected[nedge] = PETSC_TRUE;
271: for (PetscInt i = 0; i < nselected; i++) {
272: PetscInt e = selected[i], row;
274: PetscCheck(e >= 0 && e < nmesh, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Invalid mesh edge %" PetscInt_FMT, e);
275: row = eoffset[e] + orders[e] / 2;
276: requested[nrequested++] = row;
277: expected[row] = PETSC_TRUE;
278: if (copies * active > 2 && !setprimal)
279: for (PetscInt j = eoffset[e]; j < eoffset[e + 1]; j++) expected[j] = PETSC_TRUE;
280: }
281: for (PetscInt i = 0; i < nrequested; i++) {
282: PetscInt row = requested[i], localrow;
283: PetscMPIInt contributor = distributed_primal ? (local_primal ? (i % 2 ? 1 : 0) : size - 1) : 0;
285: if (rank != contributor) continue;
286: if (local_primal) {
287: PetscCall(ISGlobalToLocalMappingApply(map, IS_GTOLM_MASK, 1, &row, NULL, &localrow));
288: PetscCheck(localrow >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Contributor does not own the requested primal dof");
289: input[ninput++] = localrow;
290: } else input[ninput++] = row;
291: }
292: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ninput, input, PETSC_COPY_VALUES, &primal));
293: if (local_primal) PetscCall(PCBDDCSetPrimalVerticesLocalIS(pc, primal));
294: else PetscCall(PCBDDCSetPrimalVerticesIS(pc, primal));
295: }
296: if (setprimal && copies * active > 2)
297: for (PetscInt i = 0; i < eoffset[7]; i++) expected[i] = PETSC_TRUE;
298: PetscCall(KSPSetFromOptions(ksp));
299: PetscCall(KSPSetUp(ksp));
301: if (primal && !local_primal) {
302: PetscCall(PCBDDCGetPrimalVerticesIS(pc, &stored));
303: PetscCheck(stored, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Global user primal vertices were lost during setup");
304: PetscCall(ISEqual(primal, stored, &flg));
305: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Global user primal vertices changed during setup");
306: }
307: PetscCall(ISDestroy(&primal));
309: /* Check the stored vertices and actual point constraints in each rank's numbering. */
310: PetscCall(PCBDDCGetPrimalVerticesLocalIS(pc, &stored));
311: bddc = (PC_BDDC *)pc->data;
312: if (check_multilevel) {
313: if (bddc->coarse_ksp) {
314: PC coarsepc;
315: PetscBool isbddc;
317: PetscCall(KSPGetPC(bddc->coarse_ksp, &coarsepc));
318: PetscCall(PetscObjectTypeCompare((PetscObject)coarsepc, PCBDDC, &isbddc));
319: PetscCheck(isbddc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing coarse BDDC level");
320: PetscCheck(((PC_BDDC *)coarsepc->data)->discretegradient, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing coarse discrete gradient");
321: }
322: PetscCheck(bddc->nedcG, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing generated coarse gradient");
323: }
324: PetscCall(MatGetSize(bddc->ConstraintMatrix, &n, NULL));
325: for (PetscInt localrow = 0; localrow < nlocal; localrow++) {
326: PetscInt i = indices[localrow], pos;
327: PetscBool found = PETSC_FALSE;
329: if (i >= nedge + 3 || !expected[i]) continue;
330: PetscCall(ISLocate(stored, localrow, &pos));
331: PetscCheck(pos >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Primal dof %" PetscInt_FMT " was lost", i);
332: for (PetscInt j = 0; j < n; j++) {
333: const PetscInt *cols;
334: const PetscScalar *vals;
335: PetscInt nc;
337: PetscCall(MatGetRow(bddc->ConstraintMatrix, j, &nc, &cols, &vals));
338: if (nc == 1 && cols[0] == localrow && vals[0] == 1.0) found = PETSC_TRUE;
339: PetscCall(MatRestoreRow(bddc->ConstraintMatrix, j, &nc, &cols, &vals));
340: }
341: PetscCheck(found, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing point constraint for primal dof %" PetscInt_FMT, i);
342: }
343: PetscCheck((PetscBool)!!bddc->user_ChangeOfBasisMatrix == (PetscBool)(copies * active > 2 && !setprimal), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected Nedelec change-of-basis path");
344: if (bddc->user_ChangeOfBasisMatrix) {
345: PetscCall(MatGetOwnershipRange(bddc->user_ChangeOfBasisMatrix, &start, &end));
346: for (PetscInt i = start; i < PetscMin(end, nedge + 3); i++) {
347: const PetscInt *cols;
348: const PetscScalar *vals;
349: PetscInt nc;
350: PetscBool identity = PETSC_TRUE, diagonal = PETSC_FALSE;
352: if (!expected[i] && i != unselected) continue;
353: PetscCall(MatGetRow(bddc->user_ChangeOfBasisMatrix, i, &nc, &cols, &vals));
354: for (PetscInt j = 0; j < nc; j++) {
355: if (cols[j] == i) {
356: diagonal = PETSC_TRUE;
357: if (PetscAbsScalar(vals[j] - 1.0) > PETSC_SMALL) identity = PETSC_FALSE;
358: } else if (PetscAbsScalar(vals[j]) > PETSC_SMALL) identity = PETSC_FALSE;
359: }
360: PetscCall(MatRestoreRow(bddc->user_ChangeOfBasisMatrix, i, &nc, &cols, &vals));
361: if (expected[i]) PetscCheck(identity && diagonal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Change of basis modifies primal dof %" PetscInt_FMT, i);
362: else PetscCheck(!identity || !diagonal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Change of basis leaves unselected edge dof %" PetscInt_FMT " unchanged", i);
363: }
364: }
366: /* Two nonconstant right-hand sides verify the solve and reuse of the setup. */
367: PetscCall(MatCreateVecs(A, &x, &b));
368: PetscCall(VecDuplicate(x, &exact));
369: PetscCall(VecDuplicate(x, &residual));
370: PetscCall(VecGetOwnershipRange(exact, &start, &end));
371: for (PetscInt solve = 0; solve < 2; solve++) {
372: PetscScalar *values;
373: PetscReal error, norm;
375: PetscCall(VecGetArray(exact, &values));
376: for (PetscInt i = start; i < end; i++) values[i - start] = PetscSinReal((solve + 1) * 0.17 * (i + 1)) + 0.03 * (i % 5);
377: PetscCall(VecRestoreArray(exact, &values));
378: PetscCall(KSPSetUp(ksp));
379: if (solve && transpose) {
380: PetscCall(MatMultTranspose(A, exact, b));
381: PetscCall(KSPSolveTranspose(ksp, b, x));
382: PetscCall(MatMultTranspose(A, x, residual));
383: } else {
384: PetscCall(MatMult(A, exact, b));
385: PetscCall(KSPSolve(ksp, b, x));
386: PetscCall(MatMult(A, x, residual));
387: }
388: PetscCall(VecAXPY(residual, -1.0, b));
389: PetscCall(VecNorm(residual, NORM_INFINITY, &error));
390: PetscCall(VecNorm(b, NORM_INFINITY, &norm));
391: PetscCheck(error < PETSC_SMALL * norm, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Relative residual %g", (double)(error / norm));
392: PetscCall(VecAXPY(x, -1.0, exact));
393: PetscCall(VecNorm(x, NORM_INFINITY, &error));
394: PetscCheck(error < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Solution error %g", (double)error);
395: }
397: PetscCall(PetscFree(indices));
398: PetscCall(PetscFree3(eoffset, voffset, orders));
399: PetscCall(PetscFree(selected));
400: PetscCall(PetscFree2(requested, input));
401: PetscCall(PetscFree(expected));
402: PetscCall(ISLocalToGlobalMappingDestroy(&map));
403: PetscCall(VecDestroy(&residual));
404: PetscCall(VecDestroy(&exact));
405: PetscCall(VecDestroy(&x));
406: PetscCall(VecDestroy(&b));
407: PetscCall(KSPDestroy(&ksp));
408: PetscCall(MatDestroy(&G));
409: PetscCall(MatDestroy(&A));
410: PetscCall(PetscFinalize());
411: return 0;
412: }
414: /*TEST
416: testset:
417: nsize: 3
418: requires: double
419: args: -local_primal {{0 1}} -ksp_error_if_not_converged -ksp_rtol 1e-12
420: output_file: output/empty.out
421: test:
422: suffix: path
423: args: -unselected_edge 1 -pc_bddc_nedelec_field_primal {{0 1}}
424: test:
425: suffix: loop
426: args: -closed_loop
427: test:
428: suffix: monolithic
429: args: -field 0
430: test:
431: suffix: quadratic
432: args: -order 2 -unselected_edge 1 -periodic {{0 1}} -pc_bddc_nedelec_order {{0 2}}
433: test:
434: suffix: cubic
435: args: -order 3 -permuted -distributed_primal -primal_edges 5,3,5,1,3
436: test:
437: suffix: mixed_order
438: args: -mixed_order -periodic -permuted -explicit_fields -distributed_primal
439: test:
440: suffix: fields
441: args: -order 2 -explicit_fields -gradient_global {{0 1}} -permuted
442: test:
443: suffix: faces
444: args: -order 2 -faces -permuted -distributed_primal -pc_bddc_nedelec_field_primal {{0 1}}
445: test:
446: suffix: boundaries
447: args: -order 2 -faces -dirichlet -neumann -permuted -conforming {{0 1}}
448: test:
449: suffix: branch
450: args: -order 2 -branch -permuted
451: test:
452: suffix: all_primal
453: args: -order 2 -primal_edges 6,4,2,0,5,3,1 -permuted -nullspace
454: test:
455: suffix: other_field
456: args: -order 2 -other_only -permuted -distributed_primal
457: test:
458: suffix: near_nullspace
459: args: -order 2 -near_nullspace -faces
460: test:
461: suffix: combined_nullspaces
462: args: -order 2 -nullspace -near_nullspace -near_nullspace_explicit {{0 1}}
463: test:
464: suffix: nullspace_fields
465: args: -order 2 -nullspace -nullspace_explicit {{0 1}} -explicit_fields -gradient_global {{0 1}} -permuted -distributed_primal
466: test:
467: suffix: multilevel
468: nsize: 4
469: args: -order 2 -nullspace -nullspace_explicit {{0 1}} -check_multilevel -pc_bddc_levels 2 -pc_bddc_coarse_eqs_limit 0 -pc_bddc_coarsening_ratio 2 -pc_bddc_aggregator_0_mat_partitioning_type average
470: test:
471: suffix: no_coarse_edges
472: nsize: 2
473: args: -order 2 -faces -permuted -distributed_primal
474: test:
475: suffix: empty_rank
476: nsize: 4
477: args: -order 2 -empty_rank -permuted -distributed_primal -nullspace
478: test:
479: suffix: repeated
480: args: -order 2 -copies 2 -explicit_fields -permuted -distributed_primal
482: test:
483: suffix: no_user
484: nsize: 3
485: requires: double
486: output_file: output/empty.out
487: args: -order 2 -no_primal -explicit_fields -nullspace {{0 1}} -ksp_error_if_not_converged -ksp_rtol 1e-12
489: TEST*/