Actual source code: ex28.c

petsc-3.4.4 2014-03-13
  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]),
 64:       gradu = (u[i+1]-u[i])/hx,
 65:       g     = 1. + gradu*gradu,
 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;
 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) {
145:       MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);
146:     } else if (i == info->mx-1) {
147:       MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);
148:     } else {
149:       PetscScalar vals[] = {-k[i-1]/hx,(k[i-1]+k[i])/hx,-k[i]/hx};
150:       MatSetValuesLocal(Buu,1,&row,3,cols,vals,INSERT_VALUES);
151:     }
152:   }
153:   return(0);
154: }

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

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

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

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

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

204:   if (!Bku) return(0); /* Not assembling this block */
205:   for (i=infok->xs; i<infok->xs+infok->xm; i++) {
206:     PetscInt    row = i-infok->gxs,cols[2];
207:     PetscScalar vals[2];
208:     const PetscScalar
209:       ubar     = 0.5*(u[i]+u[i+1]),
210:       ubar_L   = 0.5,
211:       ubar_R   = 0.5,
212:       gradu    = (u[i+1]-u[i])/hx,
213:       gradu_L  = -1./hx,
214:       gradu_R  = 1./hx,
215:       g        = 1. + PetscSqr(gradu),
216:       g_gradu  = 2.*gradu,
217:       w        = 1./(1.+ubar) + 1./g,
218:       w_ubar   = -1./PetscSqr(1.+ubar),
219:       w_gradu  = -g_gradu/PetscSqr(g),
220:       iw       = 1./w,
221:       iw_ubar  = -w_ubar * PetscSqr(iw),
222:       iw_gradu = -w_gradu * PetscSqr(iw);
223:     cols[0] = i-info->gxs;         vals[0] = -hx*(iw_ubar*ubar_L + iw_gradu*gradu_L);
224:     cols[1] = i+1-info->gxs;       vals[1] = -hx*(iw_ubar*ubar_R + iw_gradu*gradu_R);
225:     MatSetValuesLocal(Bku,1,&row,2,cols,vals,INSERT_VALUES);
226:   }
227:   return(0);
228: }

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

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

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

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

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

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

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

360:   DMCompositeCreate(PETSC_COMM_WORLD,&pack);
361:   DMSetOptionsPrefix(pack,"pack_");
362:   DMCompositeAddDM(pack,dau);
363:   DMCompositeAddDM(pack,dak);
364:   DMDASetFieldName(dau,0,"u");
365:   DMDASetFieldName(dak,0,"k");
366:   DMSetFromOptions(pack);

368:   DMCreateGlobalVector(pack,&X);
369:   VecDuplicate(X,&F);

371:   PetscNew(struct _UserCtx,&user);

373:   user->pack = pack;

375:   DMCompositeGetGlobalISs(pack,&isg);
376:   DMCompositeGetLocalVectors(pack,&user->Uloc,&user->Kloc);
377:   DMCompositeScatter(pack,X,user->Uloc,user->Kloc);

379:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Coupled problem options","SNES");
380:   {
381:     user->ptype = 0; view_draw = PETSC_FALSE; pass_dm = PETSC_TRUE;

383:     PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,0);
384:     PetscOptionsBool("-view_draw","Draw the final coupled solution regardless of whether only one physics was solved",0,view_draw,&view_draw,0);
385:     PetscOptionsBool("-pass_dm","Pass the packed DM to SNES to use when determining splits and forward into splits",0,pass_dm,&pass_dm,0);
386:   }
387:   PetscOptionsEnd();

389:   FormInitial_Coupled(user,X);

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

447:     PetscOptionsGetInt(0,"-col",&col,0);
448:     PetscOptionsGetBool(0,"-mult_dup",&mult_dup,0);
449:     PetscOptionsGetBool(0,"-view_dup",&view_dup,0);

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

467:   DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc);
468:   PetscFree(user);

470:   ISDestroy(&isg[0]);
471:   ISDestroy(&isg[1]);
472:   PetscFree(isg);
473:   VecDestroy(&X);
474:   VecDestroy(&F);
475:   MatDestroy(&B);
476:   DMDestroy(&dau);
477:   DMDestroy(&dak);
478:   DMDestroy(&pack);
479:   SNESDestroy(&snes);
480:   PetscFinalize();
481:   return 0;
482: }