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;
147: const PetscScalar g = 1.0 + gradu*gradu;
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", §ionU);
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", §ionK);
181: SectionRealCreateLocalVector(sectionK, &vecK);
182: VecCopy(user->Kloc, vecK);
183: VecDestroy(&vecK);
185: SectionRealDuplicate(sectionU, §ionFu);
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(§ionFu);
193: break;
194: case 1:
195: DMMeshGetSectionReal(dmk, "default", §ionK);
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", §ionU);
202: SectionRealCreateLocalVector(sectionU, &vecU);
203: VecCopy(user->Uloc, vecU);
204: VecDestroy(&vecU);
206: SectionRealDuplicate(sectionK, §ionFk);
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(§ionFk);
214: break;
215: case 2:
216: DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc);
217: DMCompositeScatter(user->pack, X, Uloc, Kloc);
219: DMMeshGetSectionReal(dmu, "default", §ionU);
220: SectionRealCreateLocalVector(sectionU, &vecU);
221: VecCopy(Uloc, vecU);
222: VecDestroy(&vecU);
224: DMMeshGetSectionReal(dmk, "default", §ionK);
225: SectionRealCreateLocalVector(sectionK, &vecK);
226: VecCopy(Kloc, vecK);
227: VecDestroy(&vecK);
229: DMCompositeRestoreLocalVectors(user->pack, &Uloc, &Kloc);
230: SectionRealDuplicate(sectionU, §ionFu);
231: SectionRealDuplicate(sectionK, §ionFk);
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(§ionFu);
244: SectionRealDestroy(§ionFk);
245: break;
246: }
247: SectionRealDestroy(§ionU);
248: SectionRealDestroy(§ionK);
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,
412: gradu = (ur[0]-u[0])/hx,
413: gradu_L = -1.0/hx,
414: gradu_R = 1.0/hx,
415: g = 1.0 + PetscSqr(gradu),
416: g_gradu = 2.0*gradu,
417: w = 1.0/(1.0+ubar) + 1.0/g,
418: w_ubar = -1./PetscSqr(1.+ubar),
419: w_gradu = -g_gradu/PetscSqr(g),
420: iw = 1.0/w,
421: iw_ubar = -w_ubar * PetscSqr(iw),
422: iw_gradu = -w_gradu * PetscSqr(iw);
424: values[0] = -hx*(iw_ubar*ubar_L + iw_gradu*gradu_L);
425: values[1] = -hx*(iw_ubar*ubar_R + iw_gradu*gradu_R);
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", §ionU);
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", §ionK);
453: SectionRealCreateLocalVector(sectionK, &vecK);
454: VecCopy(user->Kloc, vecK);
455: VecDestroy(&vecK);
457: FormJacobianLocal_U(dmu, dmk, sectionU, sectionK, *B, user);
458: SectionRealDestroy(§ionU);
459: SectionRealDestroy(§ionK);
460: break;
461: case 1:
462: DMMeshGetSectionReal(dmk, "default", §ionK);
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", §ionU);
469: SectionRealCreateLocalVector(sectionU, &vecU);
470: VecCopy(user->Uloc, vecU);
471: VecDestroy(&vecU);
473: FormJacobianLocal_K(dmu, dmk, sectionU, sectionK, *B, user);
474: SectionRealDestroy(§ionU);
475: SectionRealDestroy(§ionK);
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", §ionU);
485: SectionRealCreateLocalVector(sectionU, &vecU);
486: VecCopy(Uloc, vecU);
487: VecDestroy(&vecU);
489: DMMeshGetSectionReal(dmk, "default", §ionK);
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(§ionU);
511: SectionRealDestroy(§ionK);
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", §ionU);
547: DMMeshGetSectionReal(dmk, "default", §ionK);
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(§ionU);
576: SectionRealDestroy(§ionK);
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);
647: DMCompositeAddDM(pack, dmu);
648: DMCompositeAddDM(pack, dmk);
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: }