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