Actual source code: ex28.c

petsc-3.3-p7 2013-05-11
```  1: static const char help[] = "1D multiphysics prototype with analytic Jacobians to solve individual problems and a coupled problem.\n\n";

3: /* Solve a PDE coupled to an algebraic system in 1D
4:  *
5:  * PDE (U):
6:  *     -(k u_x)_x = 1 on (0,1), subject to u(0) = 0, u(1) = 1
7:  * Algebraic (K):
8:  *     exp(k-1) + k = u + 1/(1/(1+u) + 1/(1+u_x^2))
9:  *
10:  * The discretization places k at staggered points, and a separate DMDA is used for each "physics".
11:  *
12:  * This example is a prototype for coupling in multi-physics problems, therefore residual evaluation and assembly for
13:  * each problem (referred to as U and K) are written separately.  This permits the same "physics" code to be used for
14:  * solving each uncoupled problem as well as the coupled system.  In particular, run with -problem_type 0 to solve only
15:  * problem U (with K fixed), -problem_type 1 to solve only K (with U fixed), and -problem_type 2 to solve both at once.
16:  *
17:  * In all cases, a fully-assembled analytic Jacobian is available, so the systems can be solved with a direct solve or
18:  * any other standard method.  Additionally, by running with
19:  *
20:  *   -pack_dm_mat_type nest
21:  *
22:  * The same code assembles a coupled matrix where each block is stored separately, which allows the use of PCFieldSplit
23:  * without copying values to extract submatrices.
24:  */

26: #include <petscsnes.h>
27: #include <petscdmda.h>
28: #include <petscdmcomposite.h>

30: typedef struct _UserCtx *User;
31: struct _UserCtx {
32:   PetscInt ptype;
33:   DM       pack;
34:   Vec      Uloc,Kloc;
35: };

39: static PetscErrorCode FormFunctionLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])
40: {
41:   PetscReal hx = 1./info->mx;
42:   PetscInt i;

45:   for (i=info->xs; i<info->xs+info->xm; i++) {
46:     if (i == 0) f[i] = 1./hx*u[i];
47:     else if (i == info->mx-1) f[i] = 1./hx*(u[i] - 1.0);
48:     else f[i] = hx*((k[i-1]*(u[i]-u[i-1]) - k[i]*(u[i+1]-u[i]))/(hx*hx) - 1.0);
49:   }
50:   return(0);
51: }

55: static PetscErrorCode FormFunctionLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])
56: {
57:   PetscReal hx = 1./info->mx;
58:   PetscInt  i;

61:   for (i=info->xs; i<info->xs+info->xm; i++) {
62:     const PetscScalar
63:       ubar = 0.5*(u[i+1]+u[i]),
66:       w = 1./(1.+ubar) + 1./g;
67:     f[i] = hx*(PetscExpScalar(k[i]-1.0) + k[i] - 1./w);
68:   }
69:   return(0);
70: }

74: static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx)
75: {
76:   User              user = (User)ctx;
77:   DM                dau,dak;
78:   DMDALocalInfo     infou,infok;
79:   PetscScalar       *u,*k;
80:   PetscScalar       *fu,*fk;
81:   PetscErrorCode    ierr;
82:   Vec               Uloc,Kloc,Fu,Fk;

85:   DMCompositeGetEntries(user->pack,&dau,&dak);
86:   DMDAGetLocalInfo(dau,&infou);
87:   DMDAGetLocalInfo(dak,&infok);
88:   DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
89:   switch (user->ptype) {
90:   case 0:
91:     DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
92:     DMGlobalToLocalEnd  (dau,X,INSERT_VALUES,Uloc);
93:     DMDAVecGetArray(dau,Uloc,&u);
94:     DMDAVecGetArray(dak,user->Kloc,&k);
95:     DMDAVecGetArray(dau,F,&fu);
96:     FormFunctionLocal_U(user,&infou,u,k,fu);
97:     DMDAVecRestoreArray(dau,F,&fu);
98:     DMDAVecRestoreArray(dau,Uloc,&u);
99:     DMDAVecRestoreArray(dak,user->Kloc,&k);
100:     break;
101:   case 1:
102:     DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
103:     DMGlobalToLocalEnd  (dak,X,INSERT_VALUES,Kloc);
104:     DMDAVecGetArray(dau,user->Uloc,&u);
105:     DMDAVecGetArray(dak,Kloc,&k);
106:     DMDAVecGetArray(dak,F,&fk);
107:     FormFunctionLocal_K(user,&infok,u,k,fk);
108:     DMDAVecRestoreArray(dak,F,&fk);
109:     DMDAVecRestoreArray(dau,user->Uloc,&u);
110:     DMDAVecRestoreArray(dak,Kloc,&k);
111:     break;
112:   case 2:
113:     DMCompositeScatter(user->pack,X,Uloc,Kloc);
114:     DMDAVecGetArray(dau,Uloc,&u);
115:     DMDAVecGetArray(dak,Kloc,&k);
116:     DMCompositeGetAccess(user->pack,F,&Fu,&Fk);
117:     DMDAVecGetArray(dau,Fu,&fu);
118:     DMDAVecGetArray(dak,Fk,&fk);
119:     FormFunctionLocal_U(user,&infou,u,k,fu);
120:     FormFunctionLocal_K(user,&infok,u,k,fk);
121:     DMDAVecRestoreArray(dau,Fu,&fu);
122:     DMDAVecRestoreArray(dak,Fk,&fk);
123:     DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk);
124:     DMDAVecRestoreArray(dau,Uloc,&u);
125:     DMDAVecRestoreArray(dak,Kloc,&k);
126:     break;
127:   }
128:   DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
129:   return(0);
130: }

134: static PetscErrorCode FormJacobianLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Buu)
135: {
136:   PetscReal      hx = 1./info->mx;
138:   PetscInt       i;

141:   for (i=info->xs; i<info->xs+info->xm; i++) {
142:     PetscInt row = i-info->gxs,cols[] = {row-1,row,row+1};
143:     PetscScalar val = 1./hx;
144:     if (i == 0) {MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);}
145:     else if (i == info->mx-1) {MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);}
146:     else {
147:       PetscScalar vals[] = {-k[i-1]/hx,(k[i-1]+k[i])/hx,-k[i]/hx};
148:       MatSetValuesLocal(Buu,1,&row,3,cols,vals,INSERT_VALUES);
149:     }
150:   }
151:   return(0);
152: }

156: static PetscErrorCode FormJacobianLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Bkk)
157: {
158:   PetscReal      hx = 1./info->mx;
160:   PetscInt       i;

163:   for (i=info->xs; i<info->xs+info->xm; i++) {
164:     PetscInt row = i-info->gxs;
165:     PetscScalar vals[] = {hx*(PetscExpScalar(k[i]-1.)+1.)};
166:     MatSetValuesLocal(Bkk,1,&row,1,&row,vals,INSERT_VALUES);
167:   }
168:   return(0);
169: }

173: static PetscErrorCode FormJacobianLocal_UK(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Buk)
174: {
175:   PetscReal hx = 1./info->mx;
177:   PetscInt       i;
178:   PetscInt row,cols[2];
179:   PetscScalar vals[2];

182:   if (!Buk) return(0); /* Not assembling this block */
183:   for (i=info->xs; i<info->xs+info->xm; i++) {
184:     if (i == 0 || i == info->mx-1) continue;
185:     row = i-info->gxs;
186:     cols[0] = i-1-infok->gxs;  vals[0] = (u[i]-u[i-1])/hx;
187:     cols[1] = i-infok->gxs;    vals[1] = (u[i]-u[i+1])/hx;
188:     MatSetValuesLocal(Buk,1,&row,2,cols,vals,INSERT_VALUES);
189:   }
190:   return(0);
191: }

195: static PetscErrorCode FormJacobianLocal_KU(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Bku)
196: {
198:   PetscInt       i;
199:   PetscReal      hx = 1./(info->mx-1);

202:   if (!Bku) return(0); /* Not assembling this block */
203:   for (i=infok->xs; i<infok->xs+infok->xm; i++) {
204:     PetscInt row = i-infok->gxs,cols[2];
205:     PetscScalar vals[2];
206:     const PetscScalar
207:       ubar     = 0.5*(u[i]+u[i+1]),
208:       ubar_L   = 0.5,
209:       ubar_R   = 0.5,
213:       g        = 1. + PetscSqr(gradu),
215:       w        = 1./(1.+ubar) + 1./g,
216:       w_ubar   = -1./PetscSqr(1.+ubar),
218:       iw       = 1./w,
219:       iw_ubar  = -w_ubar * PetscSqr(iw),
223:     MatSetValuesLocal(Bku,1,&row,2,cols,vals,INSERT_VALUES);
224:   }
225:   return(0);
226: }

230: static PetscErrorCode FormJacobian_All(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *mstr,void *ctx)
231: {
232:   User              user = (User)ctx;
233:   DM                dau,dak;
234:   DMDALocalInfo     infou,infok;
235:   PetscScalar       *u,*k;
236:   PetscErrorCode    ierr;
237:   Vec               Uloc,Kloc;

240:   DMCompositeGetEntries(user->pack,&dau,&dak);
241:   DMDAGetLocalInfo(dau,&infou);
242:   DMDAGetLocalInfo(dak,&infok);
243:   DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
244:   switch (user->ptype) {
245:   case 0:
246:     DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
247:     DMGlobalToLocalEnd  (dau,X,INSERT_VALUES,Uloc);
248:     DMDAVecGetArray(dau,Uloc,&u);
249:     DMDAVecGetArray(dak,user->Kloc,&k);
250:     FormJacobianLocal_U(user,&infou,u,k,*B);
251:     DMDAVecRestoreArray(dau,Uloc,&u);
252:     DMDAVecRestoreArray(dak,user->Kloc,&k);
253:     break;
254:   case 1:
255:     DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
256:     DMGlobalToLocalEnd  (dak,X,INSERT_VALUES,Kloc);
257:     DMDAVecGetArray(dau,user->Uloc,&u);
258:     DMDAVecGetArray(dak,Kloc,&k);
259:     FormJacobianLocal_K(user,&infok,u,k,*B);
260:     DMDAVecRestoreArray(dau,user->Uloc,&u);
261:     DMDAVecRestoreArray(dak,Kloc,&k);
262:     break;
263:   case 2: {
264:     Mat Buu,Buk,Bku,Bkk;
265:     IS  *is;
266:     DMCompositeScatter(user->pack,X,Uloc,Kloc);
267:     DMDAVecGetArray(dau,Uloc,&u);
268:     DMDAVecGetArray(dak,Kloc,&k);
269:     DMCompositeGetLocalISs(user->pack,&is);
270:     MatGetLocalSubMatrix(*B,is[0],is[0],&Buu);
271:     MatGetLocalSubMatrix(*B,is[0],is[1],&Buk);
272:     MatGetLocalSubMatrix(*B,is[1],is[0],&Bku);
273:     MatGetLocalSubMatrix(*B,is[1],is[1],&Bkk);
274:     FormJacobianLocal_U(user,&infou,u,k,Buu);
275:     FormJacobianLocal_UK(user,&infou,&infok,u,k,Buk);
276:     FormJacobianLocal_KU(user,&infou,&infok,u,k,Bku);
277:     FormJacobianLocal_K(user,&infok,u,k,Bkk);
278:     MatRestoreLocalSubMatrix(*B,is[0],is[0],&Buu);
279:     MatRestoreLocalSubMatrix(*B,is[0],is[1],&Buk);
280:     MatRestoreLocalSubMatrix(*B,is[1],is[0],&Bku);
281:     MatRestoreLocalSubMatrix(*B,is[1],is[1],&Bkk);
282:     DMDAVecRestoreArray(dau,Uloc,&u);
283:     DMDAVecRestoreArray(dak,Kloc,&k);

285:     ISDestroy(&is[0]);
286:     ISDestroy(&is[1]);
287:     PetscFree(is);
288:   } break;
289:   }
290:   DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
291:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
292:   MatAssemblyEnd  (*B,MAT_FINAL_ASSEMBLY);
293:   if (*J != *B) {
294:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
295:     MatAssemblyEnd  (*J,MAT_FINAL_ASSEMBLY);
296:   }
297:   *mstr = DIFFERENT_NONZERO_PATTERN;
298:   return(0);
299: }

303: static PetscErrorCode FormInitial_Coupled(User user,Vec X)
304: {
306:   DM             dau,dak;
307:   DMDALocalInfo  infou,infok;
308:   Vec            Xu,Xk;
309:   PetscScalar    *u,*k,hx;
310:   PetscInt       i;

313:   DMCompositeGetEntries(user->pack,&dau,&dak);
314:   DMCompositeGetAccess(user->pack,X,&Xu,&Xk);
315:   DMDAVecGetArray(dau,Xu,&u);
316:   DMDAVecGetArray(dak,Xk,&k);
317:   DMDAGetLocalInfo(dau,&infou);
318:   DMDAGetLocalInfo(dak,&infok);
319:   hx = 1./(infok.mx);
320:   for (i=infou.xs; i<infou.xs+infou.xm; i++) u[i] = (PetscScalar)i*hx * (1.-(PetscScalar)i*hx);
321:   for (i=infok.xs; i<infok.xs+infok.xm; i++) k[i] = 1.0 + 0.5*(PetscScalar)sin((double)2*PETSC_PI*i*hx);
322:   DMDAVecRestoreArray(dau,Xu,&u);
323:   DMDAVecRestoreArray(dak,Xk,&k);
324:   DMCompositeRestoreAccess(user->pack,X,&Xu,&Xk);
325:   DMCompositeScatter(user->pack,X,user->Uloc,user->Kloc);
326:   return(0);
327: }

331: int main(int argc, char *argv[])
332: {
334:   DM             dau,dak,pack;
335:   const PetscInt *lxu;
336:   PetscInt       *lxk,m,nprocs;
337:   User           user;
338:   SNES           snes;
339:   Vec            X,F,Xu,Xk,Fu,Fk;
340:   Mat            B;
341:   IS             *isg;
342:   PetscBool      view_draw,pass_dm;

344:   PetscInitialize(&argc,&argv,0,help);
345:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-10,1,1,PETSC_NULL,&dau);
346:   DMSetOptionsPrefix(dau,"u_");
347:   DMSetFromOptions(dau);
348:   DMDAGetOwnershipRanges(dau,&lxu,0,0);
349:   DMDAGetInfo(dau,0, &m,0,0, &nprocs,0,0, 0,0,0,0,0,0);
350:   PetscMalloc(nprocs*sizeof(*lxk),&lxk);
351:   PetscMemcpy(lxk,lxu,nprocs*sizeof(*lxk));
352:   lxk[0]--;
353:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,m-1,1,1,lxk,&dak);
354:   DMSetOptionsPrefix(dak,"k_");
355:   DMSetFromOptions(dau);
356:   PetscFree(lxk);

358:   DMCompositeCreate(PETSC_COMM_WORLD,&pack);
359:   DMSetOptionsPrefix(pack,"pack_");
362:   DMDASetFieldName(dau,0,"u");
363:   DMDASetFieldName(dak,0,"k");
364:   DMSetFromOptions(pack);

366:   DMCreateGlobalVector(pack,&X);
367:   VecDuplicate(X,&F);

369:   PetscNew(struct _UserCtx,&user);
370:   user->pack = pack;
371:   DMCompositeGetGlobalISs(pack,&isg);
372:   DMCompositeGetLocalVectors(pack,&user->Uloc,&user->Kloc);
373:   DMCompositeScatter(pack,X,user->Uloc,user->Kloc);

375:   PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Coupled problem options","SNES");
376:   {
377:     user->ptype = 0; view_draw = PETSC_FALSE; pass_dm = PETSC_TRUE;
378:     PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,0);
379:     PetscOptionsBool("-view_draw","Draw the final coupled solution regardless of whether only one physics was solved",0,view_draw,&view_draw,0);
380:     PetscOptionsBool("-pass_dm","Pass the packed DM to SNES to use when determining splits and forward into splits",0,pass_dm,&pass_dm,0);
381:   }
382:   PetscOptionsEnd();

384:   FormInitial_Coupled(user,X);

386:   SNESCreate(PETSC_COMM_WORLD,&snes);
387:   switch (user->ptype) {
388:   case 0:
389:     DMCompositeGetAccess(pack,X,&Xu,0);
390:     DMCompositeGetAccess(pack,F,&Fu,0);
391:     DMCreateMatrix(dau,PETSC_NULL,&B);
392:     SNESSetFunction(snes,Fu,FormFunction_All,user);
393:     SNESSetJacobian(snes,B,B,FormJacobian_All,user);
394:     SNESSetFromOptions(snes);
395:     SNESSetDM(snes,dau);
396:     SNESSolve(snes,PETSC_NULL,Xu);
397:     DMCompositeRestoreAccess(pack,X,&Xu,0);
398:     DMCompositeRestoreAccess(pack,F,&Fu,0);
399:     break;
400:   case 1:
401:     DMCompositeGetAccess(pack,X,0,&Xk);
402:     DMCompositeGetAccess(pack,F,0,&Fk);
403:     DMCreateMatrix(dak,PETSC_NULL,&B);
404:     SNESSetFunction(snes,Fk,FormFunction_All,user);
405:     SNESSetJacobian(snes,B,B,FormJacobian_All,user);
406:     SNESSetFromOptions(snes);
407:     SNESSetDM(snes,dak);
408:     SNESSolve(snes,PETSC_NULL,Xk);
409:     DMCompositeRestoreAccess(pack,X,0,&Xk);
410:     DMCompositeRestoreAccess(pack,F,0,&Fk);
411:     break;
412:   case 2:
413:     DMCreateMatrix(pack,PETSC_NULL,&B);
414:     /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */
415:     MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
416:     MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
417:     SNESSetFunction(snes,F,FormFunction_All,user);
418:     SNESSetJacobian(snes,B,B,FormJacobian_All,user);
419:     SNESSetFromOptions(snes);
420:     if (!pass_dm) {             /* Manually provide index sets and names for the splits */
421:       KSP ksp;
422:       PC  pc;
423:       SNESGetKSP(snes,&ksp);
424:       KSPGetPC(ksp,&pc);
425:       PCFieldSplitSetIS(pc,"u",isg[0]);
426:       PCFieldSplitSetIS(pc,"k",isg[1]);
427:     } else {
428:       /* The same names come from the options prefix for dau and dak. This option can support geometric multigrid inside
429:        * of splits, but it requires using a DM (perhaps your own implementation). */
430:       SNESSetDM(snes,pack);
431:     }
432:     SNESSolve(snes,PETSC_NULL,X);
433:     break;
434:   }
435:   if (view_draw) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
436:   if (0) {
437:     PetscInt col = 0;
438:     PetscBool mult_dup = PETSC_FALSE,view_dup = PETSC_FALSE;
439:     Mat D;
440:     Vec Y;

442:     PetscOptionsGetInt(0,"-col",&col,0);
443:     PetscOptionsGetBool(0,"-mult_dup",&mult_dup,0);
444:     PetscOptionsGetBool(0,"-view_dup",&view_dup,0);

446:     VecDuplicate(X,&Y);
447:     /* MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); */
448:     /* MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); */
449:     MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&D);
450:     VecZeroEntries(X);
451:     VecSetValue(X,col,1.0,INSERT_VALUES);
452:     VecAssemblyBegin(X);
453:     VecAssemblyEnd(X);
454:     MatMult(mult_dup?D:B,X,Y);
455:     MatView(view_dup?D:B,PETSC_VIEWER_STDOUT_WORLD);
456:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
457:     VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
458:     MatDestroy(&D);
459:     VecDestroy(&Y);
460:   }

462:   DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc);
463:   PetscFree(user);

465:   ISDestroy(&isg[0]);
466:   ISDestroy(&isg[1]);
467:   PetscFree(isg);
468:   VecDestroy(&X);
469:   VecDestroy(&F);
470:   MatDestroy(&B);
471:   DMDestroy(&dau);
472:   DMDestroy(&dak);
473:   DMDestroy(&pack);
474:   SNESDestroy(&snes);
475:   PetscFinalize();
476:   return 0;
477: }
```