1 static char help[] = "Tests DMSLICED operations\n\n";
2
3 #include <petscdmsliced.h>
4
main(int argc,char * argv[])5 int main(int argc, char *argv[])
6 {
7 char mat_type[256] = MATAIJ; /* default matrix type */
8 MPI_Comm comm;
9 PetscMPIInt rank, size;
10 DM slice;
11 PetscInt i, bs = 1, N = 5, n, m, rstart, ghosts[2], *d_nnz, *o_nnz, dfill[4] = {1, 0, 0, 1}, ofill[4] = {1, 1, 1, 1};
12 PetscReal alpha = 1, K = 1, rho0 = 1, u0 = 0, sigma = 0.2;
13 PetscBool useblock = PETSC_TRUE;
14 PetscScalar *xx;
15 Mat A;
16 Vec x, b, lf;
17
18 PetscFunctionBeginUser;
19 PetscCall(PetscInitialize(&argc, &argv, 0, help));
20 comm = PETSC_COMM_WORLD;
21 PetscCallMPI(MPI_Comm_size(comm, &size));
22 PetscCallMPI(MPI_Comm_rank(comm, &rank));
23
24 PetscOptionsBegin(comm, 0, "Options for DMSliced test", 0);
25 {
26 PetscCall(PetscOptionsInt("-n", "Global number of nodes", "", N, &N, NULL));
27 PetscCall(PetscOptionsInt("-bs", "Block size (1 or 2)", "", bs, &bs, NULL));
28 if (bs != 1) {
29 PetscCheck(bs == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Block size must be 1 or 2");
30 PetscCall(PetscOptionsReal("-alpha", "Inverse time step for wave operator", "", alpha, &alpha, NULL));
31 PetscCall(PetscOptionsReal("-K", "Bulk modulus of compressibility", "", K, &K, NULL));
32 PetscCall(PetscOptionsReal("-rho0", "Reference density", "", rho0, &rho0, NULL));
33 PetscCall(PetscOptionsReal("-u0", "Reference velocity", "", u0, &u0, NULL));
34 PetscCall(PetscOptionsReal("-sigma", "Width of Gaussian density perturbation", "", sigma, &sigma, NULL));
35 PetscCall(PetscOptionsBool("-block", "Use block matrix assembly", "", useblock, &useblock, NULL));
36 }
37 PetscCall(PetscOptionsString("-sliced_mat_type", "Matrix type to use (aij or baij)", "", mat_type, mat_type, sizeof(mat_type), NULL));
38 }
39 PetscOptionsEnd();
40
41 /* Split ownership, set up periodic grid in 1D */
42 n = PETSC_DECIDE;
43 PetscCall(PetscSplitOwnership(comm, &n, &N));
44 rstart = 0;
45 PetscCallMPI(MPI_Scan(&n, &rstart, 1, MPIU_INT, MPI_SUM, comm));
46 rstart -= n;
47 ghosts[0] = (N + rstart - 1) % N;
48 ghosts[1] = (rstart + n) % N;
49
50 PetscCall(PetscMalloc2(n, &d_nnz, n, &o_nnz));
51 for (i = 0; i < n; i++) {
52 if (size > 1 && (i == 0 || i == n - 1)) {
53 d_nnz[i] = 2;
54 o_nnz[i] = 1;
55 } else {
56 d_nnz[i] = 3;
57 o_nnz[i] = 0;
58 }
59 }
60 PetscCall(DMSlicedCreate(comm, bs, n, 2, ghosts, d_nnz, o_nnz, &slice)); /* Currently does not copy X_nnz so we can't free them until after DMSlicedGetMatrix */
61
62 if (!useblock) PetscCall(DMSlicedSetBlockFills(slice, dfill, ofill)); /* Irrelevant for baij formats */
63 PetscCall(DMSetMatType(slice, mat_type));
64 PetscCall(DMCreateMatrix(slice, &A));
65 PetscCall(PetscFree2(d_nnz, o_nnz));
66 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
67
68 PetscCall(DMCreateGlobalVector(slice, &x));
69 PetscCall(VecDuplicate(x, &b));
70
71 PetscCall(VecGhostGetLocalForm(x, &lf));
72 PetscCall(VecGetSize(lf, &m));
73 PetscCheck(m == (n + 2) * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size of local form %" PetscInt_FMT ", expected %" PetscInt_FMT, m, (n + 2) * bs);
74 PetscCall(VecGetArray(lf, &xx));
75 for (i = 0; i < n; i++) {
76 PetscInt row[2], col[9], im, ip;
77 PetscScalar v[12];
78 const PetscReal xref = 2.0 * (rstart + i) / N - 1; /* [-1,1] */
79 const PetscReal h = 1.0 / N; /* grid spacing */
80 im = (i == 0) ? n : i - 1;
81 ip = (i == n - 1) ? n + 1 : i + 1;
82 switch (bs) {
83 case 1: /* Laplacian with periodic boundaries */
84 col[0] = im;
85 col[1] = i;
86 col[2] = ip;
87 v[0] = -h;
88 v[1] = 2 * h;
89 v[2] = -h;
90 PetscCall(MatSetValuesLocal(A, 1, &i, 3, col, v, INSERT_VALUES));
91 xx[i] = PetscSinReal(xref * PETSC_PI);
92 break;
93 case 2: /* Linear acoustic wave operator in variables [rho, u], central differences, periodic, timestep 1/alpha */
94 v[0] = -0.5 * u0;
95 v[1] = -0.5 * K;
96 v[2] = alpha;
97 v[3] = 0;
98 v[4] = 0.5 * u0;
99 v[5] = 0.5 * K;
100 v[6] = -0.5 / rho0;
101 v[7] = -0.5 * u0;
102 v[8] = 0;
103 v[9] = alpha;
104 v[10] = 0.5 / rho0;
105 v[11] = 0.5 * u0;
106 if (useblock) {
107 row[0] = i;
108 col[0] = im;
109 col[1] = i;
110 col[2] = ip;
111 PetscCall(MatSetValuesBlockedLocal(A, 1, row, 3, col, v, INSERT_VALUES));
112 } else {
113 row[0] = 2 * i;
114 row[1] = 2 * i + 1;
115 col[0] = 2 * im;
116 col[1] = 2 * im + 1;
117 col[2] = 2 * i;
118 col[3] = 2 * ip;
119 col[4] = 2 * ip + 1;
120 v[3] = v[4];
121 v[4] = v[5]; /* pack values in first row */
122 PetscCall(MatSetValuesLocal(A, 1, row, 5, col, v, INSERT_VALUES));
123 col[2] = 2 * i + 1;
124 v[8] = v[9];
125 v[9] = v[10];
126 v[10] = v[11]; /* pack values in second row */
127 PetscCall(MatSetValuesLocal(A, 1, row + 1, 5, col, v + 6, INSERT_VALUES));
128 }
129 /* Set current state (gaussian density perturbation) */
130 xx[2 * i] = 0.2 * PetscExpReal(-PetscSqr(xref) / (2 * PetscSqr(sigma)));
131 xx[2 * i + 1] = 0;
132 break;
133 default:
134 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "not implemented for block size %" PetscInt_FMT, bs);
135 }
136 }
137 PetscCall(VecRestoreArray(lf, &xx));
138 PetscCall(VecGhostRestoreLocalForm(x, &lf));
139 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
140 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
141
142 PetscCall(MatMult(A, x, b));
143 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
144 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
145 PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));
146
147 /* Update the ghosted values, view the result on rank 0. */
148 PetscCall(VecGhostUpdateBegin(b, INSERT_VALUES, SCATTER_FORWARD));
149 PetscCall(VecGhostUpdateEnd(b, INSERT_VALUES, SCATTER_FORWARD));
150 if (rank == 0) {
151 PetscCall(VecGhostGetLocalForm(b, &lf));
152 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Local form of b on rank 0, last two nodes are ghost nodes\n"));
153 PetscCall(VecView(lf, PETSC_VIEWER_STDOUT_SELF));
154 PetscCall(VecGhostRestoreLocalForm(b, &lf));
155 }
156
157 PetscCall(DMDestroy(&slice));
158 PetscCall(VecDestroy(&x));
159 PetscCall(VecDestroy(&b));
160 PetscCall(MatDestroy(&A));
161 PetscCall(PetscFinalize());
162 return 0;
163 }
164
165 /*TEST
166
167 test:
168 nsize: 2
169 args: -bs 2 -block 0 -sliced_mat_type baij -alpha 10 -u0 0.1
170
171 test:
172 suffix: 2
173 nsize: 2
174 args: -bs 2 -block 1 -sliced_mat_type aij -alpha 10 -u0 0.1
175
176 test:
177 suffix: 3
178 nsize: 2
179 args: -bs 2 -block 0 -sliced_mat_type aij -alpha 10 -u0 0.1
180
181 TEST*/
182