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*/