Actual source code: ex69.c
2: #include <petscdt.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: typedef struct {
18: PetscInt n; /* number of nodes */
19: PetscReal *nodes; /* GLL nodes */
20: PetscReal *weights; /* GLL weights */
21: } PetscGLL;
23: PetscErrorCode ComputeSolution(DM da,PetscGLL *gll,Vec u)
24: {
25: PetscInt j,xs,xn;
26: PetscScalar *uu,*xx;
27: PetscReal xd;
28: Vec x;
30: DMDAGetCorners(da,&xs,NULL,NULL,&xn,NULL,NULL);
31: DMGetCoordinates(da,&x);
32: DMDAVecGetArray(da,x,&xx);
33: DMDAVecGetArray(da,u,&uu);
34: /* loop over local nodes */
35: for (j=xs; j<xs+xn; j++) {
36: xd = xx[j];
37: uu[j] = (xd*xd - 1.0)*PetscCosReal(5.*PETSC_PI*xd);
38: }
39: DMDAVecRestoreArray(da,x,&xx);
40: DMDAVecRestoreArray(da,u,&uu);
41: return 0;
42: }
44: /*
45: 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
46: basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
47: */
48: PetscErrorCode ComputeRhs(DM da,PetscGLL *gll,Vec b)
49: {
50: PetscInt i,j,xs,xn,n = gll->n;
51: PetscScalar *bb,*xx;
52: PetscReal xd;
53: Vec blocal,xlocal;
55: DMDAGetCorners(da,&xs,NULL,NULL,&xn,NULL,NULL);
56: xs = xs/(n-1);
57: xn = xn/(n-1);
58: DMGetLocalVector(da,&blocal);
59: VecZeroEntries(blocal);
60: DMDAVecGetArray(da,blocal,&bb);
61: DMGetCoordinatesLocal(da,&xlocal);
62: DMDAVecGetArray(da,xlocal,&xx);
63: /* loop over local spectral elements */
64: for (j=xs; j<xs+xn; j++) {
65: /* loop over GLL points in each element */
66: for (i=0; i<n; i++) {
67: xd = xx[j*(n-1) + i];
68: 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));
69: }
70: }
71: DMDAVecRestoreArray(da,xlocal,&xx);
72: DMDAVecRestoreArray(da,blocal,&bb);
73: VecZeroEntries(b);
74: DMLocalToGlobalBegin(da,blocal,ADD_VALUES,b);
75: DMLocalToGlobalEnd(da,blocal,ADD_VALUES,b);
76: DMRestoreLocalVector(da,&blocal);
77: return 0;
78: }
80: /*
81: 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)
83: -q <q> number of spectral elements to use
84: -N <N> maximum number of GLL points per element
86: */
87: int main(int argc,char **args)
88: {
89: PetscGLL gll;
90: PetscInt N = 80,n,q = 8,xs,xn,j,l;
91: PetscReal **A;
92: Mat K;
93: KSP ksp;
94: PC pc;
95: Vec x,b;
96: PetscInt *rows;
97: PetscReal norm,xc,yc,h;
98: PetscScalar *f;
99: PetscDraw draw;
100: PetscDrawLG lg;
101: PetscDrawAxis axis;
102: DM da;
103: PetscMPIInt rank,size;
105: PetscInitialize(&argc,&args,NULL,NULL);
106: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
107: MPI_Comm_size(PETSC_COMM_WORLD,&size);
108: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
109: PetscOptionsGetInt(NULL,NULL,"-q",&q,NULL);
111: PetscDrawCreate(PETSC_COMM_WORLD,NULL,"Log(Error norm) vs Number of GLL points",0,0,500,500,&draw);
112: PetscDrawSetFromOptions(draw);
113: PetscDrawLGCreate(draw,1,&lg);
114: PetscDrawLGSetUseMarkers(lg,PETSC_TRUE);
115: PetscDrawLGGetAxis(lg,&axis);
116: PetscDrawAxisSetLabels(axis,NULL,"Number of GLL points","Log(Error Norm)");
118: for (n=4; n<N; n+=2) {
120: /*
121: da contains the information about the parallel layout of the elements
122: */
123: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,q*(n-1)+1,1,1,NULL,&da);
124: DMSetFromOptions(da);
125: DMSetUp(da);
126: DMDAGetInfo(da,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
127: q = (q-1)/(n-1); /* number of spectral elements */
129: /*
130: gll simply contains the GLL node and weight values
131: */
132: PetscMalloc2(n,&gll.nodes,n,&gll.weights);
133: PetscDTGaussLobattoLegendreQuadrature(n,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,gll.nodes,gll.weights);
134: gll.n = n;
135: DMDASetGLLCoordinates(da,gll.n,gll.nodes);
137: /*
138: Creates the element stiffness matrix for the given gll
139: */
140: PetscGaussLobattoLegendreElementLaplacianCreate(gll.n,gll.nodes,gll.weights,&A);
142: /*
143: Scale the element stiffness and weights by the size of the element
144: */
145: h = 2.0/q;
146: for (j=0; j<n; j++) {
147: gll.weights[j] *= .5*h;
148: for (l=0; l<n; l++) {
149: A[j][l] = 2.*A[j][l]/h;
150: }
151: }
153: /*
154: Create the global stiffness matrix and add the element stiffness for each local element
155: */
156: DMCreateMatrix(da,&K);
157: MatSetOption(K,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
158: DMDAGetCorners(da,&xs,NULL,NULL,&xn,NULL,NULL);
159: xs = xs/(n-1);
160: xn = xn/(n-1);
161: PetscMalloc1(n,&rows);
162: /*
163: loop over local elements
164: */
165: for (j=xs; j<xs+xn; j++) {
166: for (l=0; l<n; l++) rows[l] = j*(n-1)+l;
167: MatSetValues(K,n,rows,n,rows,&A[0][0],ADD_VALUES);
168: }
169: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
172: MatCreateVecs(K,&x,&b);
173: ComputeRhs(da,&gll,b);
175: /*
176: Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
177: */
178: rows[0] = 0;
179: rows[1] = q*(n-1);
180: MatZeroRowsColumns(K,2,rows,1.0,x,b);
181: PetscFree(rows);
183: KSPCreate(PETSC_COMM_WORLD,&ksp);
184: KSPSetOperators(ksp,K,K);
185: KSPGetPC(ksp,&pc);
186: PCSetType(pc,PCLU);
187: KSPSetFromOptions(ksp);
188: KSPSolve(ksp,b,x);
190: /* compute the error to the continium problem */
191: ComputeSolution(da,&gll,b);
192: VecAXPY(x,-1.0,b);
194: /* compute the L^2 norm of the error */
195: VecGetArray(x,&f);
196: PetscGaussLobattoLegendreIntegrate(gll.n,gll.nodes,gll.weights,f,&norm);
197: VecRestoreArray(x,&f);
198: norm = PetscSqrtReal(norm);
199: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"L^2 norm of the error %D %g\n",n,(double)norm);
201: xc = (PetscReal)n;
202: yc = PetscLog10Real(norm);
203: PetscDrawLGAddPoint(lg,&xc,&yc);
204: PetscDrawLGDraw(lg);
206: VecDestroy(&b);
207: VecDestroy(&x);
208: KSPDestroy(&ksp);
209: MatDestroy(&K);
210: PetscGaussLobattoLegendreElementLaplacianDestroy(gll.n,gll.nodes,gll.weights,&A);
211: PetscFree2(gll.nodes,gll.weights);
212: DMDestroy(&da);
213: }
214: PetscDrawLGDestroy(&lg);
215: PetscDrawDestroy(&draw);
216: PetscFinalize();
217: return 0;
218: }
220: /*TEST
222: build:
223: requires: !complex
225: test:
226: requires: !single
228: test:
229: suffix: 2
230: nsize: 2
231: requires: !single superlu_dist
233: TEST*/