Actual source code: ex74.c
1: static char help[] = "Solves the constant-coefficient 1D heat equation \n\
2: with an Implicit Runge-Kutta method using MatKAIJ. \n\
3: \n\
4: du d^2 u \n\
5: -- = a ----- ; 0 <= x <= 1; \n\
6: dt dx^2 \n\
7: \n\
8: with periodic boundary conditions \n\
9: \n\
10: 2nd order central discretization in space: \n\
11: \n\
12: [ d^2 u ] u_{i+1} - 2u_i + u_{i-1} \n\
13: [ ----- ] = ------------------------ \n\
14: [ dx^2 ]i h^2 \n\
15: \n\
16: i = grid index; h = x_{i+1}-x_i (Uniform) \n\
17: 0 <= i < n h = 1.0/n \n\
18: \n\
19: Thus, \n\
20: \n\
21: du \n\
22: -- = Ju; J = (a/h^2) tridiagonal(1,-2,1)_n \n\
23: dt \n\
24: \n\
25: Implicit Runge-Kutta method: \n\
26: \n\
27: U^(k) = u^n + dt \\sum_i a_{ki} JU^{i} \n\
28: u^{n+1} = u^n + dt \\sum_i b_i JU^{i} \n\
29: \n\
30: i = 1,...,s (s -> number of stages) \n\
31: \n\
32: At each time step, we solve \n\
33: \n\
34: [ 1 ] 1 \n\
35: [ -- I \\otimes A^{-1} - J \\otimes I ] U = -- u^n \\otimes A^{-1} \n\
36: [ dt ] dt \n\
37: \n\
38: where A is the Butcher tableaux of the implicit \n\
39: Runge-Kutta method, \n\
40: \n\
41: with MATKAIJ and KSP. \n\
42: \n\
43: Available IRK Methods: \n\
44: gauss n-stage Gauss method \n\
45: \n";
47: /*T
48: Concepts: MATKAIJ
49: Concepts: MAT
50: Concepts: KSP
51: T*/
53: /*
54: Include "petscksp.h" so that we can use KSP solvers. Note that this file
55: automatically includes:
56: petscsys.h - base PETSc routines
57: petscvec.h - vectors
58: petscmat.h - matrices
59: petscis.h - index sets
60: petscviewer.h - viewers
61: petscpc.h - preconditioners
62: */
63: #include <petscksp.h>
64: #include <petscdt.h>
66: /* define the IRK methods available */
67: #define IRKGAUSS "gauss"
69: typedef enum {
70: PHYSICS_DIFFUSION,
71: PHYSICS_ADVECTION
72: } PhysicsType;
73: const char *const PhysicsTypes[] = {"DIFFUSION","ADVECTION","PhysicsType","PHYSICS_",NULL};
75: typedef struct __context__ {
76: PetscReal a; /* diffusion coefficient */
77: PetscReal xmin,xmax; /* domain bounds */
78: PetscInt imax; /* number of grid points */
79: PetscInt niter; /* number of time iterations */
80: PetscReal dt; /* time step size */
81: PhysicsType physics_type;
82: } UserContext;
84: static PetscErrorCode ExactSolution(Vec,void*,PetscReal);
85: static PetscErrorCode RKCreate_Gauss(PetscInt,PetscScalar**,PetscScalar**,PetscReal**);
86: static PetscErrorCode Assemble_AdvDiff(MPI_Comm,UserContext*,Mat*);
88: #include <petsc/private/kernels/blockinvert.h>
90: int main(int argc, char **argv)
91: {
92: PetscErrorCode ierr;
93: Vec u,uex,rhs,z;
94: UserContext ctxt;
95: PetscInt nstages,is,ie,matis,matie,*ix,*ix2;
96: PetscInt n,i,s,t,total_its;
97: PetscScalar *A,*B,*At,*b,*zvals,one = 1.0;
98: PetscReal *c,err,time;
99: Mat Identity,J,TA,SC,R;
100: KSP ksp;
101: PetscFunctionList IRKList = NULL;
102: char irktype[256] = IRKGAUSS;
104: PetscInitialize(&argc,&argv,(char*)0,help);
105: PetscFunctionListAdd(&IRKList,IRKGAUSS,RKCreate_Gauss);
107: /* default value */
108: ctxt.a = 1.0;
109: ctxt.xmin = 0.0;
110: ctxt.xmax = 1.0;
111: ctxt.imax = 20;
112: ctxt.niter = 0;
113: ctxt.dt = 0.0;
114: ctxt.physics_type = PHYSICS_DIFFUSION;
116: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"IRK options","");
117: PetscOptionsReal("-a","diffusion coefficient","<1.0>",ctxt.a,&ctxt.a,NULL);
118: PetscOptionsInt ("-imax","grid size","<20>",ctxt.imax,&ctxt.imax,NULL);
119: PetscOptionsReal("-xmin","xmin","<0.0>",ctxt.xmin,&ctxt.xmin,NULL);
120: PetscOptionsReal("-xmax","xmax","<1.0>",ctxt.xmax,&ctxt.xmax,NULL);
121: PetscOptionsInt ("-niter","number of time steps","<0>",ctxt.niter,&ctxt.niter,NULL);
122: PetscOptionsReal("-dt","time step size","<0.0>",ctxt.dt,&ctxt.dt,NULL);
123: PetscOptionsFList("-irk_type","IRK method family","",IRKList,irktype,irktype,sizeof(irktype),NULL);
124: nstages = 2;
125: PetscOptionsInt ("-irk_nstages","Number of stages in IRK method","",nstages,&nstages,NULL);
126: PetscOptionsEnum("-physics_type","Type of process to discretize","",PhysicsTypes,(PetscEnum)ctxt.physics_type,(PetscEnum*)&ctxt.physics_type,NULL);
127: PetscOptionsEnd();
129: /* allocate and initialize solution vector and exact solution */
130: VecCreate(PETSC_COMM_WORLD,&u);
131: VecSetSizes(u,PETSC_DECIDE,ctxt.imax);
132: VecSetFromOptions(u);
133: VecDuplicate(u,&uex);
134: /* initial solution */
135: ExactSolution(u ,&ctxt,0.0);
136: /* exact solution */
137: ExactSolution(uex,&ctxt,ctxt.dt*ctxt.niter);
139: { /* Create A,b,c */
140: PetscErrorCode (*irkcreate)(PetscInt,PetscScalar**,PetscScalar**,PetscReal**);
141: PetscFunctionListFind(IRKList,irktype,&irkcreate);
142: (*irkcreate)(nstages,&A,&b,&c);
143: }
144: { /* Invert A */
145: /* PETSc does not provide a routine to calculate the inverse of a general matrix.
146: * To get the inverse of A, we form a sequential BAIJ matrix from it, consisting of a single block with block size
147: * equal to the dimension of A, and then use MatInvertBlockDiagonal(). */
148: Mat A_baij;
149: PetscInt idxm[1]={0},idxn[1]={0};
150: const PetscScalar *A_inv;
151: MatCreateSeqBAIJ(PETSC_COMM_SELF,nstages,nstages,nstages,1,NULL,&A_baij);
152: MatSetOption(A_baij,MAT_ROW_ORIENTED,PETSC_FALSE);
153: MatSetValuesBlocked(A_baij,1,idxm,1,idxn,A,INSERT_VALUES);
154: MatAssemblyBegin(A_baij,MAT_FINAL_ASSEMBLY);
155: MatAssemblyEnd(A_baij,MAT_FINAL_ASSEMBLY);
156: MatInvertBlockDiagonal(A_baij,&A_inv);
157: PetscMemcpy(A,A_inv,nstages*nstages*sizeof(PetscScalar));
158: MatDestroy(&A_baij);
159: }
160: /* Scale (1/dt)*A^{-1} and (1/dt)*b */
161: for (s=0; s<nstages*nstages; s++) A[s] *= 1.0/ctxt.dt;
162: for (s=0; s<nstages; s++) b[s] *= (-ctxt.dt);
164: /* Compute row sums At and identity B */
165: PetscMalloc2(nstages,&At,PetscSqr(nstages),&B);
166: for (s=0; s<nstages; s++) {
167: At[s] = 0;
168: for (t=0; t<nstages; t++) {
169: At[s] += A[s+nstages*t]; /* Row sums of */
170: B[s+nstages*t] = 1.*(s == t); /* identity */
171: }
172: }
174: /* allocate and calculate the (-J) matrix */
175: switch (ctxt.physics_type) {
176: case PHYSICS_ADVECTION:
177: case PHYSICS_DIFFUSION:
178: Assemble_AdvDiff(PETSC_COMM_WORLD,&ctxt,&J);
179: }
180: MatCreate(PETSC_COMM_WORLD,&Identity);
181: MatSetType(Identity,MATAIJ);
182: MatGetOwnershipRange(J,&matis,&matie);
183: MatSetSizes(Identity,matie-matis,matie-matis,ctxt.imax,ctxt.imax);
184: MatSetUp(Identity);
185: for (i=matis; i<matie; i++) {
186: MatSetValues(Identity,1,&i,1,&i,&one,INSERT_VALUES);
187: }
188: MatAssemblyBegin(Identity,MAT_FINAL_ASSEMBLY);
189: MatAssemblyEnd (Identity,MAT_FINAL_ASSEMBLY);
191: /* Create the KAIJ matrix for solving the stages */
192: MatCreateKAIJ(J,nstages,nstages,A,B,&TA);
194: /* Create the KAIJ matrix for step completion */
195: MatCreateKAIJ(J,1,nstages,NULL,b,&SC);
197: /* Create the KAIJ matrix to create the R for solving the stages */
198: MatCreateKAIJ(Identity,nstages,1,NULL,At,&R);
200: /* Create and set options for KSP */
201: KSPCreate(PETSC_COMM_WORLD,&ksp);
202: KSPSetOperators(ksp,TA,TA);
203: KSPSetFromOptions(ksp);
205: /* Allocate work and right-hand-side vectors */
206: VecCreate(PETSC_COMM_WORLD,&z);
207: VecSetFromOptions(z);
208: VecSetSizes(z,PETSC_DECIDE,ctxt.imax*nstages);
209: VecDuplicate(z,&rhs);
211: VecGetOwnershipRange(u,&is,&ie);
212: PetscMalloc3(nstages,&ix,nstages,&zvals,ie-is,&ix2);
213: /* iterate in time */
214: for (n=0,time=0.,total_its=0; n<ctxt.niter; n++) {
215: PetscInt its;
217: /* compute and set the right hand side */
218: MatMult(R,u,rhs);
220: /* Solve the system */
221: KSPSolve(ksp,rhs,z);
222: KSPGetIterationNumber(ksp,&its);
223: total_its += its;
225: /* Update the solution */
226: MatMultAdd(SC,z,u,u);
228: /* time step complete */
229: time += ctxt.dt;
230: }
231: PetscFree3(ix,ix2,zvals);
233: /* Deallocate work and right-hand-side vectors */
234: VecDestroy(&z);
235: VecDestroy(&rhs);
237: /* Calculate error in final solution */
238: VecAYPX(uex,-1.0,u);
239: VecNorm(uex,NORM_2,&err);
240: err = PetscSqrtReal(err*err/((PetscReal)ctxt.imax));
241: PetscPrintf(PETSC_COMM_WORLD,"L2 norm of the numerical error = %g (time=%g)\n",(double)err,(double)time);
242: PetscPrintf(PETSC_COMM_WORLD,"Number of time steps: %D (%D Krylov iterations)\n",ctxt.niter,total_its);
244: /* Free up memory */
245: KSPDestroy(&ksp);
246: MatDestroy(&TA);
247: MatDestroy(&SC);
248: MatDestroy(&R);
249: MatDestroy(&J);
250: MatDestroy(&Identity);
251: PetscFree3(A,b,c);
252: PetscFree2(At,B);
253: VecDestroy(&uex);
254: VecDestroy(&u);
255: PetscFunctionListDestroy(&IRKList);
257: PetscFinalize();
258: return 0;
259: }
261: PetscErrorCode ExactSolution(Vec u,void *c,PetscReal t)
262: {
263: UserContext *ctxt = (UserContext*) c;
264: PetscInt i,is,ie;
265: PetscScalar *uarr;
266: PetscReal x,dx,a=ctxt->a,pi=PETSC_PI;
268: dx = (ctxt->xmax - ctxt->xmin)/((PetscReal) ctxt->imax);
269: VecGetOwnershipRange(u,&is,&ie);
270: VecGetArray(u,&uarr);
271: for (i=is; i<ie; i++) {
272: x = i * dx;
273: switch (ctxt->physics_type) {
274: case PHYSICS_DIFFUSION:
275: uarr[i-is] = PetscExpScalar(-4.0*pi*pi*a*t)*PetscSinScalar(2*pi*x);
276: break;
277: case PHYSICS_ADVECTION:
278: uarr[i-is] = PetscSinScalar(2*pi*(x - a*t));
279: break;
280: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[ctxt->physics_type]);
281: }
282: }
283: VecRestoreArray(u,&uarr);
284: return 0;
285: }
287: /* Arrays should be freed with PetscFree3(A,b,c) */
288: static PetscErrorCode RKCreate_Gauss(PetscInt nstages,PetscScalar **gauss_A,PetscScalar **gauss_b,PetscReal **gauss_c)
289: {
290: PetscScalar *A,*G0,*G1;
291: PetscReal *b,*c;
292: PetscInt i,j;
293: Mat G0mat,G1mat,Amat;
295: PetscMalloc3(PetscSqr(nstages),&A,nstages,gauss_b,nstages,&c);
296: PetscMalloc3(nstages,&b,PetscSqr(nstages),&G0,PetscSqr(nstages),&G1);
297: PetscDTGaussQuadrature(nstages,0.,1.,c,b);
298: for (i=0; i<nstages; i++) (*gauss_b)[i] = b[i]; /* copy to possibly-complex array */
300: /* A^T = G0^{-1} G1 */
301: for (i=0; i<nstages; i++) {
302: for (j=0; j<nstages; j++) {
303: G0[i*nstages+j] = PetscPowRealInt(c[i],j);
304: G1[i*nstages+j] = PetscPowRealInt(c[i],j+1)/(j+1);
305: }
306: }
307: /* The arrays above are row-aligned, but we create dense matrices as the transpose */
308: MatCreateSeqDense(PETSC_COMM_SELF,nstages,nstages,G0,&G0mat);
309: MatCreateSeqDense(PETSC_COMM_SELF,nstages,nstages,G1,&G1mat);
310: MatCreateSeqDense(PETSC_COMM_SELF,nstages,nstages,A,&Amat);
311: MatLUFactor(G0mat,NULL,NULL,NULL);
312: MatMatSolve(G0mat,G1mat,Amat);
313: MatTranspose(Amat,MAT_INPLACE_MATRIX,&Amat);
315: MatDestroy(&G0mat);
316: MatDestroy(&G1mat);
317: MatDestroy(&Amat);
318: PetscFree3(b,G0,G1);
319: *gauss_A = A;
320: *gauss_c = c;
321: return 0;
322: }
324: static PetscErrorCode Assemble_AdvDiff(MPI_Comm comm,UserContext *user,Mat *J)
325: {
326: PetscInt matis,matie,i;
327: PetscReal dx,dx2;
329: dx = (user->xmax - user->xmin)/((PetscReal)user->imax); dx2 = dx*dx;
330: MatCreate(comm,J);
331: MatSetType(*J,MATAIJ);
332: MatSetSizes(*J,PETSC_DECIDE,PETSC_DECIDE,user->imax,user->imax);
333: MatSetUp(*J);
334: MatGetOwnershipRange(*J,&matis,&matie);
335: for (i=matis; i<matie; i++) {
336: PetscScalar values[3];
337: PetscInt col[3];
338: switch (user->physics_type) {
339: case PHYSICS_DIFFUSION:
340: values[0] = -user->a*1.0/dx2;
341: values[1] = user->a*2.0/dx2;
342: values[2] = -user->a*1.0/dx2;
343: break;
344: case PHYSICS_ADVECTION:
345: values[0] = -user->a*.5/dx;
346: values[1] = 0.;
347: values[2] = user->a*.5/dx;
348: break;
349: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[user->physics_type]);
350: }
351: /* periodic boundaries */
352: if (i == 0) {
353: col[0] = user->imax-1;
354: col[1] = i;
355: col[2] = i+1;
356: } else if (i == user->imax-1) {
357: col[0] = i-1;
358: col[1] = i;
359: col[2] = 0;
360: } else {
361: col[0] = i-1;
362: col[1] = i;
363: col[2] = i+1;
364: }
365: MatSetValues(*J,1,&i,3,col,values,INSERT_VALUES);
366: }
367: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
368: MatAssemblyEnd (*J,MAT_FINAL_ASSEMBLY);
369: return 0;
370: }
372: /*TEST
373: testset:
374: suffix: 1
375: args: -a 0.1 -dt .125 -niter 5 -imax 40 -ksp_monitor_short -pc_type pbjacobi -irk_type gauss -irk_nstages 2
376: test:
377: args: -ksp_atol 1e-6
378: test:
379: requires: hpddm !single
380: suffix: hpddm
381: output_file: output/ex74_1.out
382: args: -ksp_atol 1e-6 -ksp_type hpddm
383: test:
384: requires: hpddm
385: suffix: hpddm_gcrodr
386: output_file: output/ex74_1_hpddm.out
387: args: -ksp_atol 1e-4 -ksp_view_final_residual -ksp_type hpddm -ksp_hpddm_type gcrodr -ksp_hpddm_recycle 2
388: test:
389: suffix: 2
390: 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
391: testset:
392: suffix: 3
393: requires: !single
394: 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
395: test:
396: args:
397: test:
398: requires: hpddm
399: suffix: hpddm
400: output_file: output/ex74_3.out
401: args: -ksp_type hpddm
402: test:
403: requires: hpddm
404: suffix: hpddm_gcrodr
405: output_file: output/ex74_3_hpddm.out
406: args: -ksp_view_final_residual -ksp_type hpddm -ksp_hpddm_type gcrodr -ksp_hpddm_recycle 5
408: TEST*/