Actual source code: ex68.c
petsc-3.13.6 2020-09-29
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*/