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

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

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

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

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

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

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

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

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

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

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

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

133: /*TEST

135:   build:
136:       requires: !complex

138:    test:

140: TEST*/