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