Actual source code: ex68.c

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