Actual source code: ex30.c
1: static char help[] = "Tests DMSLICED operations\n\n";
3: #include <petscdmsliced.h>
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;
19: PetscInitialize(&argc, &argv, 0, help);
20: comm = PETSC_COMM_WORLD;
21: MPI_Comm_size(comm, &size);
22: MPI_Comm_rank(comm, &rank);
24: PetscOptionsBegin(comm, 0, "Options for DMSliced test", 0);
25: {
26: PetscOptionsInt("-n", "Global number of nodes", "", N, &N, NULL);
27: PetscOptionsInt("-bs", "Block size (1 or 2)", "", bs, &bs, NULL);
28: if (bs != 1) {
30: PetscOptionsReal("-alpha", "Inverse time step for wave operator", "", alpha, &alpha, NULL);
31: PetscOptionsReal("-K", "Bulk modulus of compressibility", "", K, &K, NULL);
32: PetscOptionsReal("-rho0", "Reference density", "", rho0, &rho0, NULL);
33: PetscOptionsReal("-u0", "Reference velocity", "", u0, &u0, NULL);
34: PetscOptionsReal("-sigma", "Width of Gaussian density perturbation", "", sigma, &sigma, NULL);
35: PetscOptionsBool("-block", "Use block matrix assembly", "", useblock, &useblock, NULL);
36: }
37: PetscOptionsString("-sliced_mat_type", "Matrix type to use (aij or baij)", "", mat_type, mat_type, sizeof(mat_type), NULL);
38: }
39: PetscOptionsEnd();
41: /* Split ownership, set up periodic grid in 1D */
42: n = PETSC_DECIDE;
43: PetscSplitOwnership(comm, &n, &N);
44: rstart = 0;
45: 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;
50: 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: 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 */
62: if (!useblock) DMSlicedSetBlockFills(slice, dfill, ofill); /* Irrelevant for baij formats */
63: DMSetMatType(slice, mat_type);
64: DMCreateMatrix(slice, &A);
65: PetscFree2(d_nnz, o_nnz);
66: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
68: DMCreateGlobalVector(slice, &x);
69: VecDuplicate(x, &b);
71: VecGhostGetLocalForm(x, &lf);
72: VecGetSize(lf, &m);
74: 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: 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: 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: 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: 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: VecRestoreArray(lf, &xx);
138: VecGhostRestoreLocalForm(x, &lf);
139: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
140: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
142: MatMult(A, x, b);
143: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
144: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
145: VecView(b, PETSC_VIEWER_STDOUT_WORLD);
147: /* Update the ghosted values, view the result on rank 0. */
148: VecGhostUpdateBegin(b, INSERT_VALUES, SCATTER_FORWARD);
149: VecGhostUpdateEnd(b, INSERT_VALUES, SCATTER_FORWARD);
150: if (rank == 0) {
151: VecGhostGetLocalForm(b, &lf);
152: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Local form of b on rank 0, last two nodes are ghost nodes\n");
153: VecView(lf, PETSC_VIEWER_STDOUT_SELF);
154: VecGhostRestoreLocalForm(b, &lf);
155: }
157: DMDestroy(&slice);
158: VecDestroy(&x);
159: VecDestroy(&b);
160: MatDestroy(&A);
161: PetscFinalize();
162: return 0;
163: }
165: /*TEST
167: test:
168: nsize: 2
169: args: -bs 2 -block 0 -sliced_mat_type baij -alpha 10 -u0 0.1
171: test:
172: suffix: 2
173: nsize: 2
174: args: -bs 2 -block 1 -sliced_mat_type aij -alpha 10 -u0 0.1
176: test:
177: suffix: 3
178: nsize: 2
179: args: -bs 2 -block 0 -sliced_mat_type aij -alpha 10 -u0 0.1
181: TEST*/