Actual source code: ex10.c
petsc-3.4.5 2014-06-29
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;
149: const PetscScalar g = 1.0 + gradu*gradu;
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", §ionU);
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", §ionK);
183: SectionRealCreateLocalVector(sectionK, &vecK);
184: VecCopy(user->Kloc, vecK);
185: VecDestroy(&vecK);
187: SectionRealDuplicate(sectionU, §ionFu);
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(§ionFu);
195: break;
196: case 1:
197: DMMeshGetSectionReal(dmk, "default", §ionK);
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", §ionU);
204: SectionRealCreateLocalVector(sectionU, &vecU);
205: VecCopy(user->Uloc, vecU);
206: VecDestroy(&vecU);
208: SectionRealDuplicate(sectionK, §ionFk);
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(§ionFk);
216: break;
217: case 2:
218: DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc);
219: DMCompositeScatter(user->pack, X, Uloc, Kloc);
221: DMMeshGetSectionReal(dmu, "default", §ionU);
222: SectionRealCreateLocalVector(sectionU, &vecU);
223: VecCopy(Uloc, vecU);
224: VecDestroy(&vecU);
226: DMMeshGetSectionReal(dmk, "default", §ionK);
227: SectionRealCreateLocalVector(sectionK, &vecK);
228: VecCopy(Kloc, vecK);
229: VecDestroy(&vecK);
231: DMCompositeRestoreLocalVectors(user->pack, &Uloc, &Kloc);
232: SectionRealDuplicate(sectionU, §ionFu);
233: SectionRealDuplicate(sectionK, §ionFk);
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(§ionFu);
246: SectionRealDestroy(§ionFk);
247: break;
248: }
249: SectionRealDestroy(§ionU);
250: SectionRealDestroy(§ionK);
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,
418: gradu = (ur[0]-u[0])/hx,
419: gradu_L = -1.0/hx,
420: gradu_R = 1.0/hx,
421: g = 1.0 + PetscSqr(gradu),
422: g_gradu = 2.0*gradu,
423: w = 1.0/(1.0+ubar) + 1.0/g,
424: w_ubar = -1./PetscSqr(1.+ubar),
425: w_gradu = -g_gradu/PetscSqr(g),
426: iw = 1.0/w,
427: iw_ubar = -w_ubar * PetscSqr(iw),
428: iw_gradu = -w_gradu * PetscSqr(iw);
430: values[0] = -hx*(iw_ubar*ubar_L + iw_gradu*gradu_L);
431: values[1] = -hx*(iw_ubar*ubar_R + iw_gradu*gradu_R);
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", §ionU);
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", §ionK);
459: SectionRealCreateLocalVector(sectionK, &vecK);
460: VecCopy(user->Kloc, vecK);
461: VecDestroy(&vecK);
463: FormJacobianLocal_U(dmu, dmk, sectionU, sectionK, *B, user);
464: SectionRealDestroy(§ionU);
465: SectionRealDestroy(§ionK);
466: break;
467: case 1:
468: DMMeshGetSectionReal(dmk, "default", §ionK);
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", §ionU);
475: SectionRealCreateLocalVector(sectionU, &vecU);
476: VecCopy(user->Uloc, vecU);
477: VecDestroy(&vecU);
479: FormJacobianLocal_K(dmu, dmk, sectionU, sectionK, *B, user);
480: SectionRealDestroy(§ionU);
481: SectionRealDestroy(§ionK);
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", §ionU);
491: SectionRealCreateLocalVector(sectionU, &vecU);
492: VecCopy(Uloc, vecU);
493: VecDestroy(&vecU);
495: DMMeshGetSectionReal(dmk, "default", §ionK);
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(§ionU);
517: SectionRealDestroy(§ionK);
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", §ionU);
553: DMMeshGetSectionReal(dmk, "default", §ionK);
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(§ionU);
582: SectionRealDestroy(§ionK);
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);
653: DMCompositeAddDM(pack, dmu);
654: DMCompositeAddDM(pack, dmk);
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: }