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: {
     16:   PetscInt       i,m;
     17:   PetscScalar    *xx;
     18:   PetscReal      xd;

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

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

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

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

     71:   PetscInitialize(&argc,&args,NULL,NULL);if (ierr) return ierr;
     72:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);

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

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

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

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

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

    136: /*TEST

    138:   build:
    139:       requires: !complex

    141:    test:

    143: TEST*/