Actual source code: ex69.c

petsc-master 2018-06-24
Report Typos and Errors

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