Actual source code: ex68.c


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

  7: /*
  8:       Solves -Laplacian u = f,  u(-1) = u(1) = 0 with a single spectral element for n = 4 to N GLL points

 10:       Plots the L_2 norm of the error (evaluated via the GLL nodes and weights) as a function of n.

 12: */
 13: PetscErrorCode ComputeSolution(PetscInt n, PetscReal *nodes, PetscReal *weights, Vec x)
 14: {
 15:   PetscInt     i, m;
 16:   PetscScalar *xx;
 17:   PetscReal    xd;

 19:   VecGetSize(x, &m);
 20:   VecGetArray(x, &xx);
 21:   for (i = 0; i < m; i++) {
 22:     xd    = nodes[i];
 23:     xx[i] = (xd * xd - 1.0) * PetscCosReal(5. * PETSC_PI * xd);
 24:   }
 25:   VecRestoreArray(x, &xx);
 26:   return 0;
 27: }

 29: /*
 30:       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
 31:       basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
 32: */
 33: PetscErrorCode ComputeRhs(PetscInt n, PetscReal *nodes, PetscReal *weights, Vec b)
 34: {
 35:   PetscInt     i, m;
 36:   PetscScalar *bb;
 37:   PetscReal    xd;

 39:   VecGetSize(b, &m);
 40:   VecGetArray(b, &bb);
 41:   for (i = 0; i < m; i++) {
 42:     xd    = nodes[i];
 43:     bb[i] = -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));
 44:   }
 45:   VecRestoreArray(b, &bb);
 46:   return 0;
 47: }

 49: int main(int argc, char **args)
 50: {
 51:   PetscReal    *nodes;
 52:   PetscReal    *weights;
 53:   PetscInt      N = 80, n;
 54:   PetscReal   **A;
 55:   Mat           K;
 56:   KSP           ksp;
 57:   PC            pc;
 58:   Vec           x, b;
 59:   PetscInt      rows[2];
 60:   PetscReal     norm, xc, yc;
 61:   PetscScalar  *f;
 62:   PetscDraw     draw;
 63:   PetscDrawLG   lg;
 64:   PetscDrawAxis axis;

 67:   PetscInitialize(&argc, &args, NULL, NULL);
 68:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);

 70:   PetscDrawCreate(PETSC_COMM_SELF, NULL, "Log(Error norm) vs Number of GLL points", 0, 0, 500, 500, &draw);
 71:   PetscDrawSetFromOptions(draw);
 72:   PetscDrawLGCreate(draw, 1, &lg);
 73:   PetscDrawLGSetUseMarkers(lg, PETSC_TRUE);
 74:   PetscDrawLGGetAxis(lg, &axis);
 75:   PetscDrawAxisSetLabels(axis, NULL, "Number of GLL points", "Log(Error Norm)");

 77:   for (n = 4; n < N; n += 2) {
 78:     /*
 79:        compute GLL node and weight values
 80:     */
 81:     PetscMalloc2(n, &nodes, n, &weights);
 82:     PetscDTGaussLobattoLegendreQuadrature(n, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, nodes, weights);
 83:     /*
 84:        Creates the element stiffness matrix for the given gll
 85:     */
 86:     PetscGaussLobattoLegendreElementLaplacianCreate(n, nodes, weights, &A);
 87:     MatCreateSeqDense(PETSC_COMM_SELF, n, n, &A[0][0], &K);
 88:     rows[0] = 0;
 89:     rows[1] = n - 1;
 90:     KSPCreate(PETSC_COMM_SELF, &ksp);
 91:     MatCreateVecs(K, &x, &b);
 92:     ComputeRhs(n, nodes, weights, b);
 93:     /*
 94:         Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
 95:     */
 96:     MatZeroRowsColumns(K, 2, rows, 1.0, x, b);
 97:     KSPSetOperators(ksp, K, K);
 98:     KSPGetPC(ksp, &pc);
 99:     PCSetType(pc, PCLU);
100:     KSPSetFromOptions(ksp);
101:     KSPSolve(ksp, b, x);

103:     /* compute the error to the continium problem */
104:     ComputeSolution(n, nodes, weights, b);
105:     VecAXPY(x, -1.0, b);

107:     /* compute the L^2 norm of the error */
108:     VecGetArray(x, &f);
109:     PetscGaussLobattoLegendreIntegrate(n, nodes, weights, f, &norm);
110:     VecRestoreArray(x, &f);
111:     norm = PetscSqrtReal(norm);
112:     PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "L^2 norm of the error %" PetscInt_FMT " %g\n", n, (double)norm);
113:     xc = (PetscReal)n;
114:     yc = PetscLog10Real(norm);
115:     PetscDrawLGAddPoint(lg, &xc, &yc);
116:     PetscDrawLGDraw(lg);

118:     VecDestroy(&b);
119:     VecDestroy(&x);
120:     KSPDestroy(&ksp);
121:     MatDestroy(&K);
122:     PetscGaussLobattoLegendreElementLaplacianDestroy(n, nodes, weights, &A);
123:     PetscFree2(nodes, weights);
124:   }
125:   PetscDrawSetPause(draw, -2);
126:   PetscDrawLGDestroy(&lg);
127:   PetscDrawDestroy(&draw);
128:   PetscFinalize();
129:   return 0;
130: }

132: /*TEST

134:   build:
135:       requires: !complex

137:    test:

139: TEST*/