Actual source code: ex5.c
1: static char help[] = "Tests BDDC coarse solver reuse when the coarse MPI communicator changes.\n\n";
3: #include <petsc/private/pcbddcimpl.h>
5: int main(int argc, char **args)
6: {
7: Mat A, local, coarse;
8: MatNullSpace nsp;
9: ISLocalToGlobalMapping map;
10: KSP ksp;
11: PC pc, coarsepc;
12: PC_BDDC *bddc;
13: Vec modes[3], x, b, exact;
14: PetscScalar *values;
15: PetscInt indices[] = {0, 1, 2, 3, 4, 5, 6, 7};
16: PetscInt counts[] = {2, 2, 1, 1, 3, 3, 1, 2};
17: PetscInt nlocal, start, end, participating, previous = -1, empty_ranks = 0;
18: PetscMPIInt rank, size, active, comparison;
19: PetscObjectId id, previous_id = 0;
20: PetscBool reset = PETSC_FALSE, changed = PETSC_FALSE, redundant;
21: PetscReal error;
23: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &args, NULL, help));
25: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
26: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
27: PetscCall(PetscOptionsGetInt(NULL, NULL, "-empty_ranks", &empty_ranks, NULL));
28: PetscCall(PetscOptionsGetBool(NULL, NULL, "-reset", &reset, NULL));
29: PetscCall(PetscMPIIntCast(size - empty_ranks, &active));
30: PetscCheck(active >= 2, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "At least two nonempty subdomains are required");
31: nlocal = rank < active ? 8 : 0;
32: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, nlocal, indices, PETSC_COPY_VALUES, &map));
33: PetscCall(MatCreateIS(PETSC_COMM_WORLD, 1, PETSC_DECIDE, PETSC_DECIDE, 8, 8, map, map, &A));
34: PetscCall(ISLocalToGlobalMappingDestroy(&map));
35: PetscCall(MatISSetPreallocation(A, 8, NULL, 8, NULL));
36: PetscCall(MatISGetLocalMat(A, &local));
37: for (PetscInt i = 0; i < nlocal; i++)
38: for (PetscInt j = 0; j < nlocal; j++) PetscCall(MatSetValue(local, i, j, i == j ? 2.0 : 0.1, INSERT_VALUES));
39: PetscCall(MatISRestoreLocalMat(A, &local));
40: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
41: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
42: PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));
43: PetscCall(MatCreateVecs(A, &x, &b));
44: PetscCall(VecDuplicate(x, &exact));
45: PetscCall(VecGetOwnershipRange(x, &start, &end));
46: for (PetscInt k = 0; k < 3; k++) {
47: PetscCall(VecDuplicate(x, &modes[k]));
48: PetscCall(VecGetArray(modes[k], &values));
49: 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);
50: PetscCall(VecRestoreArray(modes[k], &values));
51: PetscCall(VecNormalize(modes[k], NULL));
52: }
53: PetscCall(VecSet(exact, 1.0));
54: PetscCall(VecAXPY(exact, 0.5, modes[1]));
55: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
56: PetscCall(KSPSetType(ksp, KSPCG));
57: PetscCall(KSPGetPC(ksp, &pc));
58: PetscCall(PCSetType(pc, PCBDDC));
59: PetscCall(KSPSetFromOptions(ksp));
60: bddc = (PC_BDDC *)pc->data;
61: for (PetscInt step = 0; step < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(counts); step++) {
62: if (reset && step) PetscCall(PCReset(pc));
63: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, counts[step], modes, &nsp));
64: PetscCall(MatSetNearNullSpace(A, nsp));
65: PetscCall(MatNullSpaceDestroy(&nsp));
66: PetscCall(MatScale(A, 1.01));
67: PetscCall(MatMult(A, exact, b));
68: PetscCall(KSPSetOperators(ksp, A, A));
69: PetscCall(KSPSolve(ksp, b, x));
70: PetscCall(VecAXPY(x, -1.0, exact));
71: PetscCall(VecNorm(x, NORM_INFINITY, &error));
72: PetscCheck(error < PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Solution error after changing the coarse space: %g", (double)error);
73: PetscCheck(bddc->coarse_size == counts[step], PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected coarse size: %" PetscInt_FMT, bddc->coarse_size);
74: participating = (PetscInt)!!bddc->coarse_ksp;
75: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &participating, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
76: if (step && participating != previous) changed = PETSC_TRUE;
77: previous = participating;
78: if (bddc->coarse_ksp) {
79: PetscCall(PetscObjectGetId((PetscObject)bddc->coarse_ksp, &id));
80: PetscCheck(reset || !step || counts[step] != counts[step - 1] || id == previous_id, PETSC_COMM_SELF, PETSC_ERR_PLIB, "The coarse KSP was replaced despite unchanged participation");
81: previous_id = id;
82: PetscCall(KSPGetOperators(bddc->coarse_ksp, &coarse, NULL));
83: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)bddc->coarse_ksp), PetscObjectComm((PetscObject)coarse), &comparison));
84: PetscCheck(comparison == MPI_IDENT || comparison == MPI_CONGRUENT, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coarse KSP and matrix communicators differ");
85: PetscCall(KSPGetPC(bddc->coarse_ksp, &coarsepc));
86: PetscCall(PetscObjectTypeCompare((PetscObject)coarsepc, PCREDUNDANT, &redundant));
87: PetscCheck(redundant, PETSC_COMM_SELF, PETSC_ERR_PLIB, "The terminal coarse PC must remain PCREDUNDANT");
88: }
89: }
90: PetscCheck(changed, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The test did not change coarse participation");
91: PetscCall(KSPDestroy(&ksp));
92: for (PetscInt k = 0; k < 3; k++) PetscCall(VecDestroy(&modes[k]));
93: PetscCall(VecDestroy(&x));
94: PetscCall(VecDestroy(&b));
95: PetscCall(VecDestroy(&exact));
96: PetscCall(MatDestroy(&A));
97: PetscCall(PetscFinalize());
98: return 0;
99: }
101: /*TEST
103: testset:
104: requires: double
105: output_file: output/empty.out
106: args: -ksp_error_if_not_converged -ksp_rtol 1e-12
107: test:
108: suffix: resize
109: nsize: {{2 3}}
110: args: -pc_bddc_use_change_of_basis {{0 1}} -reset {{0 1}}
111: test:
112: suffix: empty
113: nsize: 4
114: args: -empty_ranks 2 -pc_bddc_coarse_eqs_per_proc 2 -pc_bddc_aggregator_0_mat_partitioning_type average -reset {{0 1}}
115: test:
116: suffix: deluxe
117: nsize: 3
118: args: -pc_bddc_use_deluxe_scaling -pc_bddc_use_change_of_basis {{0 1}}
120: TEST*/