Actual source code: ex28.c
petsc-3.7.7 2017-09-25
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 = 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 <petscdm.h>
28: #include <petscdmda.h>
29: #include <petscdmcomposite.h>
31: typedef struct _UserCtx *User;
32: struct _UserCtx {
33: PetscInt ptype;
34: DM pack;
35: Vec Uloc,Kloc;
36: };
40: static PetscErrorCode FormFunctionLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])
41: {
42: PetscReal hx = 1./info->mx;
43: PetscInt i;
46: for (i=info->xs; i<info->xs+info->xm; i++) {
47: if (i == 0) f[i] = 1./hx*u[i];
48: else if (i == info->mx-1) f[i] = 1./hx*(u[i] - 1.0);
49: else f[i] = hx*((k[i-1]*(u[i]-u[i-1]) - k[i]*(u[i+1]-u[i]))/(hx*hx) - 1.0);
50: }
51: return(0);
52: }
56: static PetscErrorCode FormFunctionLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])
57: {
58: PetscReal hx = 1./info->mx;
59: PetscInt i;
62: for (i=info->xs; i<info->xs+info->xm; i++) {
63: const PetscScalar
64: ubar = 0.5*(u[i+1]+u[i]),
65: gradu = (u[i+1]-u[i])/hx,
66: g = 1. + gradu*gradu,
67: w = 1./(1.+ubar) + 1./g;
68: f[i] = hx*(PetscExpScalar(k[i]-1.0) + k[i] - 1./w);
69: }
70: return(0);
71: }
75: static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx)
76: {
77: User user = (User)ctx;
78: DM dau,dak;
79: DMDALocalInfo infou,infok;
80: PetscScalar *u,*k;
81: PetscScalar *fu,*fk;
83: Vec Uloc,Kloc,Fu,Fk;
86: DMCompositeGetEntries(user->pack,&dau,&dak);
87: DMDAGetLocalInfo(dau,&infou);
88: DMDAGetLocalInfo(dak,&infok);
89: DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
90: switch (user->ptype) {
91: case 0:
92: DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
93: DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc);
94: DMDAVecGetArray(dau,Uloc,&u);
95: DMDAVecGetArray(dak,user->Kloc,&k);
96: DMDAVecGetArray(dau,F,&fu);
97: FormFunctionLocal_U(user,&infou,u,k,fu);
98: DMDAVecRestoreArray(dau,F,&fu);
99: DMDAVecRestoreArray(dau,Uloc,&u);
100: DMDAVecRestoreArray(dak,user->Kloc,&k);
101: break;
102: case 1:
103: DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
104: DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc);
105: DMDAVecGetArray(dau,user->Uloc,&u);
106: DMDAVecGetArray(dak,Kloc,&k);
107: DMDAVecGetArray(dak,F,&fk);
108: FormFunctionLocal_K(user,&infok,u,k,fk);
109: DMDAVecRestoreArray(dak,F,&fk);
110: DMDAVecRestoreArray(dau,user->Uloc,&u);
111: DMDAVecRestoreArray(dak,Kloc,&k);
112: break;
113: case 2:
114: DMCompositeScatter(user->pack,X,Uloc,Kloc);
115: DMDAVecGetArray(dau,Uloc,&u);
116: DMDAVecGetArray(dak,Kloc,&k);
117: DMCompositeGetAccess(user->pack,F,&Fu,&Fk);
118: DMDAVecGetArray(dau,Fu,&fu);
119: DMDAVecGetArray(dak,Fk,&fk);
120: FormFunctionLocal_U(user,&infou,u,k,fu);
121: FormFunctionLocal_K(user,&infok,u,k,fk);
122: DMDAVecRestoreArray(dau,Fu,&fu);
123: DMDAVecRestoreArray(dak,Fk,&fk);
124: DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk);
125: DMDAVecRestoreArray(dau,Uloc,&u);
126: DMDAVecRestoreArray(dak,Kloc,&k);
127: break;
128: }
129: DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
130: return(0);
131: }
135: static PetscErrorCode FormJacobianLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Buu)
136: {
137: PetscReal hx = 1./info->mx;
139: PetscInt i;
142: for (i=info->xs; i<info->xs+info->xm; i++) {
143: PetscInt row = i-info->gxs,cols[] = {row-1,row,row+1};
144: PetscScalar val = 1./hx;
145: if (i == 0) {
146: MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);
147: } else if (i == info->mx-1) {
148: MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);
149: } else {
150: PetscScalar vals[] = {-k[i-1]/hx,(k[i-1]+k[i])/hx,-k[i]/hx};
151: MatSetValuesLocal(Buu,1,&row,3,cols,vals,INSERT_VALUES);
152: }
153: }
154: return(0);
155: }
159: static PetscErrorCode FormJacobianLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Bkk)
160: {
161: PetscReal hx = 1./info->mx;
163: PetscInt i;
166: for (i=info->xs; i<info->xs+info->xm; i++) {
167: PetscInt row = i-info->gxs;
168: PetscScalar vals[] = {hx*(PetscExpScalar(k[i]-1.)+1.)};
169: MatSetValuesLocal(Bkk,1,&row,1,&row,vals,INSERT_VALUES);
170: }
171: return(0);
172: }
176: static PetscErrorCode FormJacobianLocal_UK(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Buk)
177: {
178: PetscReal hx = 1./info->mx;
180: PetscInt i;
181: PetscInt row,cols[2];
182: PetscScalar vals[2];
185: if (!Buk) return(0); /* Not assembling this block */
186: for (i=info->xs; i<info->xs+info->xm; i++) {
187: if (i == 0 || i == info->mx-1) continue;
188: row = i-info->gxs;
189: cols[0] = i-1-infok->gxs; vals[0] = (u[i]-u[i-1])/hx;
190: cols[1] = i-infok->gxs; vals[1] = (u[i]-u[i+1])/hx;
191: MatSetValuesLocal(Buk,1,&row,2,cols,vals,INSERT_VALUES);
192: }
193: return(0);
194: }
198: static PetscErrorCode FormJacobianLocal_KU(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Bku)
199: {
201: PetscInt i;
202: PetscReal hx = 1./(info->mx-1);
205: if (!Bku) return(0); /* Not assembling this block */
206: for (i=infok->xs; i<infok->xs+infok->xm; i++) {
207: PetscInt row = i-infok->gxs,cols[2];
208: PetscScalar vals[2];
209: const PetscScalar
210: ubar = 0.5*(u[i]+u[i+1]),
211: ubar_L = 0.5,
212: ubar_R = 0.5,
213: gradu = (u[i+1]-u[i])/hx,
214: gradu_L = -1./hx,
215: gradu_R = 1./hx,
216: g = 1. + PetscSqr(gradu),
217: g_gradu = 2.*gradu,
218: w = 1./(1.+ubar) + 1./g,
219: w_ubar = -1./PetscSqr(1.+ubar),
220: w_gradu = -g_gradu/PetscSqr(g),
221: iw = 1./w,
222: iw_ubar = -w_ubar * PetscSqr(iw),
223: iw_gradu = -w_gradu * PetscSqr(iw);
224: cols[0] = i-info->gxs; vals[0] = -hx*(iw_ubar*ubar_L + iw_gradu*gradu_L);
225: cols[1] = i+1-info->gxs; vals[1] = -hx*(iw_ubar*ubar_R + iw_gradu*gradu_R);
226: MatSetValuesLocal(Bku,1,&row,2,cols,vals,INSERT_VALUES);
227: }
228: return(0);
229: }
233: static PetscErrorCode FormJacobian_All(SNES snes,Vec X,Mat J,Mat B,void *ctx)
234: {
235: User user = (User)ctx;
236: DM dau,dak;
237: DMDALocalInfo infou,infok;
238: PetscScalar *u,*k;
240: Vec Uloc,Kloc;
243: DMCompositeGetEntries(user->pack,&dau,&dak);
244: DMDAGetLocalInfo(dau,&infou);
245: DMDAGetLocalInfo(dak,&infok);
246: DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
247: switch (user->ptype) {
248: case 0:
249: DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
250: DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc);
251: DMDAVecGetArray(dau,Uloc,&u);
252: DMDAVecGetArray(dak,user->Kloc,&k);
253: FormJacobianLocal_U(user,&infou,u,k,B);
254: DMDAVecRestoreArray(dau,Uloc,&u);
255: DMDAVecRestoreArray(dak,user->Kloc,&k);
256: break;
257: case 1:
258: DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
259: DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc);
260: DMDAVecGetArray(dau,user->Uloc,&u);
261: DMDAVecGetArray(dak,Kloc,&k);
262: FormJacobianLocal_K(user,&infok,u,k,B);
263: DMDAVecRestoreArray(dau,user->Uloc,&u);
264: DMDAVecRestoreArray(dak,Kloc,&k);
265: break;
266: case 2: {
267: Mat Buu,Buk,Bku,Bkk;
268: PetscBool nest;
269: IS *is;
270: DMCompositeScatter(user->pack,X,Uloc,Kloc);
271: DMDAVecGetArray(dau,Uloc,&u);
272: DMDAVecGetArray(dak,Kloc,&k);
273: DMCompositeGetLocalISs(user->pack,&is);
274: MatGetLocalSubMatrix(B,is[0],is[0],&Buu);
275: MatGetLocalSubMatrix(B,is[0],is[1],&Buk);
276: MatGetLocalSubMatrix(B,is[1],is[0],&Bku);
277: MatGetLocalSubMatrix(B,is[1],is[1],&Bkk);
278: FormJacobianLocal_U(user,&infou,u,k,Buu);
279: PetscObjectTypeCompare((PetscObject)B,MATNEST,&nest);
280: if (!nest) {
281: /*
282: DMCreateMatrix_Composite() for a nested matrix does not generate off-block matrices that one can call MatSetValuesLocal() on, it just creates dummy
283: matrices with no entries; there cannot insert values into them. Commit b6480e041dd2293a65f96222772d68cdb4ed6306
284: changed Mat_Nest() from returning NULL pointers for these submatrices to dummy matrices because PCFIELDSPLIT could not
285: handle the returned null matrices.
286: */
287: FormJacobianLocal_UK(user,&infou,&infok,u,k,Buk);
288: FormJacobianLocal_KU(user,&infou,&infok,u,k,Bku);
289: }
290: FormJacobianLocal_K(user,&infok,u,k,Bkk);
291: MatRestoreLocalSubMatrix(B,is[0],is[0],&Buu);
292: MatRestoreLocalSubMatrix(B,is[0],is[1],&Buk);
293: MatRestoreLocalSubMatrix(B,is[1],is[0],&Bku);
294: MatRestoreLocalSubMatrix(B,is[1],is[1],&Bkk);
295: DMDAVecRestoreArray(dau,Uloc,&u);
296: DMDAVecRestoreArray(dak,Kloc,&k);
298: ISDestroy(&is[0]);
299: ISDestroy(&is[1]);
300: PetscFree(is);
301: } break;
302: }
303: DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
304: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
305: MatAssemblyEnd (B,MAT_FINAL_ASSEMBLY);
306: if (J != B) {
307: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd (J,MAT_FINAL_ASSEMBLY);
309: }
310: return(0);
311: }
315: static PetscErrorCode FormInitial_Coupled(User user,Vec X)
316: {
318: DM dau,dak;
319: DMDALocalInfo infou,infok;
320: Vec Xu,Xk;
321: PetscScalar *u,*k,hx;
322: PetscInt i;
325: DMCompositeGetEntries(user->pack,&dau,&dak);
326: DMCompositeGetAccess(user->pack,X,&Xu,&Xk);
327: DMDAVecGetArray(dau,Xu,&u);
328: DMDAVecGetArray(dak,Xk,&k);
329: DMDAGetLocalInfo(dau,&infou);
330: DMDAGetLocalInfo(dak,&infok);
331: hx = 1./(infok.mx);
332: for (i=infou.xs; i<infou.xs+infou.xm; i++) u[i] = (PetscScalar)i*hx * (1.-(PetscScalar)i*hx);
333: for (i=infok.xs; i<infok.xs+infok.xm; i++) k[i] = 1.0 + 0.5*PetscSinScalar(2*PETSC_PI*i*hx);
334: DMDAVecRestoreArray(dau,Xu,&u);
335: DMDAVecRestoreArray(dak,Xk,&k);
336: DMCompositeRestoreAccess(user->pack,X,&Xu,&Xk);
337: DMCompositeScatter(user->pack,X,user->Uloc,user->Kloc);
338: return(0);
339: }
343: int main(int argc, char *argv[])
344: {
346: DM dau,dak,pack;
347: const PetscInt *lxu;
348: PetscInt *lxk,m,sizes;
349: User user;
350: SNES snes;
351: Vec X,F,Xu,Xk,Fu,Fk;
352: Mat B;
353: IS *isg;
354: PetscBool view_draw,pass_dm;
356: PetscInitialize(&argc,&argv,0,help);
357: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-10,1,1,NULL,&dau);
358: DMSetOptionsPrefix(dau,"u_");
359: DMSetFromOptions(dau);
360: DMDAGetOwnershipRanges(dau,&lxu,0,0);
361: DMDAGetInfo(dau,0, &m,0,0, &sizes,0,0, 0,0,0,0,0,0);
362: PetscMalloc1(sizes,&lxk);
363: PetscMemcpy(lxk,lxu,sizes*sizeof(*lxk));
364: lxk[0]--;
365: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m-1,1,1,lxk,&dak);
366: DMSetOptionsPrefix(dak,"k_");
367: DMSetFromOptions(dak);
368: PetscFree(lxk);
370: DMCompositeCreate(PETSC_COMM_WORLD,&pack);
371: DMSetOptionsPrefix(pack,"pack_");
372: DMCompositeAddDM(pack,dau);
373: DMCompositeAddDM(pack,dak);
374: DMDASetFieldName(dau,0,"u");
375: DMDASetFieldName(dak,0,"k");
376: DMSetFromOptions(pack);
378: DMCreateGlobalVector(pack,&X);
379: VecDuplicate(X,&F);
381: PetscNew(&user);
383: user->pack = pack;
385: DMCompositeGetGlobalISs(pack,&isg);
386: DMCompositeGetLocalVectors(pack,&user->Uloc,&user->Kloc);
387: DMCompositeScatter(pack,X,user->Uloc,user->Kloc);
389: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Coupled problem options","SNES");
390: {
391: user->ptype = 0; view_draw = PETSC_FALSE; pass_dm = PETSC_TRUE;
393: PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,NULL);
394: PetscOptionsBool("-view_draw","Draw the final coupled solution regardless of whether only one physics was solved",0,view_draw,&view_draw,NULL);
395: PetscOptionsBool("-pass_dm","Pass the packed DM to SNES to use when determining splits and forward into splits",0,pass_dm,&pass_dm,NULL);
396: }
397: PetscOptionsEnd();
399: FormInitial_Coupled(user,X);
401: SNESCreate(PETSC_COMM_WORLD,&snes);
402: switch (user->ptype) {
403: case 0:
404: DMCompositeGetAccess(pack,X,&Xu,0);
405: DMCompositeGetAccess(pack,F,&Fu,0);
406: DMCreateMatrix(dau,&B);
407: SNESSetFunction(snes,Fu,FormFunction_All,user);
408: SNESSetJacobian(snes,B,B,FormJacobian_All,user);
409: SNESSetFromOptions(snes);
410: SNESSetDM(snes,dau);
411: SNESSolve(snes,NULL,Xu);
412: DMCompositeRestoreAccess(pack,X,&Xu,0);
413: DMCompositeRestoreAccess(pack,F,&Fu,0);
414: break;
415: case 1:
416: DMCompositeGetAccess(pack,X,0,&Xk);
417: DMCompositeGetAccess(pack,F,0,&Fk);
418: DMCreateMatrix(dak,&B);
419: SNESSetFunction(snes,Fk,FormFunction_All,user);
420: SNESSetJacobian(snes,B,B,FormJacobian_All,user);
421: SNESSetFromOptions(snes);
422: SNESSetDM(snes,dak);
423: SNESSolve(snes,NULL,Xk);
424: DMCompositeRestoreAccess(pack,X,0,&Xk);
425: DMCompositeRestoreAccess(pack,F,0,&Fk);
426: break;
427: case 2:
428: DMCreateMatrix(pack,&B);
429: /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */
430: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
431: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
432: SNESSetFunction(snes,F,FormFunction_All,user);
433: SNESSetJacobian(snes,B,B,FormJacobian_All,user);
434: SNESSetFromOptions(snes);
435: if (!pass_dm) { /* Manually provide index sets and names for the splits */
436: KSP ksp;
437: PC pc;
438: SNESGetKSP(snes,&ksp);
439: KSPGetPC(ksp,&pc);
440: PCFieldSplitSetIS(pc,"u",isg[0]);
441: PCFieldSplitSetIS(pc,"k",isg[1]);
442: } else {
443: /* The same names come from the options prefix for dau and dak. This option can support geometric multigrid inside
444: * of splits, but it requires using a DM (perhaps your own implementation). */
445: SNESSetDM(snes,pack);
446: }
447: SNESSolve(snes,NULL,X);
448: break;
449: }
450: if (view_draw) {VecView(X,PETSC_VIEWER_DRAW_WORLD);}
451: if (0) {
452: PetscInt col = 0;
453: PetscBool mult_dup = PETSC_FALSE,view_dup = PETSC_FALSE;
454: Mat D;
455: Vec Y;
457: PetscOptionsGetInt(NULL,0,"-col",&col,0);
458: PetscOptionsGetBool(NULL,0,"-mult_dup",&mult_dup,0);
459: PetscOptionsGetBool(NULL,0,"-view_dup",&view_dup,0);
461: VecDuplicate(X,&Y);
462: /* MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); */
463: /* MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); */
464: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&D);
465: VecZeroEntries(X);
466: VecSetValue(X,col,1.0,INSERT_VALUES);
467: VecAssemblyBegin(X);
468: VecAssemblyEnd(X);
469: MatMult(mult_dup ? D : B,X,Y);
470: MatView(view_dup ? D : B,PETSC_VIEWER_STDOUT_WORLD);
471: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
472: VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
473: MatDestroy(&D);
474: VecDestroy(&Y);
475: }
477: DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc);
478: PetscFree(user);
480: ISDestroy(&isg[0]);
481: ISDestroy(&isg[1]);
482: PetscFree(isg);
483: VecDestroy(&X);
484: VecDestroy(&F);
485: MatDestroy(&B);
486: DMDestroy(&dau);
487: DMDestroy(&dak);
488: DMDestroy(&pack);
489: SNESDestroy(&snes);
490: PetscFinalize();
491: return 0;
492: }