Actual source code: ex69.c
petsc-3.9.4 2018-09-11
2: #include <petscgll.h>
3: #include <petscdraw.h>
4: #include <petscviewer.h>
5: #include <petscksp.h>
6: #include <petscdmda.h>
8: /*
9: Solves -Laplacian u = f, u(-1) = u(1) = 0 with multiple spectral elements
11: Uses DMDA to manage the parallelization of the elements
13: This is not intended to be highly optimized in either memory usage or time, but strifes for simplicity.
15: */
17: PetscErrorCode ComputeSolution(DM da,PetscGLL *gll,Vec u)
18: {
20: PetscInt j,xs,xn;
21: PetscScalar *uu,*xx;
22: PetscReal xd;
23: Vec x;
26: DMDAGetCorners(da,&xs,NULL,NULL,&xn,NULL,NULL);
27: DMGetCoordinates(da,&x);
28: DMDAVecGetArray(da,x,&xx);
29: DMDAVecGetArray(da,u,&uu);
30: /* loop over local nodes */
31: for (j=xs; j<xs+xn; j++) {
32: xd = xx[j];
33: uu[j] = (xd*xd - 1.0)*PetscCosReal(5.*PETSC_PI*xd);
34: }
35: DMDAVecRestoreArray(da,x,&xx);
36: DMDAVecRestoreArray(da,u,&uu);
37: return(0);
38: }
40: /*
41: 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
42: basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
43: */
44: PetscErrorCode ComputeRhs(DM da,PetscGLL *gll,Vec b)
45: {
47: PetscInt i,j,xs,xn,n = gll->n;
48: PetscScalar *bb,*xx;
49: PetscReal xd;
50: Vec blocal,xlocal;
53: DMDAGetCorners(da,&xs,NULL,NULL,&xn,NULL,NULL);
54: xs = xs/(n-1);
55: xn = xn/(n-1);
56: DMGetLocalVector(da,&blocal);
57: VecZeroEntries(blocal);
58: DMDAVecGetArray(da,blocal,&bb);
59: DMGetCoordinatesLocal(da,&xlocal);
60: DMDAVecGetArray(da,xlocal,&xx);
61: /* loop over local spectral elements */
62: for (j=xs; j<xs+xn; j++) {
63: /* loop over GLL points in each element */
64: for (i=0; i<n; i++) {
65: xd = xx[j*(n-1) + i];
66: bb[j*(n-1) + 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));
67: }
68: }
69: DMDAVecRestoreArray(da,xlocal,&xx);
70: DMDAVecRestoreArray(da,blocal,&bb);
71: VecZeroEntries(b);
72: DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
73: DMLocalToGlobalEnd(da,blocal,ADD_VALUES,b);
74: DMRestoreLocalVector(da,&blocal);
75: return(0);
76: }
78: /*
79: Run with -build_twosided allreduce -pc_type bjacobi -sub_pc_type lu -q 16 -ksp_rtol 1.e-34 (or 1.e-14 for double precision)
81: -q <q> number of spectral elements to use
82: -N <N> maximum number of GLL points per element
84: */
85: int main(int argc,char **args)
86: {
88: PetscGLL gll;
89: PetscInt N = 80,n,q = 8,xs,xn,j,l;
90: PetscReal **A;
91: Mat K;
92: KSP ksp;
93: PC pc;
94: Vec x,b;
95: PetscInt *rows;
96: PetscReal norm,xc,yc,h;
97: PetscScalar *f;
98: PetscDraw draw;
99: PetscDrawLG lg;
100: PetscDrawAxis axis;
101: DM da;
102: PetscMPIInt rank,size;
104: PetscInitialize(&argc,&args,NULL,NULL);if (ierr) return ierr;
105: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
106: MPI_Comm_size(PETSC_COMM_WORLD,&size);
107: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
108: PetscOptionsGetInt(NULL,NULL,"-q",&q,NULL);
110: PetscDrawCreate(PETSC_COMM_WORLD,NULL,"Log(Error norm) vs Number of GLL points",0,0,500,500,&draw);
111: PetscDrawSetFromOptions(draw);
112: PetscDrawLGCreate(draw,1,&lg);
113: PetscDrawLGSetUseMarkers(lg,PETSC_TRUE);
114: PetscDrawLGGetAxis(lg,&axis);
115: PetscDrawAxisSetLabels(axis,NULL,"Number of GLL points","Log(Error Norm)");
117: for (n=4; n<N; n+=2) {
119: /*
120: da contains the information about the parallel layout of the elements
121: */
122: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,q*(n-1)+1,1,1,NULL,&da);
123: DMSetFromOptions(da);
124: DMSetUp(da);
125: DMDAGetInfo(da,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
126: q = (q-1)/(n-1); /* number of spectral elements */
128: /*
129: gll simply contains the GLL node and weight values
130: */
131: PetscGLLCreate(n,PETSCGLL_VIA_LINEARALGEBRA,&gll);
132: DMDASetGLLCoordinates(da,&gll);
134: /*
135: Creates the element stiffness matrix for the given gll
136: */
137: PetscGLLElementLaplacianCreate(&gll,&A);
139: /*
140: Scale the element stiffness and weights by the size of the element
141: */
142: h = 2.0/q;
143: for (j=0; j<n; j++) {
144: gll.weights[j] *= .5*h;
145: for (l=0; l<n; l++) {
146: A[j][l] = 2.*A[j][l]/h;
147: }
148: }
150: /*
151: Create the global stiffness matrix and add the element stiffness for each local element
152: */
153: DMCreateMatrix(da,&K);
154: MatSetOption(K,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
155: DMDAGetCorners(da,&xs,NULL,NULL,&xn,NULL,NULL);
156: xs = xs/(n-1);
157: xn = xn/(n-1);
158: PetscMalloc1(n,&rows);
159: /*
160: loop over local elements
161: */
162: for (j=xs; j<xs+xn; j++) {
163: for (l=0; l<n; l++) rows[l] = j*(n-1)+l;
164: MatSetValues(K,n,rows,n,rows,&A[0][0],ADD_VALUES);
165: }
166: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
167: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
169: MatCreateVecs(K,&x,&b);
170: ComputeRhs(da,&gll,b);
172: /*
173: Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
174: */
175: rows[0] = 0;
176: rows[1] = q*(n-1);
177: MatZeroRowsColumns(K,2,rows,1.0,x,b);
178: PetscFree(rows);
180: KSPCreate(PETSC_COMM_WORLD,&ksp);
181: KSPSetOperators(ksp,K,K);
182: KSPGetPC(ksp,&pc);
183: PCSetType(pc,PCLU);
184: KSPSetFromOptions(ksp);
185: KSPSolve(ksp,b,x);
187: /* compute the error to the continium problem */
188: ComputeSolution(da,&gll,b);
189: VecAXPY(x,-1.0,b);
191: /* compute the L^2 norm of the error */
192: VecGetArray(x,&f);
193: PetscGLLIntegrate(&gll,f,&norm);
194: VecRestoreArray(x,&f);
195: norm = PetscSqrtReal(norm);
196: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"L^2 norm of the error %D %g\n",n,(double)norm);
197: if (n > 10 && norm > 1.e-8) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Slower convergence than expected");
198: xc = (PetscReal)n;
199: yc = PetscLog10Real(norm);
200: PetscDrawLGAddPoint(lg,&xc,&yc);
201: PetscDrawLGDraw(lg);
203: VecDestroy(&b);
204: VecDestroy(&x);
205: KSPDestroy(&ksp);
206: MatDestroy(&K);
207: PetscGLLElementLaplacianDestroy(&gll,&A);
208: PetscGLLDestroy(&gll);
209: DMDestroy(&da);
210: }
211: PetscDrawLGDestroy(&lg);
212: PetscDrawDestroy(&draw);
213: PetscFinalize();
214: return ierr;
215: }
217: /*TEST
219: build:
220: requires: !complex
222: test:
223: requires: !single
225: test:
226: suffix: 2
227: nsize: 2
228: requires: !single superlu_dist
230: TEST*/