Actual source code: ex68.c
petsc-3.10.5 2019-03-28
2: #include <petscgll.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(PetscGLL *gll,Vec x)
14: {
16: PetscInt i,n;
17: PetscScalar *xx;
18: PetscReal xd;
21: VecGetSize(x,&n);
22: VecGetArray(x,&xx);
23: for (i=0; i<n; i++) {
24: xd = gll->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(PetscGLL *gll,Vec b)
36: {
38: PetscInt i,n;
39: PetscScalar *bb;
40: PetscReal xd;
43: VecGetSize(b,&n);
44: VecGetArray(b,&bb);
45: for (i=0; i<n; i++) {
46: xd = gll->nodes[i];
47: bb[i] = -gll->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: PetscGLL gll;
57: PetscInt N = 80,n;
58: PetscReal **A;
59: Mat K;
60: KSP ksp;
61: PC pc;
62: Vec x,b;
63: PetscInt rows[2];
64: PetscReal norm,xc,yc;
65: PetscScalar *f;
66: PetscDraw draw;
67: PetscDrawLG lg;
68: PetscDrawAxis axis;
70: PetscInitialize(&argc,&args,NULL,NULL);if (ierr) return ierr;
71: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
73: PetscDrawCreate(PETSC_COMM_SELF,NULL,"Log(Error norm) vs Number of GLL points",0,0,500,500,&draw);
74: PetscDrawSetFromOptions(draw);
75: PetscDrawLGCreate(draw,1,&lg);
76: PetscDrawLGSetUseMarkers(lg,PETSC_TRUE);
77: PetscDrawLGGetAxis(lg,&axis);
78: PetscDrawAxisSetLabels(axis,NULL,"Number of GLL points","Log(Error Norm)");
80: for (n=4; n<N; n+=2) {
81: /*
82: gll simply contains the GLL node and weight values
83: */
84: PetscGLLCreate(n,PETSCGLL_VIA_LINEARALGEBRA,&gll);
85: /*
86: Creates the element stiffness matrix for the given gll
87: */
88: PetscGLLElementLaplacianCreate(&gll,&A);
89: MatCreateSeqDense(PETSC_COMM_SELF,n,n,&A[0][0],&K);
90: rows[0] = 0;
91: rows[1] = n-1;
92: KSPCreate(PETSC_COMM_SELF,&ksp);
93: MatCreateVecs(K,&x,&b);
94: ComputeRhs(&gll,b);
95: /*
96: Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
97: */
98: MatZeroRowsColumns(K,2,rows,1.0,x,b);
99: KSPSetOperators(ksp,K,K);
100: KSPGetPC(ksp,&pc);
101: PCSetType(pc,PCLU);
102: KSPSetFromOptions(ksp);
103: KSPSolve(ksp,b,x);
105: /* compute the error to the continium problem */
106: ComputeSolution(&gll,b);
107: VecAXPY(x,-1.0,b);
109: /* compute the L^2 norm of the error */
110: VecGetArray(x,&f);
111: PetscGLLIntegrate(&gll,f,&norm);
112: VecRestoreArray(x,&f);
113: norm = PetscSqrtReal(norm);
114: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF,"L^2 norm of the error %D %g\n",n,(double)norm);
115: xc = (PetscReal)n;
116: yc = PetscLog10Real(norm);
117: PetscDrawLGAddPoint(lg,&xc,&yc);
118: PetscDrawLGDraw(lg);
120: VecDestroy(&b);
121: VecDestroy(&x);
122: KSPDestroy(&ksp);
123: MatDestroy(&K);
124: PetscGLLElementLaplacianDestroy(&gll,&A);
125: PetscGLLDestroy(&gll);
126: }
127: PetscDrawSetPause(draw,-2);
128: PetscDrawLGDestroy(&lg);
129: PetscDrawDestroy(&draw);
130: PetscFinalize();
131: return ierr;
132: }
134: /*TEST
136: build:
137: requires: !complex
139: test:
141: TEST*/