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:   PetscFunctionBegin;
 20:   PetscCall(VecGetSize(x, &m));
 21:   PetscCall(VecGetArray(x, &xx));
 22:   for (i = 0; i < m; i++) {
 23:     xd    = nodes[i];
 24:     xx[i] = (xd * xd - 1.0) * PetscCosReal(5. * PETSC_PI * xd);
 25:   }
 26:   PetscCall(VecRestoreArray(x, &xx));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

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

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

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

 68:   PetscFunctionBeginUser;
 69:   PetscCall(PetscInitialize(&argc, &args, NULL, NULL));
 70:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));

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

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

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

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

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

134: /*TEST

136:   build:
137:       requires: !complex

139:    test:

141: TEST*/