Actual source code: ex28.c
petsc-3.10.5 2019-03-28
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: /*T
27: TODO: Need to determine if deprecated
28: T*/
30: #include <petscsnes.h>
31: #include <petscdm.h>
32: #include <petscdmda.h>
33: #include <petscdmcomposite.h>
35: typedef struct _UserCtx *User;
36: struct _UserCtx {
37: PetscInt ptype;
38: DM pack;
39: Vec Uloc,Kloc;
40: };
42: static PetscErrorCode FormFunctionLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])
43: {
44: PetscReal hx = 1./info->mx;
45: PetscInt i;
48: for (i=info->xs; i<info->xs+info->xm; i++) {
49: if (i == 0) f[i] = 1./hx*u[i];
50: else if (i == info->mx-1) f[i] = 1./hx*(u[i] - 1.0);
51: else f[i] = hx*((k[i-1]*(u[i]-u[i-1]) - k[i]*(u[i+1]-u[i]))/(hx*hx) - 1.0);
52: }
53: return(0);
54: }
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: }
73: static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx)
74: {
75: User user = (User)ctx;
76: DM dau,dak;
77: DMDALocalInfo infou,infok;
78: PetscScalar *u,*k;
79: PetscScalar *fu,*fk;
81: Vec Uloc,Kloc,Fu,Fk;
84: DMCompositeGetEntries(user->pack,&dau,&dak);
85: DMDAGetLocalInfo(dau,&infou);
86: DMDAGetLocalInfo(dak,&infok);
87: DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
88: switch (user->ptype) {
89: case 0:
90: DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
91: DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc);
92: DMDAVecGetArray(dau,Uloc,&u);
93: DMDAVecGetArray(dak,user->Kloc,&k);
94: DMDAVecGetArray(dau,F,&fu);
95: FormFunctionLocal_U(user,&infou,u,k,fu);
96: DMDAVecRestoreArray(dau,F,&fu);
97: DMDAVecRestoreArray(dau,Uloc,&u);
98: DMDAVecRestoreArray(dak,user->Kloc,&k);
99: break;
100: case 1:
101: DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
102: DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc);
103: DMDAVecGetArray(dau,user->Uloc,&u);
104: DMDAVecGetArray(dak,Kloc,&k);
105: DMDAVecGetArray(dak,F,&fk);
106: FormFunctionLocal_K(user,&infok,u,k,fk);
107: DMDAVecRestoreArray(dak,F,&fk);
108: DMDAVecRestoreArray(dau,user->Uloc,&u);
109: DMDAVecRestoreArray(dak,Kloc,&k);
110: break;
111: case 2:
112: DMCompositeScatter(user->pack,X,Uloc,Kloc);
113: DMDAVecGetArray(dau,Uloc,&u);
114: DMDAVecGetArray(dak,Kloc,&k);
115: DMCompositeGetAccess(user->pack,F,&Fu,&Fk);
116: DMDAVecGetArray(dau,Fu,&fu);
117: DMDAVecGetArray(dak,Fk,&fk);
118: FormFunctionLocal_U(user,&infou,u,k,fu);
119: FormFunctionLocal_K(user,&infok,u,k,fk);
120: DMDAVecRestoreArray(dau,Fu,&fu);
121: DMDAVecRestoreArray(dak,Fk,&fk);
122: DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk);
123: DMDAVecRestoreArray(dau,Uloc,&u);
124: DMDAVecRestoreArray(dak,Kloc,&k);
125: break;
126: }
127: DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
128: return(0);
129: }
131: static PetscErrorCode FormJacobianLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Buu)
132: {
133: PetscReal hx = 1./info->mx;
135: PetscInt i;
138: for (i=info->xs; i<info->xs+info->xm; i++) {
139: PetscInt row = i-info->gxs,cols[] = {row-1,row,row+1};
140: PetscScalar val = 1./hx;
141: if (i == 0) {
142: MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);
143: } else if (i == info->mx-1) {
144: MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES);
145: } else {
146: PetscScalar vals[] = {-k[i-1]/hx,(k[i-1]+k[i])/hx,-k[i]/hx};
147: MatSetValuesLocal(Buu,1,&row,3,cols,vals,INSERT_VALUES);
148: }
149: }
150: return(0);
151: }
153: static PetscErrorCode FormJacobianLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Bkk)
154: {
155: PetscReal hx = 1./info->mx;
157: PetscInt i;
160: for (i=info->xs; i<info->xs+info->xm; i++) {
161: PetscInt row = i-info->gxs;
162: PetscScalar vals[] = {hx*(PetscExpScalar(k[i]-1.)+1.)};
163: MatSetValuesLocal(Bkk,1,&row,1,&row,vals,INSERT_VALUES);
164: }
165: return(0);
166: }
168: static PetscErrorCode FormJacobianLocal_UK(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Buk)
169: {
170: PetscReal hx = 1./info->mx;
172: PetscInt i;
173: PetscInt row,cols[2];
174: PetscScalar vals[2];
177: if (!Buk) return(0); /* Not assembling this block */
178: for (i=info->xs; i<info->xs+info->xm; i++) {
179: if (i == 0 || i == info->mx-1) continue;
180: row = i-info->gxs;
181: cols[0] = i-1-infok->gxs; vals[0] = (u[i]-u[i-1])/hx;
182: cols[1] = i-infok->gxs; vals[1] = (u[i]-u[i+1])/hx;
183: MatSetValuesLocal(Buk,1,&row,2,cols,vals,INSERT_VALUES);
184: }
185: return(0);
186: }
188: static PetscErrorCode FormJacobianLocal_KU(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Bku)
189: {
191: PetscInt i;
192: PetscReal hx = 1./(info->mx-1);
195: if (!Bku) return(0); /* Not assembling this block */
196: for (i=infok->xs; i<infok->xs+infok->xm; i++) {
197: PetscInt row = i-infok->gxs,cols[2];
198: PetscScalar vals[2];
199: const PetscScalar
200: ubar = 0.5*(u[i]+u[i+1]),
201: ubar_L = 0.5,
202: ubar_R = 0.5,
203: gradu = (u[i+1]-u[i])/hx,
204: gradu_L = -1./hx,
205: gradu_R = 1./hx,
206: g = 1. + PetscSqr(gradu),
207: g_gradu = 2.*gradu,
208: w = 1./(1.+ubar) + 1./g,
209: w_ubar = -1./PetscSqr(1.+ubar),
210: w_gradu = -g_gradu/PetscSqr(g),
211: iw = 1./w,
212: iw_ubar = -w_ubar * PetscSqr(iw),
213: iw_gradu = -w_gradu * PetscSqr(iw);
214: cols[0] = i-info->gxs; vals[0] = -hx*(iw_ubar*ubar_L + iw_gradu*gradu_L);
215: cols[1] = i+1-info->gxs; vals[1] = -hx*(iw_ubar*ubar_R + iw_gradu*gradu_R);
216: MatSetValuesLocal(Bku,1,&row,2,cols,vals,INSERT_VALUES);
217: }
218: return(0);
219: }
221: static PetscErrorCode FormJacobian_All(SNES snes,Vec X,Mat J,Mat B,void *ctx)
222: {
223: User user = (User)ctx;
224: DM dau,dak;
225: DMDALocalInfo infou,infok;
226: PetscScalar *u,*k;
228: Vec Uloc,Kloc;
231: DMCompositeGetEntries(user->pack,&dau,&dak);
232: DMDAGetLocalInfo(dau,&infou);
233: DMDAGetLocalInfo(dak,&infok);
234: DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc);
235: switch (user->ptype) {
236: case 0:
237: DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc);
238: DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc);
239: DMDAVecGetArray(dau,Uloc,&u);
240: DMDAVecGetArray(dak,user->Kloc,&k);
241: FormJacobianLocal_U(user,&infou,u,k,B);
242: DMDAVecRestoreArray(dau,Uloc,&u);
243: DMDAVecRestoreArray(dak,user->Kloc,&k);
244: break;
245: case 1:
246: DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc);
247: DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc);
248: DMDAVecGetArray(dau,user->Uloc,&u);
249: DMDAVecGetArray(dak,Kloc,&k);
250: FormJacobianLocal_K(user,&infok,u,k,B);
251: DMDAVecRestoreArray(dau,user->Uloc,&u);
252: DMDAVecRestoreArray(dak,Kloc,&k);
253: break;
254: case 2: {
255: Mat Buu,Buk,Bku,Bkk;
256: PetscBool nest;
257: IS *is;
258: DMCompositeScatter(user->pack,X,Uloc,Kloc);
259: DMDAVecGetArray(dau,Uloc,&u);
260: DMDAVecGetArray(dak,Kloc,&k);
261: DMCompositeGetLocalISs(user->pack,&is);
262: MatGetLocalSubMatrix(B,is[0],is[0],&Buu);
263: MatGetLocalSubMatrix(B,is[0],is[1],&Buk);
264: MatGetLocalSubMatrix(B,is[1],is[0],&Bku);
265: MatGetLocalSubMatrix(B,is[1],is[1],&Bkk);
266: FormJacobianLocal_U(user,&infou,u,k,Buu);
267: PetscObjectTypeCompare((PetscObject)B,MATNEST,&nest);
268: if (!nest) {
269: /*
270: DMCreateMatrix_Composite() for a nested matrix does not generate off-block matrices that one can call MatSetValuesLocal() on, it just creates dummy
271: matrices with no entries; there cannot insert values into them. Commit b6480e041dd2293a65f96222772d68cdb4ed6306
272: changed Mat_Nest() from returning NULL pointers for these submatrices to dummy matrices because PCFIELDSPLIT could not
273: handle the returned null matrices.
274: */
275: FormJacobianLocal_UK(user,&infou,&infok,u,k,Buk);
276: FormJacobianLocal_KU(user,&infou,&infok,u,k,Bku);
277: }
278: FormJacobianLocal_K(user,&infok,u,k,Bkk);
279: MatRestoreLocalSubMatrix(B,is[0],is[0],&Buu);
280: MatRestoreLocalSubMatrix(B,is[0],is[1],&Buk);
281: MatRestoreLocalSubMatrix(B,is[1],is[0],&Bku);
282: MatRestoreLocalSubMatrix(B,is[1],is[1],&Bkk);
283: DMDAVecRestoreArray(dau,Uloc,&u);
284: DMDAVecRestoreArray(dak,Kloc,&k);
286: ISDestroy(&is[0]);
287: ISDestroy(&is[1]);
288: PetscFree(is);
289: } break;
290: }
291: DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc);
292: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
293: MatAssemblyEnd (B,MAT_FINAL_ASSEMBLY);
294: if (J != B) {
295: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
296: MatAssemblyEnd (J,MAT_FINAL_ASSEMBLY);
297: }
298: return(0);
299: }
301: static PetscErrorCode FormInitial_Coupled(User user,Vec X)
302: {
304: DM dau,dak;
305: DMDALocalInfo infou,infok;
306: Vec Xu,Xk;
307: PetscScalar *u,*k,hx;
308: PetscInt i;
311: DMCompositeGetEntries(user->pack,&dau,&dak);
312: DMCompositeGetAccess(user->pack,X,&Xu,&Xk);
313: DMDAVecGetArray(dau,Xu,&u);
314: DMDAVecGetArray(dak,Xk,&k);
315: DMDAGetLocalInfo(dau,&infou);
316: DMDAGetLocalInfo(dak,&infok);
317: hx = 1./(infok.mx);
318: for (i=infou.xs; i<infou.xs+infou.xm; i++) u[i] = (PetscScalar)i*hx * (1.-(PetscScalar)i*hx);
319: for (i=infok.xs; i<infok.xs+infok.xm; i++) k[i] = 1.0 + 0.5*PetscSinScalar(2*PETSC_PI*i*hx);
320: DMDAVecRestoreArray(dau,Xu,&u);
321: DMDAVecRestoreArray(dak,Xk,&k);
322: DMCompositeRestoreAccess(user->pack,X,&Xu,&Xk);
323: DMCompositeScatter(user->pack,X,user->Uloc,user->Kloc);
324: return(0);
325: }
327: int main(int argc, char *argv[])
328: {
330: DM dau,dak,pack;
331: const PetscInt *lxu;
332: PetscInt *lxk,m,sizes;
333: User user;
334: SNES snes;
335: Vec X,F,Xu,Xk,Fu,Fk;
336: Mat B;
337: IS *isg;
338: PetscBool pass_dm;
340: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
341: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,10,1,1,NULL,&dau);
342: DMSetOptionsPrefix(dau,"u_");
343: DMSetFromOptions(dau);
344: DMSetUp(dau);
345: DMDAGetOwnershipRanges(dau,&lxu,0,0);
346: DMDAGetInfo(dau,0, &m,0,0, &sizes,0,0, 0,0,0,0,0,0);
347: PetscMalloc1(sizes,&lxk);
348: PetscMemcpy(lxk,lxu,sizes*sizeof(*lxk));
349: lxk[0]--;
350: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m-1,1,1,lxk,&dak);
351: DMSetOptionsPrefix(dak,"k_");
352: DMSetFromOptions(dak);
353: DMSetUp(dak);
354: PetscFree(lxk);
356: DMCompositeCreate(PETSC_COMM_WORLD,&pack);
357: DMSetOptionsPrefix(pack,"pack_");
358: DMCompositeAddDM(pack,dau);
359: DMCompositeAddDM(pack,dak);
360: DMDASetFieldName(dau,0,"u");
361: DMDASetFieldName(dak,0,"k");
362: DMSetFromOptions(pack);
364: DMCreateGlobalVector(pack,&X);
365: VecDuplicate(X,&F);
367: PetscNew(&user);
369: 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,NULL,"Coupled problem options","SNES");
376: {
377: user->ptype = 0; pass_dm = PETSC_TRUE;
379: PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,NULL);
380: PetscOptionsBool("-pass_dm","Pass the packed DM to SNES to use when determining splits and forward into splits",0,pass_dm,&pass_dm,NULL);
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,&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,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,&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,NULL,Xk);
409: DMCompositeRestoreAccess(pack,X,0,&Xk);
410: DMCompositeRestoreAccess(pack,F,0,&Fk);
411: break;
412: case 2:
413: DMCreateMatrix(pack,&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,NULL,X);
433: break;
434: }
435: VecViewFromOptions(X,NULL,"-view_sol");
437: if (0) {
438: PetscInt col = 0;
439: PetscBool mult_dup = PETSC_FALSE,view_dup = PETSC_FALSE;
440: Mat D;
441: Vec Y;
443: PetscOptionsGetInt(NULL,0,"-col",&col,0);
444: PetscOptionsGetBool(NULL,0,"-mult_dup",&mult_dup,0);
445: PetscOptionsGetBool(NULL,0,"-view_dup",&view_dup,0);
447: VecDuplicate(X,&Y);
448: /* MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); */
449: /* MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); */
450: MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&D);
451: VecZeroEntries(X);
452: VecSetValue(X,col,1.0,INSERT_VALUES);
453: VecAssemblyBegin(X);
454: VecAssemblyEnd(X);
455: MatMult(mult_dup ? D : B,X,Y);
456: MatView(view_dup ? D : B,PETSC_VIEWER_STDOUT_WORLD);
457: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
458: VecView(Y,PETSC_VIEWER_STDOUT_WORLD);
459: MatDestroy(&D);
460: VecDestroy(&Y);
461: }
463: DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc);
464: PetscFree(user);
466: ISDestroy(&isg[0]);
467: ISDestroy(&isg[1]);
468: PetscFree(isg);
469: VecDestroy(&X);
470: VecDestroy(&F);
471: MatDestroy(&B);
472: DMDestroy(&dau);
473: DMDestroy(&dak);
474: DMDestroy(&pack);
475: SNESDestroy(&snes);
476: PetscFinalize();
477: return ierr;
478: }
480: /*TEST
482: build:
483: requires: c99
485: test:
486: suffix: 0
487: nsize: 3
488: args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 0
490: test:
491: suffix: 1
492: nsize: 3
493: args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 1
495: test:
496: suffix: 2
497: nsize: 3
498: args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 2
500: test:
501: suffix: 3
502: nsize: 3
503: args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -ksp_monitor_short -problem_type 2 -snes_mf_operator -pack_dm_mat_type {{aij nest}} -pc_type fieldsplit -pc_fieldsplit_dm_splits -pc_fieldsplit_type additive -fieldsplit_u_ksp_type gmres -fieldsplit_k_pc_type jacobi
505: test:
506: suffix: 4
507: nsize: 6
508: args: -u_da_grid_x 257 -snes_converged_reason -snes_monitor_short -ksp_monitor_short -problem_type 2 -snes_mf_operator -pack_dm_mat_type aij -pc_type fieldsplit -pc_fieldsplit_type multiplicative -fieldsplit_u_ksp_type gmres -fieldsplit_u_ksp_pc_side right -fieldsplit_u_pc_type mg -fieldsplit_u_pc_mg_levels 4 -fieldsplit_u_mg_levels_ksp_type richardson -fieldsplit_u_mg_levels_ksp_max_it 1 -fieldsplit_u_mg_levels_pc_type sor -fieldsplit_u_pc_mg_galerkin pmat -fieldsplit_u_ksp_converged_reason -fieldsplit_k_pc_type jacobi
509: requires: !single
511: test:
512: suffix: glvis_composite_da_1d
513: args: -u_da_grid_x 20 -problem_type 0 -snes_monitor_solution glvis:
515: TEST*/