Actual source code: ex28.c
petsc-3.3-p7 2013-05-11
1: static const char help[] = "1D multiphysics prototype with analytic Jacobians to solve individual problems and a coupled problem.\n\n";
3: /* Solve a PDE coupled to an algebraic system in 1D
4: *
5: * PDE (U):
6: * -(k u_x)_x = 1 on (0,1), subject to u(0) = 0, u(1) = 1
7: * Algebraic (K):
8: * exp(k-1) + k = u + 1/(1/(1+u) + 1/(1+u_x^2))
9: *
10: * The discretization places k at staggered points, and a separate DMDA is used for each "physics".
11: *
12: * This example is a prototype for coupling in multi-physics problems, therefore residual evaluation and assembly for
13: * each problem (referred to as U and K) are written separately. This permits the same "physics" code to be used for
14: * solving each uncoupled problem as well as the coupled system. In particular, run with -problem_type 0 to solve only
15: * problem U (with K fixed), -problem_type 1 to solve only K (with U fixed), and -problem_type 2 to solve both at once.
16: *
17: * In all cases, a fully-assembled analytic Jacobian is available, so the systems can be solved with a direct solve or
18: * any other standard method. Additionally, by running with
19: *
20: * -pack_dm_mat_type nest
21: *
22: * The same code assembles a coupled matrix where each block is stored separately, which allows the use of PCFieldSplit
23: * without copying values to extract submatrices.
24: */
26: #include <petscsnes.h>
27: #include <petscdmda.h>
28: #include <petscdmcomposite.h>
30: typedef struct _UserCtx *User;
31: struct _UserCtx {
32: PetscInt ptype;
33: DM pack;
34: Vec Uloc,Kloc;
35: };
39: static PetscErrorCode FormFunctionLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])
40: {
41: PetscReal hx = 1./info->mx;
42: PetscInt i;
45: for (i=info->xs; i<info->xs+info->xm; i++) {
46: if (i == 0) f[i] = 1./hx*u[i];
47: else if (i == info->mx-1) f[i] = 1./hx*(u[i] - 1.0);
48: else f[i] = hx*((k[i-1]*(u[i]-u[i-1]) - k[i]*(u[i+1]-u[i]))/(hx*hx) - 1.0);
49: }
50: return(0);
51: }
55: static PetscErrorCode FormFunctionLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])
56: {
57: PetscReal hx = 1./info->mx;
58: PetscInt i;
61: for (i=info->xs; i<info->xs+info->xm; i++) {
62: const PetscScalar
63: ubar = 0.5*(u[i+1]+u[i]),
64: gradu = (u[i+1]-u[i])/hx,
65: g = 1. + gradu*gradu,
66: w = 1./(1.+ubar) + 1./g;
67: f[i] = hx*(PetscExpScalar(k[i]-1.0) + k[i] - 1./w);
68: }
69: return(0);
70: }
74: static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx)
75: {
76: User user = (User)ctx;
77: DM dau,dak;
78: DMDALocalInfo infou,infok;
79: PetscScalar *u,*k;
80: PetscScalar *fu,*fk;
81: PetscErrorCode ierr;
82: Vec Uloc,Kloc,Fu,Fk;
85: DMCompositeGetEntries(user->pack,&dau,&dak);
86: DMDAGetLocalInfo(dau,&infou);
87: DMDAGetLocalInfo(dak,&infok);
88: DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
89: switch (user->ptype) {
90: case 0:
91: DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
92: DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc);
93: DMDAVecGetArray(dau,Uloc,&u);
94: DMDAVecGetArray(dak,user->Kloc,&k);
95: DMDAVecGetArray(dau,F,&fu);
96: FormFunctionLocal_U(user,&infou,u,k,fu);
97: DMDAVecRestoreArray(dau,F,&fu);
98: DMDAVecRestoreArray(dau,Uloc,&u);
99: DMDAVecRestoreArray(dak,user->Kloc,&k);
100: break;
101: case 1:
102: DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
103: DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc);
104: DMDAVecGetArray(dau,user->Uloc,&u);
105: DMDAVecGetArray(dak,Kloc,&k);
106: DMDAVecGetArray(dak,F,&fk);
107: FormFunctionLocal_K(user,&infok,u,k,fk);
108: DMDAVecRestoreArray(dak,F,&fk);
109: DMDAVecRestoreArray(dau,user->Uloc,&u);
110: DMDAVecRestoreArray(dak,Kloc,&k);
111: break;
112: case 2:
113: DMCompositeScatter(user->pack,X,Uloc,Kloc);
114: DMDAVecGetArray(dau,Uloc,&u);
115: DMDAVecGetArray(dak,Kloc,&k);
116: DMCompositeGetAccess(user->pack,F,&Fu,&Fk);
117: DMDAVecGetArray(dau,Fu,&fu);
118: DMDAVecGetArray(dak,Fk,&fk);
119: FormFunctionLocal_U(user,&infou,u,k,fu);
120: FormFunctionLocal_K(user,&infok,u,k,fk);
121: DMDAVecRestoreArray(dau,Fu,&fu);
122: DMDAVecRestoreArray(dak,Fk,&fk);
123: DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk);
124: DMDAVecRestoreArray(dau,Uloc,&u);
125: DMDAVecRestoreArray(dak,Kloc,&k);
126: break;
127: }
128: DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
129: return(0);
130: }
134: static PetscErrorCode FormJacobianLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Buu)
135: {
136: PetscReal hx = 1./info->mx;
138: PetscInt i;
141: for (i=info->xs; i<info->xs+info->xm; i++) {
142: PetscInt row = i-info->gxs,cols[] = {row-1,row,row+1};
143: PetscScalar val = 1./hx;
144: if (i == 0) {MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);}
145: else if (i == info->mx-1) {MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);}
146: else {
147: PetscScalar vals[] = {-k[i-1]/hx,(k[i-1]+k[i])/hx,-k[i]/hx};
148: MatSetValuesLocal(Buu,1,&row,3,cols,vals,INSERT_VALUES);
149: }
150: }
151: return(0);
152: }
156: static PetscErrorCode FormJacobianLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Bkk)
157: {
158: PetscReal hx = 1./info->mx;
160: PetscInt i;
163: for (i=info->xs; i<info->xs+info->xm; i++) {
164: PetscInt row = i-info->gxs;
165: PetscScalar vals[] = {hx*(PetscExpScalar(k[i]-1.)+1.)};
166: MatSetValuesLocal(Bkk,1,&row,1,&row,vals,INSERT_VALUES);
167: }
168: return(0);
169: }
173: static PetscErrorCode FormJacobianLocal_UK(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Buk)
174: {
175: PetscReal hx = 1./info->mx;
177: PetscInt i;
178: PetscInt row,cols[2];
179: PetscScalar vals[2];
182: if (!Buk) return(0); /* Not assembling this block */
183: for (i=info->xs; i<info->xs+info->xm; i++) {
184: if (i == 0 || i == info->mx-1) continue;
185: row = i-info->gxs;
186: cols[0] = i-1-infok->gxs; vals[0] = (u[i]-u[i-1])/hx;
187: cols[1] = i-infok->gxs; vals[1] = (u[i]-u[i+1])/hx;
188: MatSetValuesLocal(Buk,1,&row,2,cols,vals,INSERT_VALUES);
189: }
190: return(0);
191: }
195: static PetscErrorCode FormJacobianLocal_KU(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Bku)
196: {
198: PetscInt i;
199: PetscReal hx = 1./(info->mx-1);
202: if (!Bku) return(0); /* Not assembling this block */
203: for (i=infok->xs; i<infok->xs+infok->xm; i++) {
204: PetscInt row = i-infok->gxs,cols[2];
205: PetscScalar vals[2];
206: const PetscScalar
207: ubar = 0.5*(u[i]+u[i+1]),
208: ubar_L = 0.5,
209: ubar_R = 0.5,
210: gradu = (u[i+1]-u[i])/hx,
211: gradu_L = -1./hx,
212: gradu_R = 1./hx,
213: g = 1. + PetscSqr(gradu),
214: g_gradu = 2.*gradu,
215: w = 1./(1.+ubar) + 1./g,
216: w_ubar = -1./PetscSqr(1.+ubar),
217: w_gradu = -g_gradu/PetscSqr(g),
218: iw = 1./w,
219: iw_ubar = -w_ubar * PetscSqr(iw),
220: iw_gradu = -w_gradu * PetscSqr(iw);
221: cols[0] = i-info->gxs; vals[0] = -hx*(iw_ubar*ubar_L + iw_gradu*gradu_L);
222: cols[1] = i+1-info->gxs; vals[1] = -hx*(iw_ubar*ubar_R + iw_gradu*gradu_R);
223: MatSetValuesLocal(Bku,1,&row,2,cols,vals,INSERT_VALUES);
224: }
225: return(0);
226: }
230: static PetscErrorCode FormJacobian_All(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *mstr,void *ctx)
231: {
232: User user = (User)ctx;
233: DM dau,dak;
234: DMDALocalInfo infou,infok;
235: PetscScalar *u,*k;
236: PetscErrorCode ierr;
237: Vec Uloc,Kloc;
240: DMCompositeGetEntries(user->pack,&dau,&dak);
241: DMDAGetLocalInfo(dau,&infou);
242: DMDAGetLocalInfo(dak,&infok);
243: DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
244: switch (user->ptype) {
245: case 0:
246: DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
247: DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc);
248: DMDAVecGetArray(dau,Uloc,&u);
249: DMDAVecGetArray(dak,user->Kloc,&k);
250: FormJacobianLocal_U(user,&infou,u,k,*B);
251: DMDAVecRestoreArray(dau,Uloc,&u);
252: DMDAVecRestoreArray(dak,user->Kloc,&k);
253: break;
254: case 1:
255: DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
256: DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc);
257: DMDAVecGetArray(dau,user->Uloc,&u);
258: DMDAVecGetArray(dak,Kloc,&k);
259: FormJacobianLocal_K(user,&infok,u,k,*B);
260: DMDAVecRestoreArray(dau,user->Uloc,&u);
261: DMDAVecRestoreArray(dak,Kloc,&k);
262: break;
263: case 2: {
264: Mat Buu,Buk,Bku,Bkk;
265: IS *is;
266: DMCompositeScatter(user->pack,X,Uloc,Kloc);
267: DMDAVecGetArray(dau,Uloc,&u);
268: DMDAVecGetArray(dak,Kloc,&k);
269: DMCompositeGetLocalISs(user->pack,&is);
270: MatGetLocalSubMatrix(*B,is[0],is[0],&Buu);
271: MatGetLocalSubMatrix(*B,is[0],is[1],&Buk);
272: MatGetLocalSubMatrix(*B,is[1],is[0],&Bku);
273: MatGetLocalSubMatrix(*B,is[1],is[1],&Bkk);
274: FormJacobianLocal_U(user,&infou,u,k,Buu);
275: FormJacobianLocal_UK(user,&infou,&infok,u,k,Buk);
276: FormJacobianLocal_KU(user,&infou,&infok,u,k,Bku);
277: FormJacobianLocal_K(user,&infok,u,k,Bkk);
278: MatRestoreLocalSubMatrix(*B,is[0],is[0],&Buu);
279: MatRestoreLocalSubMatrix(*B,is[0],is[1],&Buk);
280: MatRestoreLocalSubMatrix(*B,is[1],is[0],&Bku);
281: MatRestoreLocalSubMatrix(*B,is[1],is[1],&Bkk);
282: DMDAVecRestoreArray(dau,Uloc,&u);
283: DMDAVecRestoreArray(dak,Kloc,&k);
285: ISDestroy(&is[0]);
286: ISDestroy(&is[1]);
287: PetscFree(is);
288: } break;
289: }
290: DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
291: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
292: MatAssemblyEnd (*B,MAT_FINAL_ASSEMBLY);
293: if (*J != *B) {
294: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
295: MatAssemblyEnd (*J,MAT_FINAL_ASSEMBLY);
296: }
297: *mstr = DIFFERENT_NONZERO_PATTERN;
298: return(0);
299: }
303: static PetscErrorCode FormInitial_Coupled(User user,Vec X)
304: {
306: DM dau,dak;
307: DMDALocalInfo infou,infok;
308: Vec Xu,Xk;
309: PetscScalar *u,*k,hx;
310: PetscInt i;
313: DMCompositeGetEntries(user->pack,&dau,&dak);
314: DMCompositeGetAccess(user->pack,X,&Xu,&Xk);
315: DMDAVecGetArray(dau,Xu,&u);
316: DMDAVecGetArray(dak,Xk,&k);
317: DMDAGetLocalInfo(dau,&infou);
318: DMDAGetLocalInfo(dak,&infok);
319: hx = 1./(infok.mx);
320: for (i=infou.xs; i<infou.xs+infou.xm; i++) u[i] = (PetscScalar)i*hx * (1.-(PetscScalar)i*hx);
321: for (i=infok.xs; i<infok.xs+infok.xm; i++) k[i] = 1.0 + 0.5*(PetscScalar)sin((double)2*PETSC_PI*i*hx);
322: DMDAVecRestoreArray(dau,Xu,&u);
323: DMDAVecRestoreArray(dak,Xk,&k);
324: DMCompositeRestoreAccess(user->pack,X,&Xu,&Xk);
325: DMCompositeScatter(user->pack,X,user->Uloc,user->Kloc);
326: return(0);
327: }
331: int main(int argc, char *argv[])
332: {
334: DM dau,dak,pack;
335: const PetscInt *lxu;
336: PetscInt *lxk,m,nprocs;
337: User user;
338: SNES snes;
339: Vec X,F,Xu,Xk,Fu,Fk;
340: Mat B;
341: IS *isg;
342: PetscBool view_draw,pass_dm;
344: PetscInitialize(&argc,&argv,0,help);
345: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-10,1,1,PETSC_NULL,&dau);
346: DMSetOptionsPrefix(dau,"u_");
347: DMSetFromOptions(dau);
348: DMDAGetOwnershipRanges(dau,&lxu,0,0);
349: DMDAGetInfo(dau,0, &m,0,0, &nprocs,0,0, 0,0,0,0,0,0);
350: PetscMalloc(nprocs*sizeof(*lxk),&lxk);
351: PetscMemcpy(lxk,lxu,nprocs*sizeof(*lxk));
352: lxk[0]--;
353: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,m-1,1,1,lxk,&dak);
354: DMSetOptionsPrefix(dak,"k_");
355: DMSetFromOptions(dau);
356: PetscFree(lxk);
358: DMCompositeCreate(PETSC_COMM_WORLD,&pack);
359: DMSetOptionsPrefix(pack,"pack_");
360: DMCompositeAddDM(pack,dau);
361: DMCompositeAddDM(pack,dak);
362: DMDASetFieldName(dau,0,"u");
363: DMDASetFieldName(dak,0,"k");
364: DMSetFromOptions(pack);
366: DMCreateGlobalVector(pack,&X);
367: VecDuplicate(X,&F);
369: PetscNew(struct _UserCtx,&user);
370: user->pack = pack;
371: DMCompositeGetGlobalISs(pack,&isg);
372: DMCompositeGetLocalVectors(pack,&user->Uloc,&user->Kloc);
373: DMCompositeScatter(pack,X,user->Uloc,user->Kloc);
375: PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Coupled problem options","SNES");
376: {
377: user->ptype = 0; view_draw = PETSC_FALSE; pass_dm = PETSC_TRUE;
378: PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,0);
379: PetscOptionsBool("-view_draw","Draw the final coupled solution regardless of whether only one physics was solved",0,view_draw,&view_draw,0);
380: PetscOptionsBool("-pass_dm","Pass the packed DM to SNES to use when determining splits and forward into splits",0,pass_dm,&pass_dm,0);
381: }
382: PetscOptionsEnd();
384: FormInitial_Coupled(user,X);
386: SNESCreate(PETSC_COMM_WORLD,&snes);
387: switch (user->ptype) {
388: case 0:
389: DMCompositeGetAccess(pack,X,&Xu,0);
390: DMCompositeGetAccess(pack,F,&Fu,0);
391: DMCreateMatrix(dau,PETSC_NULL,&B);
392: SNESSetFunction(snes,Fu,FormFunction_All,user);
393: SNESSetJacobian(snes,B,B,FormJacobian_All,user);
394: SNESSetFromOptions(snes);
395: SNESSetDM(snes,dau);
396: SNESSolve(snes,PETSC_NULL,Xu);
397: DMCompositeRestoreAccess(pack,X,&Xu,0);
398: DMCompositeRestoreAccess(pack,F,&Fu,0);
399: break;
400: case 1:
401: DMCompositeGetAccess(pack,X,0,&Xk);
402: DMCompositeGetAccess(pack,F,0,&Fk);
403: DMCreateMatrix(dak,PETSC_NULL,&B);
404: SNESSetFunction(snes,Fk,FormFunction_All,user);
405: SNESSetJacobian(snes,B,B,FormJacobian_All,user);
406: SNESSetFromOptions(snes);
407: SNESSetDM(snes,dak);
408: SNESSolve(snes,PETSC_NULL,Xk);
409: DMCompositeRestoreAccess(pack,X,0,&Xk);
410: DMCompositeRestoreAccess(pack,F,0,&Fk);
411: break;
412: case 2:
413: DMCreateMatrix(pack,PETSC_NULL,&B);
414: /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */
415: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
416: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
417: SNESSetFunction(snes,F,FormFunction_All,user);
418: SNESSetJacobian(snes,B,B,FormJacobian_All,user);
419: SNESSetFromOptions(snes);
420: if (!pass_dm) { /* Manually provide index sets and names for the splits */
421: KSP ksp;
422: PC pc;
423: SNESGetKSP(snes,&ksp);
424: KSPGetPC(ksp,&pc);
425: PCFieldSplitSetIS(pc,"u",isg[0]);
426: PCFieldSplitSetIS(pc,"k",isg[1]);
427: } else {
428: /* The same names come from the options prefix for dau and dak. This option can support geometric multigrid inside
429: * of splits, but it requires using a DM (perhaps your own implementation). */
430: SNESSetDM(snes,pack);
431: }
432: SNESSolve(snes,PETSC_NULL,X);
433: break;
434: }
435: if (view_draw) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
436: if (0) {
437: PetscInt col = 0;
438: PetscBool mult_dup = PETSC_FALSE,view_dup = PETSC_FALSE;
439: Mat D;
440: Vec Y;
442: PetscOptionsGetInt(0,"-col",&col,0);
443: PetscOptionsGetBool(0,"-mult_dup",&mult_dup,0);
444: PetscOptionsGetBool(0,"-view_dup",&view_dup,0);
446: VecDuplicate(X,&Y);
447: /* MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); */
448: /* MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); */
449: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&D);
450: VecZeroEntries(X);
451: VecSetValue(X,col,1.0,INSERT_VALUES);
452: VecAssemblyBegin(X);
453: VecAssemblyEnd(X);
454: MatMult(mult_dup?D:B,X,Y);
455: MatView(view_dup?D:B,PETSC_VIEWER_STDOUT_WORLD);
456: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
457: VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
458: MatDestroy(&D);
459: VecDestroy(&Y);
460: }
462: DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc);
463: PetscFree(user);
465: ISDestroy(&isg[0]);
466: ISDestroy(&isg[1]);
467: PetscFree(isg);
468: VecDestroy(&X);
469: VecDestroy(&F);
470: MatDestroy(&B);
471: DMDestroy(&dau);
472: DMDestroy(&dak);
473: DMDestroy(&pack);
474: SNESDestroy(&snes);
475: PetscFinalize();
476: return 0;
477: }