Actual source code: ex71.c
1: static char help[] = "This example illustrates the use of PCBDDC/FETI-DP with 2D/3D DMDA.\n\
2: It solves the constant coefficient Poisson problem or the Elasticity problem \n\
3: on a uniform grid of [0,cells_x] x [0,cells_y] x [0,cells_z]\n\n";
5: /* Contributed by Wim Vanroose <wim@vanroo.se> */
7: #include <petscksp.h>
8: #include <petscpc.h>
9: #include <petscdm.h>
10: #include <petscdmda.h>
11: #include <petscdmplex.h>
13: static PetscScalar poiss_1D_emat[] = {1.0000000000000000e+00, -1.0000000000000000e+00, -1.0000000000000000e+00, 1.0000000000000000e+00};
14: static PetscScalar poiss_2D_emat[] = {6.6666666666666674e-01, -1.6666666666666666e-01, -1.6666666666666666e-01, -3.3333333333333337e-01, -1.6666666666666666e-01, 6.6666666666666674e-01, -3.3333333333333337e-01, -1.6666666666666666e-01,
15: -1.6666666666666666e-01, -3.3333333333333337e-01, 6.6666666666666674e-01, -1.6666666666666666e-01, -3.3333333333333337e-01, -1.6666666666666666e-01, -1.6666666666666666e-01, 6.6666666666666674e-01};
16: static PetscScalar poiss_3D_emat[] = {3.3333333333333348e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -8.3333333333333343e-02, 0.0000000000000000e+00, -8.3333333333333343e-02, -8.3333333333333343e-02, -8.3333333333333356e-02,
17: 0.0000000000000000e+00, 3.3333333333333337e-01, -8.3333333333333343e-02, 0.0000000000000000e+00, -8.3333333333333343e-02, 0.0000000000000000e+00, -8.3333333333333356e-02, -8.3333333333333343e-02,
18: 0.0000000000000000e+00, -8.3333333333333343e-02, 3.3333333333333337e-01, 0.0000000000000000e+00, -8.3333333333333343e-02, -8.3333333333333356e-02, 0.0000000000000000e+00, -8.3333333333333343e-02,
19: -8.3333333333333343e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, 3.3333333333333348e-01, -8.3333333333333356e-02, -8.3333333333333343e-02, -8.3333333333333343e-02, 0.0000000000000000e+00,
20: 0.0000000000000000e+00, -8.3333333333333343e-02, -8.3333333333333343e-02, -8.3333333333333356e-02, 3.3333333333333337e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -8.3333333333333343e-02,
21: -8.3333333333333343e-02, 0.0000000000000000e+00, -8.3333333333333356e-02, -8.3333333333333343e-02, 0.0000000000000000e+00, 3.3333333333333337e-01, -8.3333333333333343e-02, 0.0000000000000000e+00,
22: -8.3333333333333343e-02, -8.3333333333333356e-02, 0.0000000000000000e+00, -8.3333333333333343e-02, 0.0000000000000000e+00, -8.3333333333333343e-02, 3.3333333333333337e-01, 0.0000000000000000e+00,
23: -8.3333333333333356e-02, -8.3333333333333343e-02, -8.3333333333333343e-02, 0.0000000000000000e+00, -8.3333333333333343e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, 3.3333333333333337e-01};
24: static PetscScalar elast_1D_emat[] = {3.0000000000000000e+00, -3.0000000000000000e+00, -3.0000000000000000e+00, 3.0000000000000000e+00};
25: static PetscScalar elast_2D_emat[] = {1.3333333333333335e+00, 5.0000000000000000e-01, -8.3333333333333337e-01, 0.0000000000000000e+00, 1.6666666666666671e-01, 0.0000000000000000e+00, -6.6666666666666674e-01, -5.0000000000000000e-01,
26: 5.0000000000000000e-01, 1.3333333333333335e+00, 0.0000000000000000e+00, 1.6666666666666671e-01, 0.0000000000000000e+00, -8.3333333333333337e-01, -5.0000000000000000e-01, -6.6666666666666674e-01,
27: -8.3333333333333337e-01, 0.0000000000000000e+00, 1.3333333333333335e+00, -5.0000000000000000e-01, -6.6666666666666674e-01, 5.0000000000000000e-01, 1.6666666666666674e-01, 0.0000000000000000e+00,
28: 0.0000000000000000e+00, 1.6666666666666671e-01, -5.0000000000000000e-01, 1.3333333333333335e+00, 5.0000000000000000e-01, -6.6666666666666674e-01, 0.0000000000000000e+00, -8.3333333333333337e-01,
29: 1.6666666666666671e-01, 0.0000000000000000e+00, -6.6666666666666674e-01, 5.0000000000000000e-01, 1.3333333333333335e+00, -5.0000000000000000e-01, -8.3333333333333337e-01, 0.0000000000000000e+00,
30: 0.0000000000000000e+00, -8.3333333333333337e-01, 5.0000000000000000e-01, -6.6666666666666674e-01, -5.0000000000000000e-01, 1.3333333333333335e+00, 0.0000000000000000e+00, 1.6666666666666674e-01,
31: -6.6666666666666674e-01, -5.0000000000000000e-01, 1.6666666666666674e-01, 0.0000000000000000e+00, -8.3333333333333337e-01, 0.0000000000000000e+00, 1.3333333333333335e+00, 5.0000000000000000e-01,
32: -5.0000000000000000e-01, -6.6666666666666674e-01, 0.0000000000000000e+00, -8.3333333333333337e-01, 0.0000000000000000e+00, 1.6666666666666674e-01, 5.0000000000000000e-01, 1.3333333333333335e+00};
33: static PetscScalar elast_3D_emat[] =
34: {5.5555555555555558e-01, 1.6666666666666666e-01, 1.6666666666666666e-01, -2.2222222222222232e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 1.1111111111111113e-01, 0.0000000000000000e+00, 8.3333333333333356e-02,
35: -1.9444444444444442e-01, -1.6666666666666669e-01, 0.0000000000000000e+00, 1.1111111111111112e-01, 8.3333333333333356e-02, 0.0000000000000000e+00, -1.9444444444444445e-01, 0.0000000000000000e+00, -1.6666666666666669e-01,
36: -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.3888888888888887e-01, -8.3333333333333356e-02, -8.3333333333333356e-02, 1.6666666666666666e-01, 5.5555555555555558e-01, 1.6666666666666666e-01,
37: 0.0000000000000000e+00, 1.1111111111111113e-01, 8.3333333333333356e-02, 0.0000000000000000e+00, -2.2222222222222232e-01, 0.0000000000000000e+00, -1.6666666666666669e-01, -1.9444444444444442e-01, 0.0000000000000000e+00,
38: 8.3333333333333356e-02, 1.1111111111111112e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.9444444444444445e-01, -1.6666666666666669e-01,
39: -8.3333333333333356e-02, -1.3888888888888887e-01, -8.3333333333333356e-02, 1.6666666666666666e-01, 1.6666666666666666e-01, 5.5555555555555558e-01, 0.0000000000000000e+00, 8.3333333333333356e-02, 1.1111111111111112e-01,
40: 8.3333333333333356e-02, 0.0000000000000000e+00, 1.1111111111111112e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222229e-01,
41: -1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444445e-01, 0.0000000000000000e+00, -1.6666666666666669e-01, -1.9444444444444445e-01, -8.3333333333333356e-02, -8.3333333333333356e-02, -1.3888888888888887e-01,
42: -2.2222222222222232e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 5.5555555555555558e-01, -1.6666666666666666e-01, -1.6666666666666666e-01, -1.9444444444444442e-01, 1.6666666666666669e-01, 0.0000000000000000e+00,
43: 1.1111111111111113e-01, 0.0000000000000000e+00, -8.3333333333333356e-02, -1.9444444444444445e-01, 0.0000000000000000e+00, 1.6666666666666669e-01, 1.1111111111111113e-01, -8.3333333333333356e-02, 0.0000000000000000e+00,
44: -1.3888888888888887e-01, 8.3333333333333356e-02, 8.3333333333333356e-02, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 1.1111111111111113e-01, 8.3333333333333356e-02,
45: -1.6666666666666666e-01, 5.5555555555555558e-01, 1.6666666666666669e-01, 1.6666666666666669e-01, -1.9444444444444442e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222229e-01, 0.0000000000000000e+00,
46: 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00, -8.3333333333333356e-02, 1.1111111111111112e-01, 0.0000000000000000e+00, 8.3333333333333356e-02, -1.3888888888888887e-01, -8.3333333333333356e-02,
47: 0.0000000000000000e+00, -1.9444444444444448e-01, -1.6666666666666666e-01, 0.0000000000000000e+00, 8.3333333333333356e-02, 1.1111111111111112e-01, -1.6666666666666666e-01, 1.6666666666666669e-01, 5.5555555555555558e-01,
48: 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, -8.3333333333333356e-02, 0.0000000000000000e+00, 1.1111111111111112e-01, 1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444445e-01,
49: 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01, 8.3333333333333356e-02, -8.3333333333333356e-02, -1.3888888888888887e-01, 0.0000000000000000e+00, -1.6666666666666666e-01, -1.9444444444444448e-01,
50: 1.1111111111111113e-01, 0.0000000000000000e+00, 8.3333333333333356e-02, -1.9444444444444442e-01, 1.6666666666666669e-01, 0.0000000000000000e+00, 5.5555555555555569e-01, -1.6666666666666666e-01, 1.6666666666666669e-01,
51: -2.2222222222222229e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.3888888888888887e-01, 8.3333333333333356e-02, -8.3333333333333356e-02,
52: 1.1111111111111112e-01, -8.3333333333333343e-02, 0.0000000000000000e+00, -1.9444444444444448e-01, 0.0000000000000000e+00, -1.6666666666666669e-01, 0.0000000000000000e+00, -2.2222222222222232e-01, 0.0000000000000000e+00,
53: 1.6666666666666669e-01, -1.9444444444444442e-01, 0.0000000000000000e+00, -1.6666666666666666e-01, 5.5555555555555558e-01, -1.6666666666666669e-01, 0.0000000000000000e+00, 1.1111111111111113e-01, -8.3333333333333343e-02,
54: 0.0000000000000000e+00, -1.9444444444444445e-01, 1.6666666666666669e-01, 8.3333333333333356e-02, -1.3888888888888887e-01, 8.3333333333333356e-02, -8.3333333333333343e-02, 1.1111111111111113e-01, 0.0000000000000000e+00,
55: 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00, 8.3333333333333356e-02, 0.0000000000000000e+00, 1.1111111111111112e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02,
56: 1.6666666666666669e-01, -1.6666666666666669e-01, 5.5555555555555558e-01, 0.0000000000000000e+00, -8.3333333333333343e-02, 1.1111111111111112e-01, 0.0000000000000000e+00, 1.6666666666666669e-01, -1.9444444444444445e-01,
57: -8.3333333333333356e-02, 8.3333333333333356e-02, -1.3888888888888887e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01, -1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444448e-01,
58: -1.9444444444444442e-01, -1.6666666666666669e-01, 0.0000000000000000e+00, 1.1111111111111113e-01, 0.0000000000000000e+00, -8.3333333333333356e-02, -2.2222222222222229e-01, 0.0000000000000000e+00, 0.0000000000000000e+00,
59: 5.5555555555555558e-01, 1.6666666666666669e-01, -1.6666666666666666e-01, -1.3888888888888887e-01, -8.3333333333333356e-02, 8.3333333333333356e-02, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00,
60: -1.9444444444444448e-01, 0.0000000000000000e+00, 1.6666666666666669e-01, 1.1111111111111112e-01, 8.3333333333333343e-02, 0.0000000000000000e+00, -1.6666666666666669e-01, -1.9444444444444442e-01, 0.0000000000000000e+00,
61: 0.0000000000000000e+00, -2.2222222222222229e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 1.1111111111111113e-01, -8.3333333333333343e-02, 1.6666666666666669e-01, 5.5555555555555558e-01, -1.6666666666666669e-01,
62: -8.3333333333333356e-02, -1.3888888888888887e-01, 8.3333333333333356e-02, 0.0000000000000000e+00, -1.9444444444444448e-01, 1.6666666666666669e-01, 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00,
63: 8.3333333333333343e-02, 1.1111111111111112e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, -8.3333333333333356e-02, 0.0000000000000000e+00, 1.1111111111111112e-01,
64: 0.0000000000000000e+00, -8.3333333333333343e-02, 1.1111111111111112e-01, -1.6666666666666666e-01, -1.6666666666666669e-01, 5.5555555555555558e-01, 8.3333333333333356e-02, 8.3333333333333356e-02, -1.3888888888888887e-01,
65: 0.0000000000000000e+00, 1.6666666666666669e-01, -1.9444444444444448e-01, 1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444448e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01,
66: 1.1111111111111112e-01, 8.3333333333333356e-02, 0.0000000000000000e+00, -1.9444444444444445e-01, 0.0000000000000000e+00, 1.6666666666666669e-01, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00,
67: -1.3888888888888887e-01, -8.3333333333333356e-02, 8.3333333333333356e-02, 5.5555555555555569e-01, 1.6666666666666669e-01, -1.6666666666666669e-01, -2.2222222222222227e-01, 0.0000000000000000e+00, 0.0000000000000000e+00,
68: 1.1111111111111112e-01, 0.0000000000000000e+00, -8.3333333333333343e-02, -1.9444444444444448e-01, -1.6666666666666669e-01, 0.0000000000000000e+00, 8.3333333333333356e-02, 1.1111111111111112e-01, 0.0000000000000000e+00,
69: 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.9444444444444445e-01, 1.6666666666666669e-01, -8.3333333333333356e-02, -1.3888888888888887e-01, 8.3333333333333356e-02,
70: 1.6666666666666669e-01, 5.5555555555555558e-01, -1.6666666666666669e-01, 0.0000000000000000e+00, 1.1111111111111112e-01, -8.3333333333333343e-02, 0.0000000000000000e+00, -2.2222222222222227e-01, 0.0000000000000000e+00,
71: -1.6666666666666669e-01, -1.9444444444444448e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222229e-01, 1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444445e-01,
72: 0.0000000000000000e+00, 1.6666666666666669e-01, -1.9444444444444445e-01, 8.3333333333333356e-02, 8.3333333333333356e-02, -1.3888888888888887e-01, -1.6666666666666669e-01, -1.6666666666666669e-01, 5.5555555555555558e-01,
73: 0.0000000000000000e+00, -8.3333333333333343e-02, 1.1111111111111113e-01, -8.3333333333333343e-02, 0.0000000000000000e+00, 1.1111111111111113e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02,
74: -1.9444444444444445e-01, 0.0000000000000000e+00, -1.6666666666666669e-01, 1.1111111111111113e-01, -8.3333333333333356e-02, 0.0000000000000000e+00, -1.3888888888888887e-01, 8.3333333333333356e-02, -8.3333333333333356e-02,
75: -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 5.5555555555555558e-01, -1.6666666666666669e-01, 1.6666666666666669e-01,
76: -1.9444444444444448e-01, 1.6666666666666669e-01, 0.0000000000000000e+00, 1.1111111111111112e-01, 0.0000000000000000e+00, 8.3333333333333343e-02, 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00,
77: -8.3333333333333356e-02, 1.1111111111111112e-01, 0.0000000000000000e+00, 8.3333333333333356e-02, -1.3888888888888887e-01, 8.3333333333333356e-02, 0.0000000000000000e+00, -1.9444444444444448e-01, 1.6666666666666669e-01,
78: 0.0000000000000000e+00, 1.1111111111111112e-01, -8.3333333333333343e-02, -1.6666666666666669e-01, 5.5555555555555558e-01, -1.6666666666666666e-01, 1.6666666666666669e-01, -1.9444444444444448e-01, 0.0000000000000000e+00,
79: 0.0000000000000000e+00, -2.2222222222222227e-01, 0.0000000000000000e+00, -1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444445e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01,
80: -8.3333333333333356e-02, 8.3333333333333356e-02, -1.3888888888888887e-01, 0.0000000000000000e+00, 1.6666666666666669e-01, -1.9444444444444448e-01, 0.0000000000000000e+00, -8.3333333333333343e-02, 1.1111111111111113e-01,
81: 1.6666666666666669e-01, -1.6666666666666666e-01, 5.5555555555555558e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, 8.3333333333333343e-02, 0.0000000000000000e+00, 1.1111111111111113e-01,
82: -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.3888888888888887e-01, 8.3333333333333356e-02, 8.3333333333333356e-02, 1.1111111111111112e-01, -8.3333333333333343e-02, 0.0000000000000000e+00,
83: -1.9444444444444448e-01, 0.0000000000000000e+00, 1.6666666666666669e-01, 1.1111111111111112e-01, 0.0000000000000000e+00, -8.3333333333333343e-02, -1.9444444444444448e-01, 1.6666666666666669e-01, 0.0000000000000000e+00,
84: 5.5555555555555558e-01, -1.6666666666666669e-01, -1.6666666666666669e-01, -2.2222222222222227e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.9444444444444445e-01, -1.6666666666666669e-01,
85: 8.3333333333333356e-02, -1.3888888888888887e-01, -8.3333333333333356e-02, -8.3333333333333343e-02, 1.1111111111111113e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00,
86: 0.0000000000000000e+00, -2.2222222222222227e-01, 0.0000000000000000e+00, 1.6666666666666669e-01, -1.9444444444444448e-01, 0.0000000000000000e+00, -1.6666666666666669e-01, 5.5555555555555558e-01, 1.6666666666666669e-01,
87: 0.0000000000000000e+00, 1.1111111111111112e-01, 8.3333333333333343e-02, 0.0000000000000000e+00, -1.6666666666666669e-01, -1.9444444444444445e-01, 8.3333333333333356e-02, -8.3333333333333356e-02, -1.3888888888888887e-01,
88: 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01, 1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444448e-01, -8.3333333333333343e-02, 0.0000000000000000e+00, 1.1111111111111113e-01,
89: 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02, -1.6666666666666669e-01, 1.6666666666666669e-01, 5.5555555555555558e-01, 0.0000000000000000e+00, 8.3333333333333343e-02, 1.1111111111111113e-01,
90: -1.3888888888888887e-01, -8.3333333333333356e-02, -8.3333333333333356e-02, -2.7777777777777769e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, -1.9444444444444448e-01, 0.0000000000000000e+00, -1.6666666666666669e-01,
91: 1.1111111111111112e-01, 8.3333333333333343e-02, 0.0000000000000000e+00, -1.9444444444444448e-01, -1.6666666666666669e-01, 0.0000000000000000e+00, 1.1111111111111112e-01, 0.0000000000000000e+00, 8.3333333333333343e-02,
92: -2.2222222222222227e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 5.5555555555555558e-01, 1.6666666666666669e-01, 1.6666666666666669e-01, -8.3333333333333356e-02, -1.3888888888888887e-01, -8.3333333333333356e-02,
93: 0.0000000000000000e+00, -1.9444444444444448e-01, -1.6666666666666666e-01, 0.0000000000000000e+00, -2.7777777777777769e-02, 0.0000000000000000e+00, 8.3333333333333343e-02, 1.1111111111111112e-01, 0.0000000000000000e+00,
94: -1.6666666666666669e-01, -1.9444444444444448e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 1.1111111111111112e-01, 8.3333333333333343e-02,
95: 1.6666666666666669e-01, 5.5555555555555558e-01, 1.6666666666666669e-01, -8.3333333333333356e-02, -8.3333333333333356e-02, -1.3888888888888887e-01, 0.0000000000000000e+00, -1.6666666666666666e-01, -1.9444444444444448e-01,
96: -1.6666666666666669e-01, 0.0000000000000000e+00, -1.9444444444444448e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.2222222222222227e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, -2.7777777777777769e-02,
97: 8.3333333333333343e-02, 0.0000000000000000e+00, 1.1111111111111113e-01, 0.0000000000000000e+00, 8.3333333333333343e-02, 1.1111111111111113e-01, 1.6666666666666669e-01, 1.6666666666666669e-01, 5.5555555555555558e-01};
99: typedef enum {
100: PDE_POISSON,
101: PDE_ELASTICITY
102: } PDEType;
104: typedef struct {
105: PDEType pde;
106: PetscInt dim;
107: PetscInt dof;
108: PetscInt cells[3];
109: PetscBool useglobal;
110: PetscBool multi_element;
111: PetscBool dirbc;
112: PetscBool per[3];
113: PetscBool test;
114: PetscScalar *elemMat;
115: PetscBool use_composite_pc;
116: PetscBool test_reuse;
117: PetscBool random_initial_guess;
118: PetscBool random_real;
119: } AppCtx;
121: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
122: {
123: const char *pdeTypes[2] = {"Poisson", "Elasticity"};
124: PetscInt n, pde;
126: PetscFunctionBeginUser;
127: options->pde = PDE_POISSON;
128: options->elemMat = NULL;
129: options->dim = 1;
130: options->cells[0] = 8;
131: options->cells[1] = 6;
132: options->cells[2] = 4;
133: options->useglobal = PETSC_FALSE;
134: options->multi_element = PETSC_FALSE;
135: options->dirbc = PETSC_TRUE;
136: options->test = PETSC_FALSE;
137: options->per[0] = PETSC_FALSE;
138: options->per[1] = PETSC_FALSE;
139: options->per[2] = PETSC_FALSE;
140: options->use_composite_pc = PETSC_FALSE;
141: options->test_reuse = PETSC_FALSE;
142: options->random_initial_guess = PETSC_FALSE;
143: options->random_real = PETSC_FALSE;
145: PetscOptionsBegin(comm, NULL, "Problem Options", NULL);
146: pde = options->pde;
147: PetscCall(PetscOptionsEList("-pde_type", "The PDE type", __FILE__, pdeTypes, 2, pdeTypes[options->pde], &pde, NULL));
148: options->pde = (PDEType)pde;
149: PetscCall(PetscOptionsInt("-dim", "The topological mesh dimension", __FILE__, options->dim, &options->dim, NULL));
150: PetscCall(PetscOptionsIntArray("-cells", "The mesh division", __FILE__, options->cells, (n = 3, &n), NULL));
151: PetscCall(PetscOptionsBoolArray("-periodicity", "The mesh periodicity", __FILE__, options->per, (n = 3, &n), NULL));
152: PetscCall(PetscOptionsBool("-use_global", "Test MatSetValues", __FILE__, options->useglobal, &options->useglobal, NULL));
153: PetscCall(PetscOptionsBool("-multi_element", "Use multi-element BDDC", __FILE__, options->multi_element, &options->multi_element, NULL));
154: PetscCall(PetscOptionsBool("-dirichlet", "Use dirichlet BC", __FILE__, options->dirbc, &options->dirbc, NULL));
155: PetscCall(PetscOptionsBool("-use_composite_pc", "Multiplicative composite with BDDC + Richardson/Jacobi", __FILE__, options->use_composite_pc, &options->use_composite_pc, NULL));
156: PetscCall(PetscOptionsBool("-test_reuse", "Set up the preconditioner again with new matrix values", __FILE__, options->test_reuse, &options->test_reuse, NULL));
157: PetscCall(PetscOptionsBool("-random_initial_guess", "Solve A x = 0 with random initial guess, instead of A x = b with random b", __FILE__, options->random_initial_guess, &options->random_initial_guess, NULL));
158: PetscCall(PetscOptionsBool("-random_real", "Use real-valued b (or x, if -random_initial_guess) instead of default scalar type", __FILE__, options->random_real, &options->random_real, NULL));
159: PetscOptionsEnd();
161: for (n = options->dim; n < 3; n++) options->cells[n] = 0;
162: if (options->per[0]) options->dirbc = PETSC_FALSE;
164: /* element matrices */
165: switch (options->pde) {
166: case PDE_ELASTICITY:
167: options->dof = options->dim;
168: switch (options->dim) {
169: case 1:
170: options->elemMat = elast_1D_emat;
171: break;
172: case 2:
173: options->elemMat = elast_2D_emat;
174: break;
175: case 3:
176: options->elemMat = elast_3D_emat;
177: break;
178: default:
179: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, options->dim);
180: }
181: break;
182: case PDE_POISSON:
183: options->dof = 1;
184: switch (options->dim) {
185: case 1:
186: options->elemMat = poiss_1D_emat;
187: break;
188: case 2:
189: options->elemMat = poiss_2D_emat;
190: break;
191: case 3:
192: options->elemMat = poiss_3D_emat;
193: break;
194: default:
195: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, options->dim);
196: }
197: break;
198: default:
199: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported PDE %d", options->pde);
200: }
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: int main(int argc, char **args)
205: {
206: AppCtx user;
207: KSP ksp;
208: PC pc;
209: Mat A;
210: DM da;
211: Vec x, b, xcoor, xcoorl;
212: IS zero, test_primal = NULL;
213: ISLocalToGlobalMapping map;
214: MatNullSpace nullsp = NULL;
215: PetscInt i;
216: PetscInt nel, nen; /* Number of elements & element nodes */
217: const PetscInt *e_loc; /* Local indices of element nodes (in local element order) */
218: PetscInt *e_glo = NULL; /* Global indices of element nodes (in local element order) */
219: PetscInt nodes[3], test_primal_vertex = -1;
220: PetscBool ismatis, flg;
221: PetscLogStage stages[2];
223: PetscFunctionBeginUser;
224: PetscCall(PetscInitialize(&argc, &args, NULL, help));
225: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
226: for (i = 0; i < 3; i++) nodes[i] = user.cells[i] + !user.per[i];
227: switch (user.dim) {
228: case 3:
229: PetscCall(DMDACreate3d(PETSC_COMM_WORLD, user.per[0] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE, user.per[1] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE, user.per[2] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nodes[0], nodes[1], nodes[2], PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
230: user.dof, 1, NULL, NULL, NULL, &da));
231: break;
232: case 2:
233: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, user.per[0] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE, user.per[1] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nodes[0], nodes[1], PETSC_DECIDE, PETSC_DECIDE, user.dof, 1, NULL, NULL, &da));
234: break;
235: case 1:
236: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, user.per[0] ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE, nodes[0], user.dof, 1, NULL, &da));
237: break;
238: default:
239: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, user.dim);
240: }
242: PetscCall(PetscLogStageRegister("KSPSetUp", &stages[0]));
243: PetscCall(PetscLogStageRegister("KSPSolve", &stages[1]));
245: PetscCall(DMSetMatType(da, MATIS));
246: PetscCall(DMSetFromOptions(da));
247: PetscCall(DMDASetElementType(da, DMDA_ELEMENT_Q1));
248: PetscCall(DMSetUp(da));
249: {
250: PetscInt M, N, P;
251: PetscCall(DMDAGetInfo(da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0));
252: switch (user.dim) {
253: case 3:
254: user.cells[2] = P - !user.per[2]; /* fall through */
255: case 2:
256: user.cells[1] = N - !user.per[1]; /* fall through */
257: case 1:
258: user.cells[0] = M - !user.per[0];
259: break;
260: default:
261: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT, user.dim);
262: }
263: }
264: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0 * user.cells[0], 0.0, 1.0 * user.cells[1], 0.0, 1.0 * user.cells[2]));
265: PetscCall(DMGetCoordinates(da, &xcoor));
267: PetscCall(DMCreateMatrix(da, &A));
268: PetscCall(MatSetFromOptions(A));
269: PetscCall(DMGetLocalToGlobalMapping(da, &map));
270: PetscCall(DMDAGetElements(da, &nel, &nen, &e_loc));
271: if (user.useglobal) {
272: PetscCall(PetscMalloc1(nel * nen, &e_glo));
273: PetscCall(ISLocalToGlobalMappingApplyBlock(map, nen * nel, e_loc, e_glo));
274: }
276: if (user.multi_element) {
277: ISLocalToGlobalMapping mapn;
278: PetscInt *el_glo = NULL, m, n, M, N, *el_sizes;
279: Mat lA;
281: PetscCall(PetscMalloc1(nel * nen, &el_glo));
282: PetscCall(ISLocalToGlobalMappingApplyBlock(map, nen * nel, e_loc, el_glo));
283: PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)map), user.dof, nen * nel, el_glo, PETSC_OWN_POINTER, &mapn));
284: PetscCall(MatGetLocalSize(A, &m, &n));
285: PetscCall(MatGetSize(A, &M, &N));
286: PetscCall(MatDestroy(&A));
287: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
288: PetscCall(MatSetSizes(A, m, n, M, N));
289: PetscCall(MatSetBlockSize(A, user.dof));
290: PetscCall(MatSetType(A, MATIS));
291: PetscCall(MatISSetAllowRepeated(A, PETSC_TRUE));
292: PetscCall(MatSetLocalToGlobalMapping(A, mapn, mapn));
293: PetscCall(MatISSetPreallocation(A, user.dof * nen, NULL, user.dof * nen, NULL));
294: PetscCall(ISLocalToGlobalMappingViewFromOptions(mapn, NULL, "-multi_view"));
295: PetscCall(ISLocalToGlobalMappingDestroy(&mapn));
297: /* The information set with MatSetVariableBlockSizes on the local mat
298: can be used to detect the local elements instead of having to analyze
299: the sparsity pattern of the local matrix */
300: PetscCall(MatISGetLocalMat(A, &lA));
301: PetscCall(PetscMalloc1(nel, &el_sizes));
302: for (i = 0; i < nel; i++) el_sizes[i] = user.dof * nen;
303: PetscCall(MatSetVariableBlockSizes(lA, nel, el_sizes));
304: PetscCall(PetscFree(el_sizes));
305: }
307: /* we reorder the indices since the element matrices are given in lexicographic order,
308: whereas the elements indices returned by DMDAGetElements follow the usual FEM ordering
309: i.e., element matrices DMDA ordering
310: 2---3 3---2
311: / / / /
312: 0---1 0---1
313: */
314: for (i = 0; i < nel; ++i) {
315: PetscInt ord[8] = {0, 1, 3, 2, 4, 5, 7, 6};
316: PetscInt j, idxs[8];
318: PetscCheck(nen <= 8, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not coded");
319: if (!user.useglobal) {
320: if (user.multi_element) {
321: for (j = 0; j < nen; j++) idxs[j] = i * nen + ord[j];
322: } else {
323: for (j = 0; j < nen; j++) idxs[j] = e_loc[i * nen + ord[j]];
324: }
325: PetscCall(MatSetValuesBlockedLocal(A, nen, idxs, nen, idxs, user.elemMat, ADD_VALUES));
326: } else {
327: for (j = 0; j < nen; j++) idxs[j] = e_glo[i * nen + ord[j]];
328: PetscCall(MatSetValuesBlocked(A, nen, idxs, nen, idxs, user.elemMat, ADD_VALUES));
329: }
330: }
331: PetscCall(DMDARestoreElements(da, &nel, &nen, &e_loc));
332: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
333: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
334: PetscCall(MatViewFromOptions(A, NULL, "-A_mat_view"));
335: PetscCall(MatSetOption(A, MAT_SPD, PETSC_TRUE));
336: PetscCall(MatSetOption(A, MAT_SPD_ETERNAL, PETSC_TRUE));
338: /* Boundary conditions */
339: zero = NULL;
340: if (user.dirbc) { /* fix one side of DMDA */
341: Vec nat, glob;
342: PetscScalar *vals;
343: PetscInt n, *idx, j, st;
345: n = PetscGlobalRank ? 0 : (user.cells[1] + 1) * (user.cells[2] + 1);
346: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, user.cells[0] + 1, &zero));
347: if (user.dof > 1) { /* zero all components */
348: const PetscInt *idx;
349: IS bzero;
351: PetscCall(ISGetIndices(zero, &idx));
352: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, user.dof, n, idx, PETSC_COPY_VALUES, &bzero));
353: PetscCall(ISRestoreIndices(zero, &idx));
354: PetscCall(ISDestroy(&zero));
355: zero = bzero;
356: }
357: /* map indices from natural to global */
358: PetscCall(DMDACreateNaturalVector(da, &nat));
359: PetscCall(ISGetLocalSize(zero, &n));
360: PetscCall(PetscMalloc1(n, &vals));
361: for (i = 0; i < n; i++) vals[i] = 1.0;
362: PetscCall(ISGetIndices(zero, (const PetscInt **)&idx));
363: PetscCall(VecSetValues(nat, n, idx, vals, INSERT_VALUES));
364: PetscCall(ISRestoreIndices(zero, (const PetscInt **)&idx));
365: PetscCall(PetscFree(vals));
366: PetscCall(VecAssemblyBegin(nat));
367: PetscCall(VecAssemblyEnd(nat));
368: PetscCall(DMCreateGlobalVector(da, &glob));
369: PetscCall(DMDANaturalToGlobalBegin(da, nat, INSERT_VALUES, glob));
370: PetscCall(DMDANaturalToGlobalEnd(da, nat, INSERT_VALUES, glob));
371: PetscCall(VecDestroy(&nat));
372: PetscCall(ISDestroy(&zero));
373: PetscCall(VecGetLocalSize(glob, &n));
374: PetscCall(PetscMalloc1(n, &idx));
375: PetscCall(VecGetOwnershipRange(glob, &st, NULL));
376: PetscCall(VecGetArray(glob, &vals));
377: for (i = 0, j = 0; i < n; i++)
378: if (PetscRealPart(vals[i]) == 1.0) idx[j++] = i + st;
379: PetscCall(VecRestoreArray(glob, &vals));
380: PetscCall(VecDestroy(&glob));
381: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, j, idx, PETSC_OWN_POINTER, &zero));
382: PetscCall(MatZeroRowsColumnsIS(A, zero, 1.0, NULL, NULL));
383: } else {
384: switch (user.pde) {
385: case PDE_POISSON:
386: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nullsp));
387: break;
388: case PDE_ELASTICITY:
389: PetscCall(MatNullSpaceCreateRigidBody(xcoor, &nullsp));
390: break;
391: }
392: /* with periodic BC and Elasticity, just the displacements are in the nullspace
393: this is no harm since we eliminate all the components of the rhs */
394: PetscCall(MatSetNullSpace(A, nullsp));
395: }
397: PetscCall(PetscOptionsHasName(NULL, NULL, "-assembled_view", &flg));
398: if (flg) {
399: Mat AA;
401: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &AA));
402: PetscCall(MatViewFromOptions(AA, NULL, "-assembled_view"));
403: PetscCall(MatDestroy(&AA));
404: }
406: /* Attach near null space for elasticity */
407: if (user.pde == PDE_ELASTICITY) {
408: MatNullSpace nearnullsp;
410: PetscCall(MatNullSpaceCreateRigidBody(xcoor, &nearnullsp));
411: PetscCall(MatSetNearNullSpace(A, nearnullsp));
412: PetscCall(MatNullSpaceDestroy(&nearnullsp));
413: }
415: /* we may want to use MG for the local solvers: attach local nearnullspace to the local matrices */
416: PetscCall(DMGetCoordinatesLocal(da, &xcoorl));
417: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
418: if (ismatis) {
419: MatNullSpace lnullsp = NULL;
420: Mat lA;
422: PetscCall(MatISGetLocalMat(A, &lA));
423: if (user.pde == PDE_ELASTICITY) {
424: Vec lc;
425: ISLocalToGlobalMapping l2l;
426: IS is;
427: const PetscScalar *a;
428: const PetscInt *idxs;
429: PetscInt n, bs;
431: /* when using a DMDA, the local matrices have an additional local-to-local map
432: that maps from the DA local ordering to the ordering induced by the elements */
433: PetscCall(MatGetLocalToGlobalMapping(lA, &l2l, NULL));
434: if (l2l) {
435: PetscCall(MatCreateVecs(lA, &lc, NULL));
436: PetscCall(VecSetLocalToGlobalMapping(lc, l2l));
438: PetscCall(VecSetOption(lc, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE));
439: PetscCall(VecGetLocalSize(xcoorl, &n));
440: PetscCall(VecGetBlockSize(xcoorl, &bs));
441: PetscCall(ISCreateStride(PETSC_COMM_SELF, n / bs, 0, 1, &is));
442: PetscCall(ISGetIndices(is, &idxs));
443: PetscCall(VecGetArrayRead(xcoorl, &a));
444: PetscCall(VecSetValuesBlockedLocal(lc, n / bs, idxs, a, INSERT_VALUES));
445: PetscCall(VecAssemblyBegin(lc));
446: PetscCall(VecAssemblyEnd(lc));
447: PetscCall(VecRestoreArrayRead(xcoorl, &a));
448: PetscCall(ISRestoreIndices(is, &idxs));
449: PetscCall(ISDestroy(&is));
450: PetscCall(MatNullSpaceCreateRigidBody(lc, &lnullsp));
451: PetscCall(VecDestroy(&lc));
452: }
453: } else if (user.pde == PDE_POISSON) {
454: PetscCall(MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &lnullsp));
455: }
456: PetscCall(MatSetNearNullSpace(lA, lnullsp));
457: PetscCall(MatNullSpaceDestroy(&lnullsp));
458: PetscCall(MatISRestoreLocalMat(A, &lA));
459: }
461: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
462: PetscCall(KSPSetOperators(ksp, A, A));
463: PetscCall(KSPSetType(ksp, KSPCG));
464: PetscCall(KSPGetPC(ksp, &pc));
465: if (ismatis) {
466: if (user.use_composite_pc) {
467: PC pcksp, pcjacobi;
468: KSP ksprich;
469: PetscCall(PCSetType(pc, PCCOMPOSITE));
470: PetscCall(PCCompositeSetType(pc, PC_COMPOSITE_MULTIPLICATIVE));
471: PetscCall(PCCompositeAddPCType(pc, PCBDDC));
472: PetscCall(PCCompositeAddPCType(pc, PCKSP));
473: PetscCall(PCCompositeGetPC(pc, 1, &pcksp));
474: PetscCall(PCKSPGetKSP(pcksp, &ksprich));
475: PetscCall(KSPSetType(ksprich, KSPRICHARDSON));
476: PetscCall(KSPSetTolerances(ksprich, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 1));
477: PetscCall(KSPSetNormType(ksprich, KSP_NORM_NONE));
478: PetscCall(KSPSetConvergenceTest(ksprich, KSPConvergedSkip, NULL, NULL));
479: PetscCall(KSPGetPC(ksprich, &pcjacobi));
480: PetscCall(PCSetType(pcjacobi, PCJACOBI));
481: } else {
482: PetscCall(PCSetType(pc, PCBDDC));
483: }
484: }
485: PetscCall(KSPSetFromOptions(ksp));
487: /* test user-defined primal vertices API */
488: PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_primal_vertex", &test_primal_vertex, NULL));
489: if (test_primal_vertex >= 0) {
490: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, PetscGlobalRank ? 0 : 1, &test_primal_vertex, PETSC_COPY_VALUES, &test_primal));
491: PetscCall(PCBDDCSetPrimalVerticesIS(pc, test_primal));
492: }
494: PetscCall(PetscLogStagePush(stages[0]));
495: PetscCall(KSPSetUp(ksp));
496: PetscCall(PetscLogStagePop());
497: if (user.test_reuse) {
498: PetscCall(MatScale(A, 2.0));
499: PetscCall(KSPSetUp(ksp));
500: }
502: /* test that user-defined primal vertices are respected */
503: if (test_primal) {
504: IS stored, corners;
505: Mat lA;
506: ISLocalToGlobalMapping l2g, l2l;
507: const PetscInt *idx;
508: PetscInt n, localvertex, pos;
510: PetscCall(PCBDDCGetPrimalVerticesIS(pc, &stored));
511: PetscCheck(stored, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Global user primal vertices were lost during setup");
512: PetscCall(ISEqual(test_primal, stored, &flg));
513: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Global user primal vertices changed during setup");
514: PetscCall(ISDestroy(&test_primal));
515: PetscCall(PCBDDCGetPrimalVerticesLocalIS(pc, &stored));
516: PetscCheck(stored, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Missing local primal vertices");
517: PetscCall(MatGetLocalToGlobalMapping(A, &l2g, NULL));
518: PetscCall(ISGlobalToLocalMappingApply(l2g, IS_GTOLM_MASK, 1, &test_primal_vertex, NULL, &localvertex));
519: if (localvertex >= 0) {
520: PetscCall(ISLocate(stored, localvertex, &pos));
521: PetscCheck(pos >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing user primal vertex in local numbering");
522: }
523: /* Automatic DMDA corners must be retained alongside the user vertex. */
524: PetscCall(DMDAGetSubdomainCornersIS(da, &corners));
525: PetscCall(MatISGetLocalMat(A, &lA));
526: PetscCall(MatGetLocalToGlobalMapping(lA, &l2l, NULL));
527: PetscCheck(corners && l2l, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing DMDA corner data");
528: PetscCall(ISGetLocalSize(corners, &n));
529: PetscCall(ISGetIndices(corners, &idx));
530: for (PetscInt c = 0; c < n; c++) {
531: for (PetscInt d = 0; d < user.dof; d++) {
532: PetscInt corner = user.dof * idx[c] + d;
534: PetscCall(ISLocalToGlobalMappingApply(l2l, 1, &corner, &localvertex));
535: PetscCall(ISLocate(stored, localvertex, &pos));
536: PetscCheck(pos >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing automatic DMDA corner in local primal vertices");
537: }
538: }
539: PetscCall(ISRestoreIndices(corners, &idx));
540: PetscCall(MatISRestoreLocalMat(A, &lA));
541: PetscCall(DMDARestoreSubdomainCornersIS(da, &corners));
542: }
544: PetscCall(MatCreateVecs(A, &x, &b));
545: if (user.random_initial_guess) {
546: /* Solving A x = 0 with random initial guess allows Arnoldi to run for more iterations, thereby yielding a more
547: * complete Hessenberg matrix and more accurate eigenvalues. */
548: PetscCall(VecSetRandom(x, NULL));
549: if (user.random_real) PetscCall(VecRealPart(x));
550: if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, x));
551: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
552: PetscCall(KSPSetComputeEigenvalues(ksp, PETSC_TRUE));
553: PetscCall(KSPGMRESSetRestart(ksp, 100));
554: } else {
555: PetscCall(VecSetRandom(b, NULL));
556: if (user.random_real) PetscCall(VecRealPart(x));
557: if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, b));
558: }
559: PetscCall(PetscLogStagePush(stages[1]));
560: PetscCall(KSPSolve(ksp, b, x));
561: PetscCall(PetscLogStagePop());
563: /* cleanup */
564: PetscCall(VecDestroy(&x));
565: PetscCall(VecDestroy(&b));
566: PetscCall(ISDestroy(&zero));
567: PetscCall(PetscFree(e_glo));
568: PetscCall(MatNullSpaceDestroy(&nullsp));
569: PetscCall(KSPDestroy(&ksp));
570: PetscCall(MatDestroy(&A));
571: PetscCall(DMDestroy(&da));
572: PetscCall(PetscFinalize());
573: return 0;
574: }
576: /*TEST
578: test:
579: suffix: bddc_global_primal
580: nsize: 4
581: output_file: output/empty.out
582: args: -dim 2 -cells 4,4 -pde_type Poisson -da_processors_x 2 -da_processors_y 2 -test_primal_vertex 5 -test_reuse -ksp_error_if_not_converged -pc_bddc_coarse_redundant_pc_type svd
584: test:
585: nsize: 8
586: filter: grep -v "variant HERMITIAN"
587: suffix: bddc_1
588: args: -pde_type Poisson -dim 3 -dirichlet 0 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged
589: test:
590: nsize: 8
591: filter: grep -v "variant HERMITIAN"
592: suffix: bddc_2
593: args: -pde_type Poisson -dim 3 -dirichlet 0 -ksp_view -use_global -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged
594: test:
595: nsize: 8
596: filter: grep -v "variant HERMITIAN"
597: suffix: bddc_elast
598: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic
599: test:
600: nsize: 8
601: filter: grep -v "variant HERMITIAN"
602: suffix: bddc_elast_3lev
603: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_levels 1 -pc_bddc_coarsening_ratio 1 -ksp_error_if_not_converged -pc_bddc_monolithic -pc_bddc_use_faces -pc_bddc_coarse_pc_bddc_corner_selection
604: testset:
605: nsize: 8
606: requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
607: # on some architectures, this test will converge in slightly different numbers of iterations
608: filter: grep -v "variant HERMITIAN" | grep -v " tolerance" | sed -E "s/CONVERGED_RTOL iterations (1[8-9]|2[0-1])$/CONVERGED_RTOL iterations 20/g"
609: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_levels 1 -pc_bddc_coarsening_ratio 1 -ksp_error_if_not_converged -pc_bddc_monolithic -pc_bddc_use_faces -pc_bddc_coarse_pc_type hpddm -prefix_push pc_bddc_coarse_ -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 6 -pc_hpddm_levels_1_st_pc_factor_shift_type INBLOCKS -prefix_pop -ksp_type fgmres -ksp_max_it 50 -ksp_converged_reason
610: test:
611: args: -pc_bddc_coarse_pc_hpddm_coarse_mat_type baij -options_left no
612: suffix: bddc_elast_3lev_hpddm_baij
613: test:
614: requires: !complex
615: suffix: bddc_elast_3lev_hpddm
616: test:
617: nsize: 8
618: requires: !single
619: filter: grep -v "variant HERMITIAN"
620: suffix: bddc_elast_4lev
621: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_levels 2 -pc_bddc_coarsening_ratio 2 -ksp_error_if_not_converged -pc_bddc_monolithic -pc_bddc_use_faces -pc_bddc_coarse_pc_bddc_corner_selection -pc_bddc_coarse_l1_pc_bddc_corner_selection -pc_bddc_aggregator_mat_partitioning_type average -options_left 0 -pc_bddc_aggregator_mat_partitioning_view
622: test:
623: nsize: 8
624: filter: grep -v "variant HERMITIAN"
625: suffix: bddc_elast_deluxe_layers
626: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -pc_bddc_use_deluxe_scaling -pc_bddc_schur_layers 1
627: test:
628: nsize: 8
629: filter: grep -v "variant HERMITIAN" | sed -e "s/iterations 1[0-9]/iterations 10/g"
630: suffix: bddc_elast_dir_approx
631: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -pc_bddc_dirichlet_pc_type gamg -pc_bddc_dirichlet_pc_gamg_esteig_ksp_max_it 10 -ksp_converged_reason -pc_bddc_dirichlet_approximate
632: test:
633: nsize: 8
634: filter: grep -v "variant HERMITIAN" | sed -e "s/iterations 1[0-9]/iterations 10/g"
635: suffix: bddc_elast_neu_approx
636: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -pc_bddc_neumann_pc_type gamg -pc_bddc_neumann_pc_gamg_esteig_ksp_max_it 10 -ksp_converged_reason -pc_bddc_neumann_approximate
637: test:
638: nsize: 8
639: filter: grep -v "variant HERMITIAN" | sed -e "s/iterations 1[0-9]/iterations 10/g"
640: suffix: bddc_elast_both_approx
641: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -pc_bddc_dirichlet_pc_type gamg -pc_bddc_dirichlet_pc_gamg_esteig_ksp_max_it 10 -pc_bddc_neumann_pc_type gamg -pc_bddc_neumann_pc_gamg_esteig_ksp_max_it 10 -ksp_converged_reason -pc_bddc_neumann_approximate -pc_bddc_dirichlet_approximate
642: test:
643: nsize: 8
644: filter: grep -v "variant HERMITIAN"
645: suffix: fetidp_1
646: args: -pde_type Poisson -dim 3 -dirichlet 0 -ksp_view -ksp_type fetidp -fetidp_ksp_type cg -fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd -ksp_fetidp_fullyredundant -ksp_error_if_not_converged
647: test:
648: nsize: 8
649: filter: grep -v "variant HERMITIAN"
650: suffix: fetidp_2
651: args: -pde_type Poisson -dim 3 -dirichlet 0 -ksp_view -use_global -ksp_type fetidp -fetidp_ksp_type cg -fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd -ksp_fetidp_fullyredundant -ksp_error_if_not_converged
652: test:
653: nsize: 8
654: filter: grep -v "variant HERMITIAN"
655: suffix: fetidp_elast
656: args: -pde_type Elasticity -cells 9,7,8 -dim 3 -ksp_view -ksp_type fetidp -fetidp_ksp_type cg -fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd -ksp_fetidp_fullyredundant -ksp_error_if_not_converged -fetidp_bddc_pc_bddc_monolithic
657: testset:
658: nsize: 8
659: requires: hpddm slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
660: args: -pde_type Elasticity -cells 12,12 -dim 2 -ksp_converged_reason -pc_type hpddm -pc_hpddm_coarse_correction balanced -pc_hpddm_levels_1_pc_type asm -pc_hpddm_levels_1_pc_asm_overlap 1 -pc_hpddm_levels_1_pc_asm_type basic -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 10 -pc_hpddm_levels_1_st_pc_factor_shift_type INBLOCKS
661: test:
662: args: -mat_is_localmat_type {{aij baij sbaij}shared output} -pc_hpddm_coarse_mat_type {{baij sbaij}shared output}
663: suffix: hpddm
664: output_file: output/ex71_hpddm.out
665: filter: sed -e "s/CONVERGED_RTOL iterations 15/CONVERGED_RTOL iterations 14/g"
666: test:
667: args: -mat_is_localmat_type sbaij -pc_hpddm_coarse_mat_type sbaij -pc_hpddm_levels_1_st_share_sub_ksp -pc_hpddm_levels_1_eps_type lapack -pc_hpddm_levels_1_eps_smallest_magnitude -pc_hpddm_levels_1_st_type shift
668: suffix: hpddm_lapack
669: output_file: output/ex71_hpddm.out
670: testset:
671: nsize: 9
672: args: -assembled_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged
673: test:
674: args: -dim 1 -cells 12 -pde_type Poisson
675: suffix: dmda_matis_poiss_1d_loc
676: output_file: output/ex71_dmda_matis_poiss_1d.out
677: test:
678: args: -dim 1 -cells 12 -pde_type Poisson -use_global
679: suffix: dmda_matis_poiss_1d_glob
680: output_file: output/ex71_dmda_matis_poiss_1d.out
681: test:
682: args: -dim 1 -cells 12 -pde_type Elasticity
683: suffix: dmda_matis_elast_1d_loc
684: output_file: output/ex71_dmda_matis_elast_1d.out
685: test:
686: args: -dim 1 -cells 12 -pde_type Elasticity -use_global
687: suffix: dmda_matis_elast_1d_glob
688: output_file: output/ex71_dmda_matis_elast_1d.out
689: test:
690: args: -dim 2 -cells 5,7 -pde_type Poisson
691: suffix: dmda_matis_poiss_2d_loc
692: output_file: output/ex71_dmda_matis_poiss_2d.out
693: test:
694: args: -dim 2 -cells 5,7 -pde_type Poisson -use_global
695: suffix: dmda_matis_poiss_2d_glob
696: output_file: output/ex71_dmda_matis_poiss_2d.out
697: test:
698: args: -dim 2 -cells 5,7 -pde_type Elasticity
699: suffix: dmda_matis_elast_2d_loc
700: output_file: output/ex71_dmda_matis_elast_2d.out
701: test:
702: args: -dim 2 -cells 5,7 -pde_type Elasticity -use_global
703: suffix: dmda_matis_elast_2d_glob
704: output_file: output/ex71_dmda_matis_elast_2d.out
705: test:
706: args: -dim 3 -cells 3,3,3 -pde_type Poisson
707: suffix: dmda_matis_poiss_3d_loc
708: output_file: output/ex71_dmda_matis_poiss_3d.out
709: test:
710: args: -dim 3 -cells 3,3,3 -pde_type Poisson -use_global
711: suffix: dmda_matis_poiss_3d_glob
712: output_file: output/ex71_dmda_matis_poiss_3d.out
713: test:
714: args: -dim 3 -cells 3,3,3 -pde_type Elasticity
715: suffix: dmda_matis_elast_3d_loc
716: output_file: output/ex71_dmda_matis_elast_3d.out
717: test:
718: args: -dim 3 -cells 3,3,3 -pde_type Elasticity -use_global
719: suffix: dmda_matis_elast_3d_glob
720: output_file: output/ex71_dmda_matis_elast_3d.out
721: test:
722: nsize: 8
723: filter: grep -v "variant HERMITIAN" | sed -e "s/CONVERGED_RTOL iterations 1[0-9]/CONVERGED_RTOL iterations 13/g"
724: suffix: bddc_elast_deluxe_layers_adapt
725: requires: mumps !complex
726: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_converged_reason -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -sub_schurs_mat_solver_type mumps -pc_bddc_use_deluxe_scaling -pc_bddc_adaptive_threshold 2.0 -pc_bddc_schur_layers {{1 10}separate_output} -pc_bddc_adaptive_userdefined {{0 1}separate output} -sub_schurs_schur_mat_type seqdense
727: # gitlab runners have a quite old MKL (2016) which interacts badly with AMD machines (not Intel-based ones!)
728: # this is the reason behind the filtering rule
729: test:
730: nsize: 8
731: suffix: bddc_elast_deluxe_layers_adapt_mkl_pardiso
732: filter: sed -e "s/CONVERGED_RTOL iterations [1-2][0-9]/CONVERGED_RTOL iterations 13/g" | sed -e "s/CONVERGED_RTOL iterations 6/CONVERGED_RTOL iterations 5/g"
733: requires: mkl_pardiso !complex
734: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_converged_reason -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -sub_schurs_mat_solver_type mkl_pardiso -sub_schurs_mat_mkl_pardiso_65 1 -pc_bddc_use_deluxe_scaling -pc_bddc_adaptive_threshold 2.0 -pc_bddc_schur_layers {{1 10}separate_output} -pc_bddc_adaptive_userdefined {{0 1}separate output} -sub_schurs_schur_mat_type seqdense
735: test:
736: nsize: 8
737: filter: grep -v "variant HERMITIAN"
738: suffix: bddc_cusparse
739: # no kokkos since it seems kokkos's resource demand is too much with 8 ranks and the test will fail on cuda related initialization.
740: requires: cuda !kokkos
741: args: -pde_type Poisson -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_dirichlet_pc_type cholesky -pc_bddc_dirichlet_pc_factor_mat_solver_type cusparse -pc_bddc_dirichlet_pc_factor_mat_ordering_type nd -pc_bddc_neumann_pc_type cholesky -pc_bddc_neumann_pc_factor_mat_solver_type cusparse -pc_bddc_neumann_pc_factor_mat_ordering_type nd -mat_is_localmat_type aijcusparse
742: test:
743: nsize: 8
744: filter: grep -v "variant HERMITIAN"
745: suffix: bddc_elast_deluxe_layers_adapt_cuda
746: requires: !complex mumps cuda defined(PETSC_HAVE_CUSOLVERDNDPOTRI)
747: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_converged_reason -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -sub_schurs_mat_solver_type mumps -pc_bddc_use_deluxe_scaling -pc_bddc_adaptive_threshold 2.0 -pc_bddc_schur_layers {{1 10}separate_output} -pc_bddc_adaptive_userdefined {{0 1}separate output} -mat_is_localmat_type seqaijcusparse -sub_schurs_schur_mat_type {{seqdensecuda seqdense}}
748: test:
749: nsize: 8
750: filter: grep -v "variant HERMITIAN" | grep -v "I-node routines" | sed -e "s/seqaijcusparse/seqaij/g"
751: suffix: bddc_elast_deluxe_layers_adapt_cuda_approx
752: requires: !complex mumps cuda defined(PETSC_HAVE_CUSOLVERDNDPOTRI)
753: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_view -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -sub_schurs_mat_solver_type mumps -pc_bddc_use_deluxe_scaling -pc_bddc_adaptive_threshold 2.0 -pc_bddc_schur_layers 1 -mat_is_localmat_type {{seqaij seqaijcusparse}separate output} -sub_schurs_schur_mat_type {{seqdensecuda seqdense}} -pc_bddc_dirichlet_pc_type gamg -pc_bddc_dirichlet_approximate -pc_bddc_neumann_pc_type gamg -pc_bddc_neumann_approximate -pc_bddc_dirichlet_pc_gamg_esteig_ksp_max_it 10 -pc_bddc_neumann_pc_gamg_esteig_ksp_max_it 10
754: test:
755: nsize: 8
756: suffix: bddc_elast_deluxe_layers_adapt_mkl_pardiso_cuda
757: requires: !complex mkl_pardiso cuda defined(PETSC_HAVE_CUSOLVERDNDPOTRI)
758: filter: sed -e "s/CONVERGED_RTOL iterations 6/CONVERGED_RTOL iterations 5/g"
759: args: -pde_type Elasticity -cells 7,9,8 -dim 3 -ksp_converged_reason -pc_bddc_coarse_redundant_pc_type svd -ksp_error_if_not_converged -pc_bddc_monolithic -sub_schurs_mat_solver_type mkl_pardiso -sub_schurs_mat_mkl_pardiso_65 1 -pc_bddc_use_deluxe_scaling -pc_bddc_adaptive_threshold 2.0 -pc_bddc_schur_layers {{1 10}separate_output} -pc_bddc_adaptive_userdefined {{0 1}separate output} -mat_is_localmat_type seqaijcusparse -sub_schurs_schur_mat_type {{seqdensecuda seqdense}}
761: testset:
762: nsize: 2
763: output_file: output/ex71_aij_dmda_preall.out
764: filter: sed -e "s/CONVERGED_RTOL iterations 7/CONVERGED_RTOL iterations 6/g"
765: args: -pde_type Poisson -dim 1 -cells 6 -pc_type none -ksp_converged_reason
766: test:
767: suffix: aijviennacl_dmda_preall
768: requires: viennacl
769: args: -dm_mat_type aijviennacl -dm_preallocate_only {{0 1}} -dirichlet {{0 1}}
770: # -dm_preallocate_only 0 is broken
771: test:
772: suffix: aijcusparse_dmda_preall
773: requires: cuda
774: args: -dm_mat_type aijcusparse -dm_preallocate_only -dirichlet {{0 1}}
775: test:
776: suffix: aij_dmda_preall
777: args: -dm_mat_type aij -dm_preallocate_only {{0 1}} -dirichlet {{0 1}}
778: testset:
779: nsize: 4
780: args: -dim 2 -cells 16,16 -periodicity 1,1 -random_initial_guess -random_real -sub_0_pc_bddc_switch_static -use_composite_pc -ksp_monitor -ksp_converged_reason -ksp_type gmres -ksp_view_singularvalues -ksp_view_eigenvalues -sub_0_pc_bddc_use_edges 0 -sub_0_pc_bddc_coarse_pc_type svd -sub_1_ksp_ksp_max_it 1 -sub_1_ksp_ksp_richardson_scale 2.3
781: test:
782: args: -sub_0_pc_bddc_interface_ext_type lump
783: suffix: composite_bddc_lumped
784: test:
785: requires: !single
786: args: -sub_0_pc_bddc_interface_ext_type dirichlet
787: suffix: composite_bddc_dirichlet
789: # GDSW tests
790: testset:
791: nsize: 8
792: filter: grep -v "variant HERMITIAN"
793: args: -cells 7,9,8 -dim 3 -ksp_view -ksp_error_if_not_converged -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space gdsw -mg_levels_pc_type asm -mg_levels_sub_pc_type icc -mg_coarse_redundant_pc_type cholesky
794: test:
795: suffix: gdsw_poisson
796: args: -pde_type Poisson
797: test:
798: requires: mumps !complex
799: suffix: gdsw_poisson_adaptive
800: args: -pde_type Poisson -mg_levels_gdsw_tolerance 0.01 -ksp_monitor_singular_value -mg_levels_gdsw_userdefined {{0 1}separate output} -mg_levels_gdsw_pseudo_pc_type qr
801: test:
802: suffix: gdsw_elast
803: args: -pde_type Elasticity
804: test:
805: requires: hpddm
806: suffix: gdsw_elast_hpddm
807: args: -pde_type Elasticity -mg_levels_gdsw_ksp_type hpddm -mg_levels_gdsw_ksp_hpddm_type cg
808: test:
809: requires: mumps !complex
810: suffix: gdsw_elast_adaptive
811: args: -pde_type Elasticity -mg_levels_gdsw_tolerance 0.01 -ksp_monitor_singular_value -mg_levels_gdsw_userdefined {{0 1}separate output}
813: # Multi-Element tests
814: test:
815: suffix: bddc_multielement_empty
816: output_file: output/ex71_bddc_multi_element.out
817: nsize: 6
818: args: -cells 2,2 -dim 2 -ksp_error_if_not_converged -multi_element -pde_type Poisson -pc_bddc_levels 1 -pc_bddc_coarsening_ratio 2 -pc_bddc_coarse_eqs_limit 0 -ksp_converged_reason
820: test:
821: requires: mumps
822: suffix: bddc_multielement_adaptive_reuse
823: filter: sed -e "s/CONVERGED_RTOL iterations 3/CONVERGED_RTOL iterations 2/g"
824: nsize: 2
825: args: -test_reuse -cells 4,4 -dim 2 -ksp_error_if_not_converged -multi_element -pde_type Poisson -pc_bddc_levels 1 -pc_bddc_coarsening_ratio 2 -pc_bddc_aggregator_petscpartitioner_type simple -pc_bddc_use_deluxe_scaling -pc_bddc_adaptive_threshold 2 -sub_schurs_mat_solver_type mumps -ksp_converged_reason
827: test:
828: nsize: {{1 2 3}}
829: suffix: bddc_multi_element
830: args: -cells 3,3,3 -dim 3 -ksp_error_if_not_converged -multi_element -pde_type {{Poisson Elasticity}} -ksp_converged_reason
832: test:
833: suffix: bddc_multi_square
834: output_file: output/ex71_bddc_multi_element.out
835: args: -cells 2,2 -dim 2 -ksp_error_if_not_converged -multi_element -pc_bddc_local_mat_graph_square 4 -ksp_converged_reason
837: test:
838: suffix: bddc_multielement_5lev_2d
839: filter: sed -e "s/CONVERGED_RTOL iterations 4/CONVERGED_RTOL iterations 3/g"
840: args: -test_reuse -cells 4,4 -dim 2 -ksp_error_if_not_converged -multi_element -pde_type Poisson -pc_bddc_levels 3 -pc_bddc_coarsening_ratio 2 -pc_bddc_aggregator_petscpartitioner_type simple -ksp_converged_reason -pc_bddc_aggregator_petscpartitioner_view -pc_bddc_aggregator_petscpartitioner_view_graph
842: test:
843: suffix: bddc_multielement_5lev_3d
844: args: -test_reuse -cells 3,3,3 -dim 3 -ksp_error_if_not_converged -multi_element -pde_type Poisson -pc_bddc_levels 3 -pc_bddc_coarsening_ratio 3 -pc_bddc_aggregator_petscpartitioner_type simple -ksp_converged_reason -pc_bddc_aggregator_petscpartitioner_view -pc_bddc_aggregator_petscpartitioner_view_graph
846: TEST*/