Actual source code: ex69.c


  2: #include <petscdt.h>
  3: #include <petscdraw.h>
  4: #include <petscviewer.h>
  5: #include <petscksp.h>
  6: #include <petscdmda.h>

  8: /*
  9:       Solves -Laplacian u = f,  u(-1) = u(1) = 0 with multiple spectral elements

 11:       Uses DMDA to manage the parallelization of the elements

 13:       This is not intended to be highly optimized in either memory usage or time, but strifes for simplicity.

 15: */

 17: typedef struct {
 18:   PetscInt   n;       /* number of nodes */
 19:   PetscReal *nodes;   /* GLL nodes */
 20:   PetscReal *weights; /* GLL weights */
 21: } PetscGLL;

 23: PetscErrorCode ComputeSolution(DM da, PetscGLL *gll, Vec u)
 24: {
 25:   PetscInt     j, xs, xn;
 26:   PetscScalar *uu, *xx;
 27:   PetscReal    xd;
 28:   Vec          x;

 30:   DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL);
 31:   DMGetCoordinates(da, &x);
 32:   DMDAVecGetArray(da, x, &xx);
 33:   DMDAVecGetArray(da, u, &uu);
 34:   /* loop over local nodes */
 35:   for (j = xs; j < xs + xn; j++) {
 36:     xd    = xx[j];
 37:     uu[j] = (xd * xd - 1.0) * PetscCosReal(5. * PETSC_PI * xd);
 38:   }
 39:   DMDAVecRestoreArray(da, x, &xx);
 40:   DMDAVecRestoreArray(da, u, &uu);
 41:   return 0;
 42: }

 44: /*
 45:       Evaluates \integral_{-1}^{1} f*v_i  where v_i is the ith basis polynomial via the GLL nodes and weights, since the v_i
 46:       basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
 47: */
 48: PetscErrorCode ComputeRhs(DM da, PetscGLL *gll, Vec b)
 49: {
 50:   PetscInt     i, j, xs, xn, n = gll->n;
 51:   PetscScalar *bb, *xx;
 52:   PetscReal    xd;
 53:   Vec          blocal, xlocal;

 55:   DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL);
 56:   xs = xs / (n - 1);
 57:   xn = xn / (n - 1);
 58:   DMGetLocalVector(da, &blocal);
 59:   VecZeroEntries(blocal);
 60:   DMDAVecGetArray(da, blocal, &bb);
 61:   DMGetCoordinatesLocal(da, &xlocal);
 62:   DMDAVecGetArray(da, xlocal, &xx);
 63:   /* loop over local spectral elements */
 64:   for (j = xs; j < xs + xn; j++) {
 65:     /* loop over GLL points in each element */
 66:     for (i = 0; i < n; i++) {
 67:       xd = xx[j * (n - 1) + i];
 68:       bb[j * (n - 1) + i] += -gll->weights[i] * (-20. * PETSC_PI * xd * PetscSinReal(5. * PETSC_PI * xd) + (2. - (5. * PETSC_PI) * (5. * PETSC_PI) * (xd * xd - 1.)) * PetscCosReal(5. * PETSC_PI * xd));
 69:     }
 70:   }
 71:   DMDAVecRestoreArray(da, xlocal, &xx);
 72:   DMDAVecRestoreArray(da, blocal, &bb);
 73:   VecZeroEntries(b);
 74:   DMLocalToGlobalBegin(da, blocal, ADD_VALUES, b);
 75:   DMLocalToGlobalEnd(da, blocal, ADD_VALUES, b);
 76:   DMRestoreLocalVector(da, &blocal);
 77:   return 0;
 78: }

 80: /*
 81:      Run with -build_twosided allreduce -pc_type bjacobi -sub_pc_type lu -q 16 -ksp_rtol 1.e-34 (or 1.e-14 for double precision)

 83:      -q <q> number of spectral elements to use
 84:      -N <N> maximum number of GLL points per element

 86: */
 87: int main(int argc, char **args)
 88: {
 89:   PetscGLL      gll;
 90:   PetscInt      N = 80, n, q = 8, xs, xn, j, l;
 91:   PetscReal   **A;
 92:   Mat           K;
 93:   KSP           ksp;
 94:   PC            pc;
 95:   Vec           x, b;
 96:   PetscInt     *rows;
 97:   PetscReal     norm, xc, yc, h;
 98:   PetscScalar  *f;
 99:   PetscDraw     draw;
100:   PetscDrawLG   lg;
101:   PetscDrawAxis axis;
102:   DM            da;
103:   PetscMPIInt   rank, size;

106:   PetscInitialize(&argc, &args, NULL, NULL);
107:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
108:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
109:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
110:   PetscOptionsGetInt(NULL, NULL, "-q", &q, NULL);

112:   PetscDrawCreate(PETSC_COMM_WORLD, NULL, "Log(Error norm) vs Number of GLL points", 0, 0, 500, 500, &draw);
113:   PetscDrawSetFromOptions(draw);
114:   PetscDrawLGCreate(draw, 1, &lg);
115:   PetscDrawLGSetUseMarkers(lg, PETSC_TRUE);
116:   PetscDrawLGGetAxis(lg, &axis);
117:   PetscDrawAxisSetLabels(axis, NULL, "Number of GLL points", "Log(Error Norm)");

119:   for (n = 4; n < N; n += 2) {
120:     /*
121:        da contains the information about the parallel layout of the elements
122:     */
123:     DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, q * (n - 1) + 1, 1, 1, NULL, &da);
124:     DMSetFromOptions(da);
125:     DMSetUp(da);
126:     DMDAGetInfo(da, NULL, &q, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
127:     q = (q - 1) / (n - 1); /* number of spectral elements */

129:     /*
130:        gll simply contains the GLL node and weight values
131:     */
132:     PetscMalloc2(n, &gll.nodes, n, &gll.weights);
133:     PetscDTGaussLobattoLegendreQuadrature(n, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, gll.nodes, gll.weights);
134:     gll.n = n;
135:     DMDASetGLLCoordinates(da, gll.n, gll.nodes);

137:     /*
138:        Creates the element stiffness matrix for the given gll
139:     */
140:     PetscGaussLobattoLegendreElementLaplacianCreate(gll.n, gll.nodes, gll.weights, &A);

142:     /*
143:       Scale the element stiffness and weights by the size of the element
144:     */
145:     h = 2.0 / q;
146:     for (j = 0; j < n; j++) {
147:       gll.weights[j] *= .5 * h;
148:       for (l = 0; l < n; l++) A[j][l] = 2. * A[j][l] / h;
149:     }

151:     /*
152:         Create the global stiffness matrix and add the element stiffness for each local element
153:     */
154:     DMCreateMatrix(da, &K);
155:     MatSetOption(K, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
156:     DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL);
157:     xs = xs / (n - 1);
158:     xn = xn / (n - 1);
159:     PetscMalloc1(n, &rows);
160:     /*
161:         loop over local elements
162:     */
163:     for (j = xs; j < xs + xn; j++) {
164:       for (l = 0; l < n; l++) rows[l] = j * (n - 1) + l;
165:       MatSetValues(K, n, rows, n, rows, &A[0][0], ADD_VALUES);
166:     }
167:     MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY);
168:     MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY);

170:     MatCreateVecs(K, &x, &b);
171:     ComputeRhs(da, &gll, b);

173:     /*
174:         Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
175:     */
176:     rows[0] = 0;
177:     rows[1] = q * (n - 1);
178:     MatZeroRowsColumns(K, 2, rows, 1.0, x, b);
179:     PetscFree(rows);

181:     KSPCreate(PETSC_COMM_WORLD, &ksp);
182:     KSPSetOperators(ksp, K, K);
183:     KSPGetPC(ksp, &pc);
184:     PCSetType(pc, PCLU);
185:     KSPSetFromOptions(ksp);
186:     KSPSolve(ksp, b, x);

188:     /* compute the error to the continium problem */
189:     ComputeSolution(da, &gll, b);
190:     VecAXPY(x, -1.0, b);

192:     /* compute the L^2 norm of the error */
193:     VecGetArray(x, &f);
194:     PetscGaussLobattoLegendreIntegrate(gll.n, gll.nodes, gll.weights, f, &norm);
195:     VecRestoreArray(x, &f);
196:     norm = PetscSqrtReal(norm);
197:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "L^2 norm of the error %" PetscInt_FMT " %g\n", n, (double)norm);
199:     xc = (PetscReal)n;
200:     yc = PetscLog10Real(norm);
201:     PetscDrawLGAddPoint(lg, &xc, &yc);
202:     PetscDrawLGDraw(lg);

204:     VecDestroy(&b);
205:     VecDestroy(&x);
206:     KSPDestroy(&ksp);
207:     MatDestroy(&K);
208:     PetscGaussLobattoLegendreElementLaplacianDestroy(gll.n, gll.nodes, gll.weights, &A);
209:     PetscFree2(gll.nodes, gll.weights);
210:     DMDestroy(&da);
211:   }
212:   PetscDrawLGDestroy(&lg);
213:   PetscDrawDestroy(&draw);
214:   PetscFinalize();
215:   return 0;
216: }

218: /*TEST

220:   build:
221:       requires: !complex

223:    test:
224:      requires: !single

226:    test:
227:      suffix: 2
228:      nsize: 2
229:      requires: superlu_dist

231: TEST*/