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