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*/