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

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

 66:   PetscInitialize(&argc,&args,NULL,NULL);
 67:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);

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

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

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

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

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

131: /*TEST

133:   build:
134:       requires: !complex

136:    test:

138: TEST*/