Actual source code: ex10.c

petsc-3.3-p7 2013-05-11
```  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;

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);
105:       values[0] = hx*((kl[0]*(u[0]-ul[0]) - k[0]*(ur[0]-u[0]))/(hx*hx) - 1.0);
106:       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]);
107:       urp  = *(++vur_iter);
108:     }
109:     SectionRealUpdate(sectionF, up, values, INSERT_VALUES);
110:     ulp  = up;
111:     klp  = kp;
112:   }
113:   PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Ending U residual\n");
114:   PetscSynchronizedFlush(PETSC_COMM_WORLD);
115:   return(0);
116: }

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

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

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

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

150:     values[0] = hx*(PetscExpScalar(k[0]-1.0) + k[0] - 1.0/w);
151:     //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]);
152:     SectionRealUpdate(sectionF, *vk_iter, values, INSERT_VALUES);

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

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

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

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

185:     SectionRealDuplicate(sectionU, &sectionFu);

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

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

206:     SectionRealDuplicate(sectionK, &sectionFk);

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

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

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

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

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

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

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

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

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

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

307: static PetscErrorCode FormJacobianLocal_K(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, Mat Bkk, User user)
308: {
309:   ALE::Obj<PETSC_MESH_TYPE> meshU;
310:   ALE::Obj<PETSC_MESH_TYPE> meshK;

314:   DMMeshGetMesh(dmu, meshU);
315:   DMMeshGetMesh(dmk, meshK);
316:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);
317:   PetscReal hx = user->hxk;

319:   for(PETSC_MESH_TYPE::label_sequence::iterator vk_iter = verticesK->begin(); vk_iter != verticesK->end(); ++vk_iter) {
320:     PETSC_MESH_TYPE::point_type kp = *vk_iter;
321:     PetscScalar                 values[1];
322:     PetscScalar                *k;

324:     SectionRealRestrict(sectionK, kp,  &k);
325:     values[0] = hx*(PetscExpScalar(k[0] - 1.0) + 1.0);
326:     MatSetValuesTopology(Bkk, dmk, 1, &kp, dmk, 1, &kp, values, INSERT_VALUES);
327:   }
328:   return(0);
329: }

333: static PetscErrorCode FormJacobianLocal_UK(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, Mat Buk, User user)
334: {
335:   ALE::Obj<PETSC_MESH_TYPE> meshU;
336:   ALE::Obj<PETSC_MESH_TYPE> meshK;

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

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

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

367:     SectionRealRestrict(sectionU, up,  &u);
368:     SectionRealRestrict(sectionU, ulp, &ul);
369:     SectionRealRestrict(sectionU, urp, &ur);
370:     values[0] = (u[0]-ul[0])/hx;
371:     values[1] = (u[0]-ur[0])/hx;
372:     MatSetValuesTopology(Buk, dmu, 1, &up, dmk, 2, cols, values, INSERT_VALUES);
373:     ulp  = up;
374:     urp  = *(++vur_iter);
375:     klp  = kp;
376:   }
377:   return(0);
378: }

382: static PetscErrorCode FormJacobianLocal_KU(DM dmu, DM dmk, SectionReal sectionU, SectionReal sectionK, Mat Bku, User user)
383: {
384:   ALE::Obj<PETSC_MESH_TYPE> meshU;
385:   ALE::Obj<PETSC_MESH_TYPE> meshK;

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

399:   for(PETSC_MESH_TYPE::label_sequence::iterator vk_iter = verticesK->begin(), vu_iter = verticesU->begin(); vk_iter != verticesK->end(); ++vk_iter,  ++vu_iter) {
400:     PETSC_MESH_TYPE::point_type up      = *vu_iter;
401:     PETSC_MESH_TYPE::point_type kp      = *vk_iter;
402:     PetscInt                    cols[2] = {up, urp};
403:     PetscScalar                 values[2];
404:     PetscScalar                *u, *ur;

406:     SectionRealRestrict(sectionU, up,  &u);
407:     SectionRealRestrict(sectionU, urp, &ur);
408:     const PetscScalar
409:       ubar     = 0.5*(u[0]+ur[0]),
410:       ubar_L   = 0.5,
411:       ubar_R   = 0.5,
415:       g        = 1.0 + PetscSqr(gradu),
417:       w        = 1.0/(1.0+ubar) + 1.0/g,
418:       w_ubar   = -1./PetscSqr(1.+ubar),
420:       iw       = 1.0/w,
421:       iw_ubar  = -w_ubar * PetscSqr(iw),

426:     MatSetValuesTopology(Bku, dmk, 1, &kp, dmu, 2, cols, values, INSERT_VALUES);
427:     urp  = *(++vur_iter);
428:   }
429:   return(0);
430: }

434: static PetscErrorCode FormJacobian_All(SNES snes, Vec X, Mat *J, Mat *B, MatStructure *mstr, void *ctx)
435: {
436:   User           user = (User) ctx;
437:   DM             dmu, dmk;
438:   Vec            Uloc, Kloc, vecU, vecK;
439:   SectionReal    sectionU, sectionK;

443:   DMCompositeGetEntries(user->pack, &dmu, &dmk);
444:   switch (user->ptype) {
445:   case 0:
446:     DMMeshGetSectionReal(dmu, "default", &sectionU);
447:     SectionRealCreateLocalVector(sectionU, &vecU);
448:     DMGlobalToLocalBegin(dmu, X, INSERT_VALUES, vecU);
449:     DMGlobalToLocalEnd(dmu, X, INSERT_VALUES, vecU);
450:     VecDestroy(&vecU);

452:     DMMeshGetSectionReal(dmk, "default", &sectionK);
453:     SectionRealCreateLocalVector(sectionK, &vecK);
454:     VecCopy(user->Kloc, vecK);
455:     VecDestroy(&vecK);

457:     FormJacobianLocal_U(dmu, dmk, sectionU, sectionK, *B, user);
458:     SectionRealDestroy(&sectionU);
459:     SectionRealDestroy(&sectionK);
460:     break;
461:   case 1:
462:     DMMeshGetSectionReal(dmk, "default", &sectionK);
463:     SectionRealCreateLocalVector(sectionK, &vecK);
464:     DMGlobalToLocalBegin(dmk, X, INSERT_VALUES, vecK);
465:     DMGlobalToLocalEnd(dmk, X, INSERT_VALUES, vecK);
466:     VecDestroy(&vecK);

468:     DMMeshGetSectionReal(dmu, "default", &sectionU);
469:     SectionRealCreateLocalVector(sectionU, &vecU);
470:     VecCopy(user->Uloc, vecU);
471:     VecDestroy(&vecU);

473:     FormJacobianLocal_K(dmu, dmk, sectionU, sectionK, *B, user);
474:     SectionRealDestroy(&sectionU);
475:     SectionRealDestroy(&sectionK);
476:     break;
477:   case 2: {
478:     Mat Buu,Buk,Bku,Bkk;
479:     IS *is;

481:     DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc);
482:     DMCompositeScatter(user->pack, X, Uloc, Kloc);

484:     DMMeshGetSectionReal(dmu, "default", &sectionU);
485:     SectionRealCreateLocalVector(sectionU, &vecU);
486:     VecCopy(Uloc, vecU);
487:     VecDestroy(&vecU);

489:     DMMeshGetSectionReal(dmk, "default", &sectionK);
490:     SectionRealCreateLocalVector(sectionK, &vecK);
491:     VecCopy(Kloc, vecK);
492:     VecDestroy(&vecK);

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

496:     DMCompositeGetLocalISs(user->pack, &is);
497:     MatGetLocalSubMatrix(*B, is[0], is[0], &Buu);
498:     MatGetLocalSubMatrix(*B, is[0], is[1], &Buk);
499:     MatGetLocalSubMatrix(*B, is[1], is[0], &Bku);
500:     MatGetLocalSubMatrix(*B, is[1], is[1], &Bkk);
501:     FormJacobianLocal_U (dmu, dmk, sectionU, sectionK, Buu, user);
502:     FormJacobianLocal_UK(dmu, dmk, sectionU, sectionK, Buk, user);
503:     FormJacobianLocal_KU(dmu, dmk, sectionU, sectionK, Bku, user);
504:     FormJacobianLocal_K (dmu, dmk, sectionU, sectionK, Bkk, user);
505:     MatRestoreLocalSubMatrix(*B, is[0], is[0], &Buu);
506:     MatRestoreLocalSubMatrix(*B, is[0], is[1], &Buk);
507:     MatRestoreLocalSubMatrix(*B, is[1], is[0], &Bku);
508:     MatRestoreLocalSubMatrix(*B, is[1], is[1], &Bkk);

510:     SectionRealDestroy(&sectionU);
511:     SectionRealDestroy(&sectionK);
512:     ISDestroy(&is[0]);
513:     ISDestroy(&is[1]);
514:     PetscFree(is);
515:   } break;
516:   }
517:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
518:   MatAssemblyEnd  (*B, MAT_FINAL_ASSEMBLY);
519:   if (*J != *B) {
520:     MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);
521:     MatAssemblyEnd  (*J, MAT_FINAL_ASSEMBLY);
522:   }
523:   //MatView(*B, PETSC_VIEWER_STDOUT_WORLD);
524:   *mstr = DIFFERENT_NONZERO_PATTERN;
525:   return(0);
526: }

530: static PetscErrorCode FormInitial_Coupled(User user,Vec X)
531: {
532:   ALE::Obj<PETSC_MESH_TYPE> meshU;
533:   ALE::Obj<PETSC_MESH_TYPE> meshK;
534:   DM             dmu,          dmk;
535:   SectionReal    coordinatesU, coordinatesK;
536:   SectionReal    sectionU,     sectionK;
537:   Vec            vecU,         vecK;

541:   DMCompositeGetEntries(user->pack, &dmu, &dmk);
542:   DMMeshGetMesh(dmu, meshU);
543:   DMMeshGetMesh(dmk, meshK);
544:   DMMeshGetSectionReal(dmu, "coordinates", &coordinatesU);
545:   DMMeshGetSectionReal(dmk, "coordinates", &coordinatesK);
546:   DMMeshGetSectionReal(dmu, "default", &sectionU);
547:   DMMeshGetSectionReal(dmk, "default", &sectionK);
548:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesU = meshU->depthStratum(0);
549:   const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& verticesK = meshK->depthStratum(0);

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

555:     SectionRealRestrict(coordinatesU, *v_iter, &coords);
556:     values[0] = coords[0]*(1.0 - coords[0]);
557:     SectionRealUpdate(sectionU, *v_iter, values, INSERT_VALUES);
558:   }
559:   for(PETSC_MESH_TYPE::label_sequence::iterator v_iter = verticesK->begin(); v_iter != verticesK->end(); ++v_iter) {
560:     PetscScalar  values[1];
561:     PetscScalar *coords;

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

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

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

617: #if 1
618:   {
619:     ALE::Obj<PETSC_MESH_TYPE> m;
620:     DMMeshGetMesh(dmu, m);
621:     m->view("Mesh");
622:     m->getLabel("marker")->view("Marker");
623:   }
624: #endif

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

639:   DMMeshGetSectionReal(dmu, "default", &defaultSection);
640:   DMMeshSetupSection(dmu, defaultSection);
641:   SectionRealDestroy(&defaultSection);
642:   DMMeshGetSectionReal(dmk, "default", &defaultSection);
643:   DMMeshSetupSection(dmk, defaultSection);
644:   SectionRealDestroy(&defaultSection);

646:   DMCompositeCreate(PETSC_COMM_WORLD, &pack);
649:   DMSetOptionsPrefix(pack, "pack_");

651:   DMCreateGlobalVector(pack, &X);
652:   VecDuplicate(X, &F);

654:   user->pack = pack;
655:   DMCompositeGetGlobalISs(pack, &isg);
656:   DMCompositeGetLocalVectors(pack, &user->Uloc, &user->Kloc);
657:   DMCompositeScatter(pack, X, user->Uloc, user->Kloc);

659:   PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Coupled problem options","SNES");
660:   {
661:     user->ptype = 0; view_draw = PETSC_FALSE;
662:     PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,0);
663:     PetscOptionsBool("-view_draw","Draw the final coupled solution regardless of whether only one physics was solved",0,view_draw,&view_draw,0);
664:   }
665:   PetscOptionsEnd();

667:   FormInitial_Coupled(user,X);
668:   VecView(X, PETSC_VIEWER_STDOUT_WORLD);

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

713:     PetscOptionsGetInt(0,"-col",&col,0);
714:     PetscOptionsGetBool(0,"-mult_dup",&mult_dup,0);
715:     PetscOptionsGetBool(0,"-view_dup",&view_dup,0);

717:     VecDuplicate(X,&Y);
718:     /* MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); */
719:     /* MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); */
720:     MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&D);
721:     VecZeroEntries(X);
722:     VecSetValue(X,col,1.0,INSERT_VALUES);
723:     VecAssemblyBegin(X);
724:     VecAssemblyEnd(X);
725:     MatMult(mult_dup?D:B,X,Y);
726:     MatView(view_dup?D:B,PETSC_VIEWER_STDOUT_WORLD);
727:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
728:     VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
729:     MatDestroy(&D);
730:     VecDestroy(&Y);
731:   }

733:   DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc);
734:   PetscFree(user);

736:   ISDestroy(&isg[0]);
737:   ISDestroy(&isg[1]);
738:   PetscFree(isg);
739:   VecDestroy(&X);
740:   VecDestroy(&F);
741:   MatDestroy(&B);
742:   DMDestroy(&dmu);
743:   DMDestroy(&dmk);
744:   DMDestroy(&pack);
745:   SNESDestroy(&snes);
746:   PetscFree(user);
747:   PetscFinalize();
748:   return 0;
749: }
```