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