Actual source code: ex74.c

petsc-3.12.3 2020-01-03
Report Typos and Errors
  1: static char help[] = "\
2: Solves the constant-coefficient 1D Heat equation with an Implicit   \n\
3: Runge-Kutta method using MatKAIJ.                                   \n\
4:                                                                     \n\
5:     du      d^2 u                                                   \n\
6:     --  = a ----- ; 0 <= x <= 1;                                    \n\
7:     dt      dx^2                                                    \n\
8:                                                                     \n\
9:   with periodic boundary conditions                                 \n\
10:                                                                     \n\
11: 2nd order central discretization in space:                          \n\
12:                                                                     \n\
13:    [ d^2 u ]     u_{i+1} - 2u_i + u_{i-1}                           \n\
14:    [ ----- ]  =  ------------------------                           \n\
15:    [ dx^2  ]i              h^2                                      \n\
16:                                                                     \n\
17:     i = grid index;    h = x_{i+1}-x_i (Uniform)                    \n\
18:     0 <= i < n         h = 1.0/n                                    \n\
19:                                                                     \n\
20: Thus,                                                               \n\
21:                                                                     \n\
22:    du                                                               \n\
23:    --  = Ju;  J = (a/h^2) tridiagonal(1,-2,1)_n                     \n\
24:    dt                                                               \n\
25:                                                                     \n\
26: Implicit Runge-Kutta method:                                        \n\
27:                                                                     \n\
28:   U^(k)   = u^n + dt \\sum_i a_{ki} JU^{i}                          \n\
29:   u^{n+1} = u^n + dt \\sum_i b_i JU^{i}                             \n\
30:                                                                     \n\
31:   i = 1,...,s (s -> number of stages)                               \n\
32:                                                                     \n\
33: At each time step, we solve                                         \n\
34:                                                                     \n\
35:  [  1                                  ]     1                      \n\
36:  [ -- I \\otimes A^{-1} - J \\otimes I ] U = -- u^n \\otimes A^{-1} \n\
37:  [ dt                                  ]     dt                     \n\
38:                                                                     \n\
39:   where A is the Butcher tableaux of the implicit                   \n\
40:   Runge-Kutta method,                                               \n\
41:                                                                     \n\
42: with MATKAIJ and KSP.                                               \n\
43:                                                                     \n\
44: Available IRK Methods:                                              \n\
45:   gauss       n-stage Gauss method                                  \n\
46:                                                                     \n";

48: /*T
49:   Concepts: MATKAIJ
50:   Concepts: MAT
51:   Concepts: KSP
52: T*/

54: /*
55:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
56:   automatically includes:
57:   petscsys.h      - base PETSc routines
58:   petscvec.h      - vectors
59:   petscmat.h      - matrices
60:   petscis.h       - index sets
61:   petscviewer.h   - viewers
62:   petscpc.h       - preconditioners
63: */
64:  #include <petscksp.h>
65:  #include <petscdt.h>

67: /* define the IRK methods available */
68: #define IRKGAUSS      "gauss"

70: typedef enum {
71:   PHYSICS_DIFFUSION,
72:   PHYSICS_ADVECTION
73: } PhysicsType;
74: const char *const PhysicsTypes[] = {"DIFFUSION","ADVECTION","PhysicsType","PHYSICS_",NULL};

76: typedef struct __context__ {
77:   PetscReal     a;              /* diffusion coefficient      */
78:   PetscReal     xmin,xmax;      /* domain bounds              */
79:   PetscInt      imax;           /* number of grid points      */
80:   PetscInt      niter;          /* number of time iterations  */
81:   PetscReal     dt;             /* time step size             */
82:   PhysicsType   physics_type;
83: } UserContext;

85: static PetscErrorCode ExactSolution(Vec,void*,PetscReal);
86: static PetscErrorCode RKCreate_Gauss(PetscInt,PetscScalar**,PetscScalar**,PetscReal**);
87: static PetscErrorCode Assemble_AdvDiff(MPI_Comm,UserContext*,Mat*);

89:  #include <petsc/private/kernels/blockinvert.h>

91: int main(int argc, char **argv)
92: {
93:   PetscErrorCode    ierr;
94:   Vec               u,uex,rhs,z;
95:   UserContext       ctxt;
96:   PetscInt          nstages,is,ie,matis,matie,*ix,*ix2;
97:   PetscInt          n,i,s,t,total_its;
98:   PetscScalar       *A,*B,*At,*b,*zvals,one = 1.0;
99:   PetscReal         *c,err,time;
100:   Mat               Identity,J,TA,SC,R;
101:   KSP               ksp;
102:   PetscFunctionList IRKList = NULL;
103:   char              irktype[256] = IRKGAUSS;

105:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
106:   PetscFunctionListAdd(&IRKList,IRKGAUSS,RKCreate_Gauss);

108:   /* default value */
109:   ctxt.a       = 1.0;
110:   ctxt.xmin    = 0.0;
111:   ctxt.xmax    = 1.0;
112:   ctxt.imax    = 20;
113:   ctxt.niter   = 0;
114:   ctxt.dt      = 0.0;
115:   ctxt.physics_type = PHYSICS_DIFFUSION;

117:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"IRK options","");
118:   PetscOptionsReal("-a","diffusion coefficient","<1.0>",ctxt.a,&ctxt.a,NULL);
119:   PetscOptionsInt ("-imax","grid size","<20>",ctxt.imax,&ctxt.imax,NULL);
120:   PetscOptionsReal("-xmin","xmin","<0.0>",ctxt.xmin,&ctxt.xmin,NULL);
121:   PetscOptionsReal("-xmax","xmax","<1.0>",ctxt.xmax,&ctxt.xmax,NULL);
122:   PetscOptionsInt ("-niter","number of time steps","<0>",ctxt.niter,&ctxt.niter,NULL);
123:   PetscOptionsReal("-dt","time step size","<0.0>",ctxt.dt,&ctxt.dt,NULL);
124:   PetscOptionsFList("-irk_type","IRK method family","",IRKList,irktype,irktype,sizeof(irktype),NULL);
125:   nstages = 2;
126:   PetscOptionsInt ("-irk_nstages","Number of stages in IRK method","",nstages,&nstages,NULL);
127:   PetscOptionsEnum("-physics_type","Type of process to discretize","",PhysicsTypes,(PetscEnum)ctxt.physics_type,(PetscEnum*)&ctxt.physics_type,NULL);
128:   PetscOptionsEnd();

130:   /* allocate and initialize solution vector and exact solution */
131:   VecCreate(PETSC_COMM_WORLD,&u);
132:   VecSetSizes(u,PETSC_DECIDE,ctxt.imax);
133:   VecSetFromOptions(u);
134:   VecDuplicate(u,&uex);
135:   /* initial solution */
136:   ExactSolution(u  ,&ctxt,0.0);
137:   /* exact   solution */
138:   ExactSolution(uex,&ctxt,ctxt.dt*ctxt.niter);

140:   {                             /* Create A,b,c */
141:     PetscErrorCode (*irkcreate)(PetscInt,PetscScalar**,PetscScalar**,PetscReal**);
142:     PetscFunctionListFind(IRKList,irktype,&irkcreate);
143:     (*irkcreate)(nstages,&A,&b,&c);
144:   }
145:   {                             /* Invert A */
146:     /* PETSc does not provide a routine to calculate the inverse of a general matrix.
147:      * To get the inverse of A, we form a sequential BAIJ matrix from it, consisting of a single block with block size
148:      * equal to the dimension of A, and then use MatInvertBlockDiagonal(). */
149:     Mat               A_baij;
150:     PetscInt          idxm[1]={0},idxn[1]={0};
151:     const PetscScalar *A_inv;
152:     MatCreateSeqBAIJ(PETSC_COMM_SELF,nstages,nstages,nstages,1,NULL,&A_baij);
153:     MatSetOption(A_baij,MAT_ROW_ORIENTED,PETSC_FALSE);
154:     MatSetValuesBlocked(A_baij,1,idxm,1,idxn,A,INSERT_VALUES);
155:     MatAssemblyBegin(A_baij,MAT_FINAL_ASSEMBLY);
156:     MatAssemblyEnd(A_baij,MAT_FINAL_ASSEMBLY);
157:     MatInvertBlockDiagonal(A_baij,&A_inv);
158:     PetscMemcpy(A,A_inv,nstages*nstages*sizeof(PetscScalar));
159:     MatDestroy(&A_baij);
160:   }
161:   /* Scale (1/dt)*A^{-1} and (1/dt)*b */
162:   for (s=0; s<nstages*nstages; s++) A[s] *= 1.0/ctxt.dt;
163:   for (s=0; s<nstages; s++) b[s] *= (-ctxt.dt);

165:   /* Compute row sums At and identity B */
166:   PetscMalloc2(nstages,&At,PetscSqr(nstages),&B);
167:   for (s=0; s<nstages; s++) {
168:     At[s] = 0;
169:     for (t=0; t<nstages; t++) {
170:       At[s] += A[s+nstages*t];      /* Row sums of  */
171:       B[s+nstages*t] = 1.*(s == t); /* identity */
172:     }
173:   }

175:   /* allocate and calculate the (-J) matrix */
176:   switch (ctxt.physics_type) {
177:   case PHYSICS_ADVECTION:
178:   case PHYSICS_DIFFUSION:
179:     Assemble_AdvDiff(PETSC_COMM_WORLD,&ctxt,&J);
180:   }
181:   MatCreate(PETSC_COMM_WORLD,&Identity);
182:   MatSetType(Identity,MATAIJ);
183:   MatGetOwnershipRange(J,&matis,&matie);
184:   MatSetSizes(Identity,matie-matis,matie-matis,ctxt.imax,ctxt.imax);
185:   MatSetUp(Identity);
186:   for (i=matis; i<matie; i++) {
187:     ierr= MatSetValues(Identity,1,&i,1,&i,&one,INSERT_VALUES);
188:   }
189:   MatAssemblyBegin(Identity,MAT_FINAL_ASSEMBLY);
190:   MatAssemblyEnd  (Identity,MAT_FINAL_ASSEMBLY);

192:   /* Create the KAIJ matrix for solving the stages */
193:   MatCreateKAIJ(J,nstages,nstages,A,B,&TA);

195:   /* Create the KAIJ matrix for step completion */
196:   MatCreateKAIJ(J,1,nstages,NULL,b,&SC);

198:   /* Create the KAIJ matrix to create the R for solving the stages */
199:   MatCreateKAIJ(Identity,nstages,1,NULL,At,&R);

201:   /* Create and set options for KSP */
202:   KSPCreate(PETSC_COMM_WORLD,&ksp);
203:   KSPSetOperators(ksp,TA,TA);
204:   KSPSetFromOptions(ksp);

206:   /* Allocate work and right-hand-side vectors */
207:   VecCreate(PETSC_COMM_WORLD,&z);
208:   VecSetFromOptions(z);
209:   VecSetSizes(z,PETSC_DECIDE,ctxt.imax*nstages);
210:   VecDuplicate(z,&rhs);

212:   VecGetOwnershipRange(u,&is,&ie);
213:   PetscMalloc3(nstages,&ix,nstages,&zvals,ie-is,&ix2);
214:   /* iterate in time */
215:   for (n=0,time=0.,total_its=0; n<ctxt.niter; n++) {
216:     PetscInt its;

218:     /* compute and set the right hand side */
219:     MatMult(R,u,rhs);

221:     /* Solve the system */
222:     KSPSolve(ksp,rhs,z);
223:     KSPGetIterationNumber(ksp,&its);
224:     total_its += its;

226:     /* Update the solution */
227:     MatMultAdd(SC,z,u,u);

229:     /* time step complete */
230:     time += ctxt.dt;
231:   }
232:   PetscFree3(ix,ix2,zvals);

234:   /* Deallocate work and right-hand-side vectors */
235:   VecDestroy(&z);
236:   VecDestroy(&rhs);

238:   /* Calculate error in final solution */
239:   VecAYPX(uex,-1.0,u);
240:   VecNorm(uex,NORM_2,&err);
241:   err  = PetscSqrtReal(err*err/((PetscReal)ctxt.imax));
242:   PetscPrintf(PETSC_COMM_WORLD,"L2 norm of the numerical error = %g (time=%g)\n",(double)err,(double)time);
243:   PetscPrintf(PETSC_COMM_WORLD,"Number of time steps: %D (%D Krylov iterations)\n",ctxt.niter,total_its);

245:   /* Free up memory */
246:   KSPDestroy(&ksp);
247:   MatDestroy(&TA);
248:   MatDestroy(&SC);
249:   MatDestroy(&R);
250:   MatDestroy(&J);
251:   MatDestroy(&Identity);
252:   PetscFree3(A,b,c);
253:   PetscFree2(At,B);
254:   VecDestroy(&uex);
255:   VecDestroy(&u);
256:   PetscFunctionListDestroy(&IRKList);

258:   PetscFinalize();
259:   return ierr;
260: }

262: PetscErrorCode ExactSolution(Vec u,void *c,PetscReal t)
263: {
264:   UserContext     *ctxt = (UserContext*) c;
265:   PetscErrorCode  ierr;
266:   PetscInt        i,is,ie;
267:   PetscScalar     *uarr;
268:   PetscReal       x,dx,a=ctxt->a,pi=PETSC_PI;

271:   dx = (ctxt->xmax - ctxt->xmin)/((PetscReal) ctxt->imax);
272:   VecGetOwnershipRange(u,&is,&ie);
273:   VecGetArray(u,&uarr);
274:   for(i=is; i<ie; i++) {
275:     x          = i * dx;
276:     switch (ctxt->physics_type) {
277:     case PHYSICS_DIFFUSION:
278:       uarr[i-is] = PetscExpScalar(-4.0*pi*pi*a*t)*PetscSinScalar(2*pi*x);
279:       break;
280:     case PHYSICS_ADVECTION:
281:       uarr[i-is] = PetscSinScalar(2*pi*(x - a*t));
282:       break;
283:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[ctxt->physics_type]);
284:     }
285:   }
286:   VecRestoreArray(u,&uarr);
287:   return(0);
288: }

290: /* Arrays should be freed with PetscFree3(A,b,c) */
291: static PetscErrorCode RKCreate_Gauss(PetscInt nstages,PetscScalar **gauss_A,PetscScalar **gauss_b,PetscReal **gauss_c)
292: {
293:   PetscErrorCode    ierr;
294:   PetscScalar       *A,*G0,*G1;
295:   PetscReal         *b,*c;
296:   PetscInt          i,j;
297:   Mat               G0mat,G1mat,Amat;

300:   PetscMalloc3(PetscSqr(nstages),&A,nstages,gauss_b,nstages,&c);
301:   PetscMalloc3(nstages,&b,PetscSqr(nstages),&G0,PetscSqr(nstages),&G1);
302:   PetscDTGaussQuadrature(nstages,0.,1.,c,b);
303:   for (i=0; i<nstages; i++) (*gauss_b)[i] = b[i]; /* copy to possibly-complex array */

305:   /* A^T = G0^{-1} G1 */
306:   for (i=0; i<nstages; i++) {
307:     for (j=0; j<nstages; j++) {
308:       G0[i*nstages+j] = PetscPowRealInt(c[i],j);
309:       G1[i*nstages+j] = PetscPowRealInt(c[i],j+1)/(j+1);
310:     }
311:   }
312:   /* The arrays above are row-aligned, but we create dense matrices as the transpose */
313:   MatCreateSeqDense(PETSC_COMM_SELF,nstages,nstages,G0,&G0mat);
314:   MatCreateSeqDense(PETSC_COMM_SELF,nstages,nstages,G1,&G1mat);
315:   MatCreateSeqDense(PETSC_COMM_SELF,nstages,nstages,A,&Amat);
316:   MatLUFactor(G0mat,NULL,NULL,NULL);
317:   MatMatSolve(G0mat,G1mat,Amat);
318:   MatTranspose(Amat,MAT_INPLACE_MATRIX,&Amat);

320:   MatDestroy(&G0mat);
321:   MatDestroy(&G1mat);
322:   MatDestroy(&Amat);
323:   PetscFree3(b,G0,G1);
324:   *gauss_A = A;
325:   *gauss_c = c;
326:   return(0);
327: }

329: static PetscErrorCode Assemble_AdvDiff(MPI_Comm comm,UserContext *user,Mat *J)
330: {
332:   PetscInt       matis,matie,i;
333:   PetscReal      dx,dx2;

336:   dx = (user->xmax - user->xmin)/((PetscReal)user->imax); dx2 = dx*dx;
337:   MatCreate(comm,J);
338:   MatSetType(*J,MATAIJ);
339:   MatSetSizes(*J,PETSC_DECIDE,PETSC_DECIDE,user->imax,user->imax);
340:   MatSetUp(*J);
341:   MatGetOwnershipRange(*J,&matis,&matie);
342:   for (i=matis; i<matie; i++) {
343:     PetscScalar values[3];
344:     PetscInt    col[3];
345:     switch (user->physics_type) {
346:     case PHYSICS_DIFFUSION:
347:       values[0] = -user->a*1.0/dx2;
348:       values[1] = user->a*2.0/dx2;
349:       values[2] = -user->a*1.0/dx2;
350:       break;
351:     case PHYSICS_ADVECTION:
352:       values[0] = -user->a*.5/dx;
353:       values[1] = 0.;
354:       values[2] = user->a*.5/dx;
355:       break;
356:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[user->physics_type]);
357:     }
358:     /* periodic boundaries */
359:     if (i == 0) {
360:       col[0] = user->imax-1;
361:       col[1] = i;
362:       col[2] = i+1;
363:     } else if (i == user->imax-1) {
364:       col[0] = i-1;
365:       col[1] = i;
366:       col[2] = 0;
367:     } else {
368:       col[0] = i-1;
369:       col[1] = i;
370:       col[2] = i+1;
371:     }
372:     ierr= MatSetValues(*J,1,&i,3,col,values,INSERT_VALUES);
373:   }
374:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
375:   MatAssemblyEnd  (*J,MAT_FINAL_ASSEMBLY);
376:   return(0);
377: }

379: /*TEST
380:  test:
381:    suffix: 1
382:    args: -a 0.1 -dt .125 -niter 5 -imax 40 -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -irk_type gauss -irk_nstages 2
383:  test:
384:    suffix: 2
385:    args: -a 0.1 -dt .125 -niter 5 -imax 40 -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -irk_type gauss -irk_nstages 4 -ksp_gmres_restart 100
386:  test:
387:    suffix: 3
388:    requires: !single
389:    args: -a 1 -dt .33 -niter 3 -imax 40 -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -irk_type gauss -irk_nstages 4 -ksp_gmres_restart 100 -physics_type advection

391: TEST*/