Actual source code: ex10.c

petsc-3.4.3 2013-10-15
```  1: static const char help[] = "Uses 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 <petscdmmesh.h>
28: #include <petscdmcomposite.h>
29: #include <petscdmda.h>

33: PetscErrorCode CreateProblem_gen_0(DM dm, const char *name)
34: {
35:   ALE::Obj<PETSC_MESH_TYPE> m;
36:   PetscErrorCode            0;

39:   DMMeshGetMesh(dm, m);
40:   {
41:     const ALE::Obj<ALE::Discretization>& d = new ALE::Discretization(m->comm(), m->debug());

43:     d->setNumDof(0, 1);
44:     d->setNumDof(1, 0);
45:     m->setDiscretization(name, d);
46:   }
47:   return(0);
48: }

50: typedef enum {NEUMANN, DIRICHLET} BCType;

52: typedef struct _UserCtx *User;
53: struct _UserCtx {
54:   PetscInt  ptype;
55:   DM        pack;
56:   Vec       Uloc,Kloc;
57:   PetscReal hxu, hxk;
58: };

62: static PetscErrorCode FormFunctionLocal_U(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, SectionReal sectionF, User user)
63: {
64:   ALE::Obj<PETSC_MESH_TYPE> meshU;
65:   ALE::Obj<PETSC_MESH_TYPE> meshK;
66:   PetscErrorCode            ierr;

69:   DMMeshGetMesh(dmu, meshU);
70:   DMMeshGetMesh(dmk, meshK);
71:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesU = meshU->depthStratum(0);
72:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);
73:   PETSC_MESH_TYPE::label_sequence::iterator        vur_iter  = verticesU->begin();
74:   PETSC_MESH_TYPE::point_type                      ulp       = -1;
75:   PETSC_MESH_TYPE::point_type                      urp       = *(++vur_iter);
76:   PETSC_MESH_TYPE::point_type                      klp       = -1;
77:   PetscReal                                        hx        = user->hxu;

79:   PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Starting U residual\n");
80:   for (PETSC_MESH_TYPE::label_sequence::iterator vu_iter = verticesU->begin(), vk_iter = verticesK->begin(); vu_iter != verticesU->end(); ++vu_iter,  ++vk_iter) {
81:     PETSC_MESH_TYPE::point_type up     = *vu_iter;
82:     PETSC_MESH_TYPE::point_type kp     = *vk_iter;
83:     const PetscInt              marker = meshU->getValue(meshU->getLabel("marker"), *vu_iter, 0);
84:     PetscScalar                 values[1];
85:     PetscScalar                 *u;

87:     SectionRealRestrict(sectionU, up, &u);
88:     if (marker == 1) { /* Left end */
89:       values[0] = hx*u[0];
90:       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "  Left  End vu %d hx: %g f %g\n", up, hx, values[0]);
91:       urp       = *(++vur_iter);
92:     } else if (marker == 2) { /* Right end */
93:       values[0] = hx*(u[0] - 1.0);
94:       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "  Right End vu %d hx: %g f %g\n", up, hx, values[0]);
95:     } else if (marker == 3) { /* Left Ghost */
96:       urp = *(++vur_iter);
97:     } else if (marker == 4) { /* Right Ghost */
98:     } else {
99:       PetscScalar *ul, *ur, *k, *kl;

101:       SectionRealRestrict(sectionU, urp, &ur);
102:       SectionRealRestrict(sectionU, ulp, &ul);
103:       SectionRealRestrict(sectionK, kp,  &k);
104:       SectionRealRestrict(sectionK, klp, &kl);

106:       values[0] = hx*((kl[0]*(u[0]-ul[0]) - k[0]*(ur[0]-u[0]))/(hx*hx) - 1.0);

108:       PetscSynchronizedPrintf(PETSC_COMM_WORLD, "  vu %d hx: %g ul %g u %g ur %g kl %g k %g f %g\n", up, hx, ul[0], u[0], ur[0], kl[0], k[0], values[0]);
109:       urp  = *(++vur_iter);
110:     }
111:     SectionRealUpdate(sectionF, up, values, INSERT_VALUES);
112:     ulp  = up;
113:     klp  = kp;
114:   }
115:   PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Ending U residual\n");
116:   PetscSynchronizedFlush(PETSC_COMM_WORLD);
117:   return(0);
118: }

122: static PetscErrorCode FormFunctionLocal_K(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, SectionReal sectionF, User user)
123: {
124:   ALE::Obj<PETSC_MESH_TYPE> meshU;
125:   ALE::Obj<PETSC_MESH_TYPE> meshK;
126:   PetscErrorCode            ierr;

129:   DMMeshGetMesh(dmu, meshU);
130:   DMMeshGetMesh(dmk, meshK);
131:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesU = meshU->depthStratum(0);
132:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);
133:   PETSC_MESH_TYPE::label_sequence::iterator        vu_iter   = verticesU->begin();
134:   PETSC_MESH_TYPE::point_type                      up        = *vu_iter;
135:   PETSC_MESH_TYPE::point_type                      urp;
136:   PetscReal                                        hx = user->hxk;

138:   /*PetscPrintf(PETSC_COMM_WORLD, "Starting K residual\n");*/
139:   for (PETSC_MESH_TYPE::label_sequence::iterator vk_iter = verticesK->begin(); vk_iter != verticesK->end(); ++vk_iter) {
140:     PetscScalar values[1];
141:     PetscScalar *u, *ur, *k;

143:     urp  = *(++vu_iter);
144:     SectionRealRestrict(sectionU, up, &u);
145:     SectionRealRestrict(sectionU, urp, &ur);
146:     SectionRealRestrict(sectionK, *vk_iter, &k);
147:     const PetscScalar ubar  = 0.5*(ur[0] + u[0]);
148:     const PetscScalar gradu = (ur[0] - u[0])/hx;
150:     const PetscScalar w     = 1.0/(1.0 + ubar) + 1.0/g;

152:     values[0] = hx*(PetscExpScalar(k[0]-1.0) + k[0] - 1.0/w);
153:     /*PetscPrintf(PETSC_COMM_WORLD, "  vk %d vu %d vur %d: ubar %g gradu %g g %g w %g f %g\n", *vk_iter, up, urp, ubar, gradu, g, w, values[0]); */
154:     SectionRealUpdate(sectionF, *vk_iter, values, INSERT_VALUES);

156:     up = urp;
157:   }
158:   /*PetscPrintf(PETSC_COMM_WORLD, "Ending K residual\n");*/
159:   return(0);
160: }

164: static PetscErrorCode FormFunction_All(SNES snes, Vec X, Vec F, void *ctx)
165: {
166:   User           user = (User) ctx;
167:   DM             dmu, dmk;
168:   Vec            Uloc, Kloc, Fu, Fk, vecFu, vecFk, vecU, vecK;
169:   SectionReal    sectionU, sectionK, sectionFu, sectionFk;

173:   DMCompositeGetEntries(user->pack, &dmu, &dmk);
174:   switch (user->ptype) {
175:   case 0:
176:     DMMeshGetSectionReal(dmu, "default", &sectionU);
177:     SectionRealCreateLocalVector(sectionU, &vecU);
178:     DMGlobalToLocalBegin(dmu, X, INSERT_VALUES, vecU);
179:     DMGlobalToLocalEnd(dmu, X, INSERT_VALUES, vecU);
180:     VecDestroy(&vecU);

182:     DMMeshGetSectionReal(dmk, "default", &sectionK);
183:     SectionRealCreateLocalVector(sectionK, &vecK);
184:     VecCopy(user->Kloc, vecK);
185:     VecDestroy(&vecK);

187:     SectionRealDuplicate(sectionU, &sectionFu);

189:     FormFunctionLocal_U(dmu, dmk, sectionU, sectionK, sectionFu, user);
190:     SectionRealCreateLocalVector(sectionFu, &vecFu);
191:     DMLocalToGlobalBegin(dmu, vecFu, INSERT_VALUES, F);
192:     DMLocalToGlobalEnd(dmu, vecFu, INSERT_VALUES, F);
193:     VecDestroy(&vecFu);
194:     SectionRealDestroy(&sectionFu);
195:     break;
196:   case 1:
197:     DMMeshGetSectionReal(dmk, "default", &sectionK);
198:     SectionRealCreateLocalVector(sectionK, &vecK);
199:     DMGlobalToLocalBegin(dmk, X, INSERT_VALUES, vecK);
200:     DMGlobalToLocalEnd(dmk, X, INSERT_VALUES, vecK);
201:     VecDestroy(&vecK);

203:     DMMeshGetSectionReal(dmu, "default", &sectionU);
204:     SectionRealCreateLocalVector(sectionU, &vecU);
205:     VecCopy(user->Uloc, vecU);
206:     VecDestroy(&vecU);

208:     SectionRealDuplicate(sectionK, &sectionFk);

210:     FormFunctionLocal_K(dmu, dmk, sectionU, sectionK, sectionFk, user);
211:     SectionRealCreateLocalVector(sectionFk, &vecFk);
212:     DMLocalToGlobalBegin(dmk, vecFk, INSERT_VALUES, F);
213:     DMLocalToGlobalEnd(dmk, vecFk, INSERT_VALUES, F);
214:     VecDestroy(&vecFk);
215:     SectionRealDestroy(&sectionFk);
216:     break;
217:   case 2:
218:     DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc);
219:     DMCompositeScatter(user->pack, X, Uloc, Kloc);

221:     DMMeshGetSectionReal(dmu, "default", &sectionU);
222:     SectionRealCreateLocalVector(sectionU, &vecU);
223:     VecCopy(Uloc, vecU);
224:     VecDestroy(&vecU);

226:     DMMeshGetSectionReal(dmk, "default", &sectionK);
227:     SectionRealCreateLocalVector(sectionK, &vecK);
228:     VecCopy(Kloc, vecK);
229:     VecDestroy(&vecK);

231:     DMCompositeRestoreLocalVectors(user->pack, &Uloc, &Kloc);
232:     SectionRealDuplicate(sectionU, &sectionFu);
233:     SectionRealDuplicate(sectionK, &sectionFk);

235:     FormFunctionLocal_U(dmu, dmk, sectionU, sectionK, sectionFu, user);
236:     FormFunctionLocal_K(dmu, dmk, sectionU, sectionK, sectionFk, user);
237:     DMCompositeGetLocalVectors(user->pack, &Fu, &Fk);
238:     SectionRealCreateLocalVector(sectionFu, &vecFu);
239:     VecCopy(vecFu, Fu);
240:     VecDestroy(&vecFu);
241:     SectionRealCreateLocalVector(sectionFk, &vecFk);
242:     VecCopy(vecFk, Fk);
243:     VecDestroy(&vecFk);
244:     DMCompositeGather(user->pack, F, INSERT_VALUES, Fu, Fk);
245:     SectionRealDestroy(&sectionFu);
246:     SectionRealDestroy(&sectionFk);
247:     break;
248:   }
249:   SectionRealDestroy(&sectionU);
250:   SectionRealDestroy(&sectionK);
251:   return(0);
252: }

256: static PetscErrorCode FormJacobianLocal_U(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, Mat Buu, User user)
257: {
258:   ALE::Obj<PETSC_MESH_TYPE> meshU;
259:   ALE::Obj<PETSC_MESH_TYPE> meshK;
260:   PetscErrorCode            ierr;

263:   DMMeshGetMesh(dmu, meshU);
264:   DMMeshGetMesh(dmk, meshK);
265:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesU = meshU->depthStratum(0);
266:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);
267:   PETSC_MESH_TYPE::label_sequence::iterator        vur_iter  = verticesU->begin();
268:   PETSC_MESH_TYPE::point_type                      klp       = -1;
269:   PETSC_MESH_TYPE::point_type                      ulp       = -1;
270:   PETSC_MESH_TYPE::point_type                      urp       = *(++vur_iter);
271:   PetscReal                                        hx        = user->hxu;

273:   for (PETSC_MESH_TYPE::label_sequence::iterator vu_iter = verticesU->begin(), vk_iter = verticesK->begin(); vu_iter != verticesU->end(); ++vu_iter,  ++vk_iter) {
274:     PETSC_MESH_TYPE::point_type up     = *vu_iter;
275:     PETSC_MESH_TYPE::point_type kp     = *vk_iter;
276:     const PetscInt              marker = meshU->getValue(meshU->getLabel("marker"), *vu_iter, 0);
277:     PetscScalar                 values[3];

279:     if (marker == 1) {
280:       values[0] = hx;
281:       std::cout << "["<<meshU->commRank()<<"]: row " << up << " Left BC" << std::endl;
282:       MatSetValuesTopology(Buu, dmu, 1, &up, dmu, 1, &up, values, INSERT_VALUES);
283:       urp  = *(++vur_iter);
284:     } else if (marker == 2) {
285:       values[0] = hx;
286:       std::cout << "["<<meshU->commRank()<<"]: row " << up << " Right BC" << std::endl;
287:       MatSetValuesTopology(Buu, dmu, 1, &up, dmu, 1, &up, values, INSERT_VALUES);
288:     } else {
289:       PetscScalar *k, *kl;
290:       PetscInt    cols[3] = {ulp, up, urp};

292:       SectionRealRestrict(sectionK, kp,  &k);
293:       SectionRealRestrict(sectionK, klp, &kl);

295:       values[0] = -kl[0]/hx;
296:       values[1] = (kl[0]+k[0])/hx;
297:       values[2] = -k[0]/hx;
298:       std::cout << "["<<meshU->commRank()<<"]: row " << up << " cols " << cols[0] <<", "<< cols[1] <<", "<< cols[2] << std::endl;
299:       MatSetValuesTopology(Buu, dmu, 1, &up, dmu, 3, cols, values, INSERT_VALUES);
300:       urp  = *(++vur_iter);
301:     }
302:     ulp = up;
303:     klp = kp;
304:   }
305:   return(0);
306: }

310: static PetscErrorCode FormJacobianLocal_K(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, Mat Bkk, User user)
311: {
312:   ALE::Obj<PETSC_MESH_TYPE> meshU;
313:   ALE::Obj<PETSC_MESH_TYPE> meshK;
314:   PetscErrorCode            ierr;

317:   DMMeshGetMesh(dmu, meshU);
318:   DMMeshGetMesh(dmk, meshK);
319:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);

321:   PetscReal hx = user->hxk;

323:   for (PETSC_MESH_TYPE::label_sequence::iterator vk_iter = verticesK->begin(); vk_iter != verticesK->end(); ++vk_iter) {
324:     PETSC_MESH_TYPE::point_type kp = *vk_iter;
325:     PetscScalar                 values[1];
326:     PetscScalar                 *k;

328:     SectionRealRestrict(sectionK, kp,  &k);
329:     values[0] = hx*(PetscExpScalar(k[0] - 1.0) + 1.0);
330:     MatSetValuesTopology(Bkk, dmk, 1, &kp, dmk, 1, &kp, values, INSERT_VALUES);
331:   }
332:   return(0);
333: }

337: static PetscErrorCode FormJacobianLocal_UK(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, Mat Buk, User user)
338: {
339:   ALE::Obj<PETSC_MESH_TYPE> meshU;
340:   ALE::Obj<PETSC_MESH_TYPE> meshK;
341:   PetscErrorCode            ierr;

344:   if (!Buk) return(0); /* Not assembling this block */
345:   DMMeshGetMesh(dmu, meshU);
346:   DMMeshGetMesh(dmk, meshK);
347:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesU = meshU->depthStratum(0);
348:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);
349:   PETSC_MESH_TYPE::label_sequence::iterator        vur_iter  = verticesU->begin();
350:   PETSC_MESH_TYPE::point_type                      ulp       = -1;
351:   PETSC_MESH_TYPE::point_type                      urp       = *(++vur_iter);
352:   PETSC_MESH_TYPE::point_type                      klp       = -1;
353:   PetscReal                                        hx        = user->hxu;

355:   for (PETSC_MESH_TYPE::label_sequence::iterator vu_iter = verticesU->begin(), vk_iter = verticesK->begin(); vu_iter != verticesU->end(); ++vu_iter, ++vk_iter) {
356:     PETSC_MESH_TYPE::point_type up     = *vu_iter;
357:     PETSC_MESH_TYPE::point_type kp     = *vk_iter;
358:     const PetscInt              marker = meshU->getValue(meshU->getLabel("marker"), *vu_iter, 0);
359:     PetscScalar                 values[3];

361:     if (marker == 1) {
362:       ulp = up;
363:       urp = *(++vur_iter);
364:       klp = kp;
365:       continue;
366:     }
367:     if (marker == 2) continue;
368:     PetscInt    cols[3] = {klp, kp};
369:     PetscScalar *u, *ul, *ur;

371:     SectionRealRestrict(sectionU, up,  &u);
372:     SectionRealRestrict(sectionU, ulp, &ul);
373:     SectionRealRestrict(sectionU, urp, &ur);

375:     values[0] = (u[0]-ul[0])/hx;
376:     values[1] = (u[0]-ur[0])/hx;

378:     MatSetValuesTopology(Buk, dmu, 1, &up, dmk, 2, cols, values, INSERT_VALUES);
379:     ulp  = up;
380:     urp  = *(++vur_iter);
381:     klp  = kp;
382:   }
383:   return(0);
384: }

388: static PetscErrorCode FormJacobianLocal_KU(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, Mat Bku, User user)
389: {
390:   ALE::Obj<PETSC_MESH_TYPE> meshU;
391:   ALE::Obj<PETSC_MESH_TYPE> meshK;
392:   PetscErrorCode            ierr;

395:   if (!Bku) return(0); /* Not assembling this block */
396:   DMMeshGetMesh(dmu, meshU);
397:   DMMeshGetMesh(dmk, meshK);
398:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesU = meshU->depthStratum(0);
399:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);
400:   PETSC_MESH_TYPE::label_sequence::iterator        vur_iter  = verticesU->begin();
401:   PETSC_MESH_TYPE::label_sequence::iterator        vkr_iter  = verticesK->begin();
402:   PETSC_MESH_TYPE::point_type                      urp       = *(++vur_iter);
403:   PetscReal                                        hx        = user->hxk;

405:   for (PETSC_MESH_TYPE::label_sequence::iterator vk_iter = verticesK->begin(), vu_iter = verticesU->begin(); vk_iter != verticesK->end(); ++vk_iter,  ++vu_iter) {
406:     PETSC_MESH_TYPE::point_type up      = *vu_iter;
407:     PETSC_MESH_TYPE::point_type kp      = *vk_iter;
408:     PetscInt                    cols[2] = {up, urp};
409:     PetscScalar                 values[2];
410:     PetscScalar                 *u, *ur;

412:     SectionRealRestrict(sectionU, up,  &u);
413:     SectionRealRestrict(sectionU, urp, &ur);
414:     const PetscScalar
415:       ubar     = 0.5*(u[0]+ur[0]),
416:       ubar_L   = 0.5,
417:       ubar_R   = 0.5,
421:       g        = 1.0 + PetscSqr(gradu),
423:       w        = 1.0/(1.0+ubar) + 1.0/g,
424:       w_ubar   = -1./PetscSqr(1.+ubar),
426:       iw       = 1.0/w,
427:       iw_ubar  = -w_ubar * PetscSqr(iw),

432:     MatSetValuesTopology(Bku, dmk, 1, &kp, dmu, 2, cols, values, INSERT_VALUES);
433:     urp       = *(++vur_iter);
434:   }
435:   return(0);
436: }

440: static PetscErrorCode FormJacobian_All(SNES snes, Vec X, Mat *J, Mat *B, MatStructure *mstr, void *ctx)
441: {
442:   User           user = (User) ctx;
443:   DM             dmu, dmk;
444:   Vec            Uloc, Kloc, vecU, vecK;
445:   SectionReal    sectionU, sectionK;

449:   DMCompositeGetEntries(user->pack, &dmu, &dmk);
450:   switch (user->ptype) {
451:   case 0:
452:     DMMeshGetSectionReal(dmu, "default", &sectionU);
453:     SectionRealCreateLocalVector(sectionU, &vecU);
454:     DMGlobalToLocalBegin(dmu, X, INSERT_VALUES, vecU);
455:     DMGlobalToLocalEnd(dmu, X, INSERT_VALUES, vecU);
456:     VecDestroy(&vecU);

458:     DMMeshGetSectionReal(dmk, "default", &sectionK);
459:     SectionRealCreateLocalVector(sectionK, &vecK);
460:     VecCopy(user->Kloc, vecK);
461:     VecDestroy(&vecK);

463:     FormJacobianLocal_U(dmu, dmk, sectionU, sectionK, *B, user);
464:     SectionRealDestroy(&sectionU);
465:     SectionRealDestroy(&sectionK);
466:     break;
467:   case 1:
468:     DMMeshGetSectionReal(dmk, "default", &sectionK);
469:     SectionRealCreateLocalVector(sectionK, &vecK);
470:     DMGlobalToLocalBegin(dmk, X, INSERT_VALUES, vecK);
471:     DMGlobalToLocalEnd(dmk, X, INSERT_VALUES, vecK);
472:     VecDestroy(&vecK);

474:     DMMeshGetSectionReal(dmu, "default", &sectionU);
475:     SectionRealCreateLocalVector(sectionU, &vecU);
476:     VecCopy(user->Uloc, vecU);
477:     VecDestroy(&vecU);

479:     FormJacobianLocal_K(dmu, dmk, sectionU, sectionK, *B, user);
480:     SectionRealDestroy(&sectionU);
481:     SectionRealDestroy(&sectionK);
482:     break;
483:   case 2: {
484:     Mat Buu,Buk,Bku,Bkk;
485:     IS  *is;

487:     DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc);
488:     DMCompositeScatter(user->pack, X, Uloc, Kloc);

490:     DMMeshGetSectionReal(dmu, "default", &sectionU);
491:     SectionRealCreateLocalVector(sectionU, &vecU);
492:     VecCopy(Uloc, vecU);
493:     VecDestroy(&vecU);

495:     DMMeshGetSectionReal(dmk, "default", &sectionK);
496:     SectionRealCreateLocalVector(sectionK, &vecK);
497:     VecCopy(Kloc, vecK);
498:     VecDestroy(&vecK);

500:     DMCompositeRestoreLocalVectors(user->pack, &Uloc, &Kloc);

502:     DMCompositeGetLocalISs(user->pack, &is);
503:     MatGetLocalSubMatrix(*B, is[0], is[0], &Buu);
504:     MatGetLocalSubMatrix(*B, is[0], is[1], &Buk);
505:     MatGetLocalSubMatrix(*B, is[1], is[0], &Bku);
506:     MatGetLocalSubMatrix(*B, is[1], is[1], &Bkk);
507:     FormJacobianLocal_U (dmu, dmk, sectionU, sectionK, Buu, user);
508:     FormJacobianLocal_UK(dmu, dmk, sectionU, sectionK, Buk, user);
509:     FormJacobianLocal_KU(dmu, dmk, sectionU, sectionK, Bku, user);
510:     FormJacobianLocal_K (dmu, dmk, sectionU, sectionK, Bkk, user);
511:     MatRestoreLocalSubMatrix(*B, is[0], is[0], &Buu);
512:     MatRestoreLocalSubMatrix(*B, is[0], is[1], &Buk);
513:     MatRestoreLocalSubMatrix(*B, is[1], is[0], &Bku);
514:     MatRestoreLocalSubMatrix(*B, is[1], is[1], &Bkk);

516:     SectionRealDestroy(&sectionU);
517:     SectionRealDestroy(&sectionK);
518:     ISDestroy(&is[0]);
519:     ISDestroy(&is[1]);
520:     PetscFree(is);
521:   } break;
522:   }
523:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
524:   MatAssemblyEnd  (*B, MAT_FINAL_ASSEMBLY);
525:   if (*J != *B) {
526:     MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);
527:     MatAssemblyEnd  (*J, MAT_FINAL_ASSEMBLY);
528:   }
529:   /*MatView(*B, PETSC_VIEWER_STDOUT_WORLD);*/
530:   *mstr = DIFFERENT_NONZERO_PATTERN;
531:   return(0);
532: }

536: static PetscErrorCode FormInitial_Coupled(User user,Vec X)
537: {
538:   ALE::Obj<PETSC_MESH_TYPE> meshU;
539:   ALE::Obj<PETSC_MESH_TYPE> meshK;
540:   DM                        dmu,          dmk;
541:   SectionReal               coordinatesU, coordinatesK;
542:   SectionReal               sectionU,     sectionK;
543:   Vec                       vecU,         vecK;
544:   PetscErrorCode            ierr;

547:   DMCompositeGetEntries(user->pack, &dmu, &dmk);
548:   DMMeshGetMesh(dmu, meshU);
549:   DMMeshGetMesh(dmk, meshK);
550:   DMMeshGetSectionReal(dmu, "coordinates", &coordinatesU);
551:   DMMeshGetSectionReal(dmk, "coordinates", &coordinatesK);
552:   DMMeshGetSectionReal(dmu, "default", &sectionU);
553:   DMMeshGetSectionReal(dmk, "default", &sectionK);
554:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesU = meshU->depthStratum(0);
555:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);

557:   for (PETSC_MESH_TYPE::label_sequence::iterator v_iter = verticesU->begin(); v_iter != verticesU->end(); ++v_iter) {
558:     PetscScalar values[1];
559:     PetscScalar *coords;

561:     SectionRealRestrict(coordinatesU, *v_iter, &coords);
562:     values[0] = coords[0]*(1.0 - coords[0]);
563:     SectionRealUpdate(sectionU, *v_iter, values, INSERT_VALUES);
564:   }
565:   for (PETSC_MESH_TYPE::label_sequence::iterator v_iter = verticesK->begin(); v_iter != verticesK->end(); ++v_iter) {
566:     PetscScalar values[1];
567:     PetscScalar *coords;

569:     SectionRealRestrict(coordinatesK, *v_iter, &coords);
570:     values[0] = (PetscScalar) (1.0 + 0.5*sin((double) 2*PETSC_PI*coords[0]));
571:     SectionRealUpdate(sectionK, *v_iter, values, INSERT_VALUES);
572:   }
573:   SectionRealCreateLocalVector(sectionU, &vecU);
574:   SectionRealCreateLocalVector(sectionK, &vecK);
575:   VecCopy(vecU, user->Uloc);
576:   VecCopy(vecK, user->Kloc);
577:   VecDestroy(&vecU);
578:   VecDestroy(&vecK);
579:   SectionRealDestroy(&coordinatesU);
580:   SectionRealDestroy(&coordinatesK);
581:   SectionRealDestroy(&sectionU);
582:   SectionRealDestroy(&sectionK);
583:   DMCompositeGather(user->pack, X, INSERT_VALUES, user->Uloc, user->Kloc);
584:   return(0);
585: }

589: int main(int argc, char *argv[])
590: {
591:   DM             dau, dak, dmu, dmk, pack;
592:   User           user;
593:   SNES           snes;
594:   KSP            ksp;
595:   PC             pc;
596:   Mat            B;
597:   Vec            X,F,Xu,Xk,Fu,Fk;
598:   IS             *isg;
599:   const PetscInt *lxu;
600:   PetscInt       *lxk, m, nprocs;
601:   PetscBool      view_draw;

604:   PetscInitialize(&argc,&argv,0,help);
605:   /* Create meshes */
606:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-10,1,1,NULL,&dau);
607:   DMDAGetOwnershipRanges(dau,&lxu,0,0);
608:   DMDAGetInfo(dau,0, &m,0,0, &nprocs,0,0, 0,0,0,0,0,0);
609:   PetscMalloc(nprocs*sizeof(*lxk),&lxk);
610:   PetscMemcpy(lxk,lxu,nprocs*sizeof(*lxk));
611:   lxk[0]--;
612:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,m-1,1,1,lxk,&dak);
613:   PetscFree(lxk);
614:   DMDASetUniformCoordinates(dau, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
615:   DMDASetUniformCoordinates(dak, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
616:   DMConvert(dau, DMMESH, &dmu);
617:   DMConvert(dak, DMMESH, &dmk);
618:   DMSetOptionsPrefix(dmu, "u_");
619:   DMSetOptionsPrefix(dmk, "k_");
620:   DMDestroy(&dau);
621:   DMDestroy(&dak);

623: #if 1
624:   {
625:     ALE::Obj<PETSC_MESH_TYPE> m;
626:     DMMeshGetMesh(dmu, m);
627:     m->view("Mesh");
628:     m->getLabel("marker")->view("Marker");
629:   }
630: #endif

632:   PetscNew(struct _UserCtx, &user);
633:   user->hxu = 1.0/m;
634:   user->hxk = 1.0/(m-1);
635:   /* Setup dof layout.
636:    For a DMDA, this is automatic given the number of dof at each vertex. However, for a DMMesh, we need to specify this.
637:   */
638:   {
639:     /* There is perhaps a better way to do this that does not rely on the Discretization/BoundaryCondition objects in Mesh.hh */
640:     CreateProblem_gen_0(dmu, "u");
641:     CreateProblem_gen_0(dmk, "k");
642:   }
643:   SectionReal defaultSection;

645:   DMMeshGetSectionReal(dmu, "default", &defaultSection);
646:   DMMeshSetupSection(dmu, defaultSection);
647:   SectionRealDestroy(&defaultSection);
648:   DMMeshGetSectionReal(dmk, "default", &defaultSection);
649:   DMMeshSetupSection(dmk, defaultSection);
650:   SectionRealDestroy(&defaultSection);

652:   DMCompositeCreate(PETSC_COMM_WORLD, &pack);
655:   DMSetOptionsPrefix(pack, "pack_");

657:   DMCreateGlobalVector(pack, &X);
658:   VecDuplicate(X, &F);

660:   user->pack = pack;

662:   DMCompositeGetGlobalISs(pack, &isg);
663:   DMCompositeGetLocalVectors(pack, &user->Uloc, &user->Kloc);
664:   DMCompositeScatter(pack, X, user->Uloc, user->Kloc);

666:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Coupled problem options","SNES");
667:   {
668:     user->ptype = 0; view_draw = PETSC_FALSE;

670:     PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,0);
671:     PetscOptionsBool("-view_draw","Draw the final coupled solution regardless of whether only one physics was solved",0,view_draw,&view_draw,0);
672:   }
673:   PetscOptionsEnd();

675:   FormInitial_Coupled(user,X);
676:   VecView(X, PETSC_VIEWER_STDOUT_WORLD);

678:   SNESCreate(PETSC_COMM_WORLD,&snes);
679:   switch (user->ptype) {
680:   case 0:
681:     DMCompositeGetAccess(pack,X,&Xu,0);
682:     DMCompositeGetAccess(pack,F,&Fu,0);
683:     DMCreateMatrix(dmu,NULL,&B);
684:     SNESSetFunction(snes,Fu,FormFunction_All,user);
685:     SNESSetJacobian(snes,B,B,FormJacobian_All,user);
686:     SNESSetFromOptions(snes);
687:     SNESSolve(snes,NULL,Xu);
688:     DMCompositeRestoreAccess(pack,X,&Xu,0);
689:     DMCompositeRestoreAccess(pack,F,&Fu,0);
690:     break;
691:   case 1:
692:     DMCompositeGetAccess(pack,X,0,&Xk);
693:     DMCompositeGetAccess(pack,F,0,&Fk);
694:     DMCreateMatrix(dmk,NULL,&B);
695:     SNESSetFunction(snes,Fk,FormFunction_All,user);
696:     SNESSetJacobian(snes,B,B,FormJacobian_All,user);
697:     SNESSetFromOptions(snes);
698:     SNESSolve(snes,NULL,Xk);
699:     DMCompositeRestoreAccess(pack,X,0,&Xk);
700:     DMCompositeRestoreAccess(pack,F,0,&Fk);
701:     break;
702:   case 2:
703:     DMCreateMatrix(pack,NULL,&B);
704:     SNESSetFunction(snes,F,FormFunction_All,user);
705:     SNESSetJacobian(snes,B,B,FormJacobian_All,user);
706:     SNESSetFromOptions(snes);
707:     SNESGetKSP(snes,&ksp);
708:     KSPGetPC(ksp,&pc);
709:     PCFieldSplitSetIS(pc,"u",isg[0]);
710:     PCFieldSplitSetIS(pc,"k",isg[1]);
711:     SNESSolve(snes,NULL,X);
712:     break;
713:   }
714:   if (view_draw) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
715:   if (0) {
716:     PetscInt  col      = 0;
717:     PetscBool mult_dup = PETSC_FALSE,view_dup = PETSC_FALSE;
718:     Mat       D;
719:     Vec       Y;

721:     PetscOptionsGetInt(0,"-col",&col,0);
722:     PetscOptionsGetBool(0,"-mult_dup",&mult_dup,0);
723:     PetscOptionsGetBool(0,"-view_dup",&view_dup,0);

725:     VecDuplicate(X,&Y);
726:     /* MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); */
727:     /* MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); */
728:     MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&D);
729:     VecZeroEntries(X);
730:     VecSetValue(X,col,1.0,INSERT_VALUES);
731:     VecAssemblyBegin(X);
732:     VecAssemblyEnd(X);
733:     MatMult(mult_dup ? D : B,X,Y);
734:     MatView(view_dup ? D : B,PETSC_VIEWER_STDOUT_WORLD);
735:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
736:     VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
737:     MatDestroy(&D);
738:     VecDestroy(&Y);
739:   }

741:   DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc);
742:   PetscFree(user);

744:   ISDestroy(&isg[0]);
745:   ISDestroy(&isg[1]);
746:   PetscFree(isg);
747:   VecDestroy(&X);
748:   VecDestroy(&F);
749:   MatDestroy(&B);
750:   DMDestroy(&dmu);
751:   DMDestroy(&dmk);
752:   DMDestroy(&pack);
753:   SNESDestroy(&snes);
754:   PetscFree(user);
755:   PetscFinalize();
756:   return 0;
757: }
```