Actual source code: ex28.c

  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: };

 38: static PetscErrorCode FormFunctionLocal_U(User user, DMDALocalInfo *info, const PetscScalar u[], const PetscScalar k[], PetscScalar f[])
 39: {
 40:   PetscReal hx = 1. / info->mx;
 41:   PetscInt  i;

 43:   PetscFunctionBeginUser;
 44:   for (i = info->xs; i < info->xs + info->xm; i++) {
 45:     if (i == 0) f[i] = 1. / hx * u[i];
 46:     else if (i == info->mx - 1) f[i] = 1. / hx * (u[i] - 1.0);
 47:     else f[i] = hx * ((k[i - 1] * (u[i] - u[i - 1]) - k[i] * (u[i + 1] - u[i])) / (hx * hx) - 1.0);
 48:   }
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

 52: static PetscErrorCode FormFunctionLocal_K(User user, DMDALocalInfo *info, const PetscScalar u[], const PetscScalar k[], PetscScalar f[])
 53: {
 54:   PetscReal hx = 1. / info->mx;
 55:   PetscInt  i;

 57:   PetscFunctionBeginUser;
 58:   for (i = info->xs; i < info->xs + info->xm; i++) {
 59:     const PetscScalar ubar = 0.5 * (u[i + 1] + u[i]), gradu = (u[i + 1] - u[i]) / hx, g = 1. + gradu * gradu, w = 1. / (1. + ubar) + 1. / g;
 60:     f[i] = hx * (PetscExpScalar(k[i] - 1.0) + k[i] - 1. / w);
 61:   }
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

 65: static PetscErrorCode FormFunction_All(SNES snes, Vec X, Vec F, void *ctx)
 66: {
 67:   User          user = (User)ctx;
 68:   DM            dau, dak;
 69:   DMDALocalInfo infou, infok;
 70:   PetscScalar  *u, *k;
 71:   PetscScalar  *fu, *fk;
 72:   Vec           Uloc, Kloc, Fu, Fk;

 74:   PetscFunctionBeginUser;
 75:   PetscCall(DMCompositeGetEntries(user->pack, &dau, &dak));
 76:   PetscCall(DMDAGetLocalInfo(dau, &infou));
 77:   PetscCall(DMDAGetLocalInfo(dak, &infok));
 78:   PetscCall(DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc));
 79:   switch (user->ptype) {
 80:   case 0:
 81:     PetscCall(DMGlobalToLocalBegin(dau, X, INSERT_VALUES, Uloc));
 82:     PetscCall(DMGlobalToLocalEnd(dau, X, INSERT_VALUES, Uloc));
 83:     PetscCall(DMDAVecGetArray(dau, Uloc, &u));
 84:     PetscCall(DMDAVecGetArray(dak, user->Kloc, &k));
 85:     PetscCall(DMDAVecGetArray(dau, F, &fu));
 86:     PetscCall(FormFunctionLocal_U(user, &infou, u, k, fu));
 87:     PetscCall(DMDAVecRestoreArray(dau, F, &fu));
 88:     PetscCall(DMDAVecRestoreArray(dau, Uloc, &u));
 89:     PetscCall(DMDAVecRestoreArray(dak, user->Kloc, &k));
 90:     break;
 91:   case 1:
 92:     PetscCall(DMGlobalToLocalBegin(dak, X, INSERT_VALUES, Kloc));
 93:     PetscCall(DMGlobalToLocalEnd(dak, X, INSERT_VALUES, Kloc));
 94:     PetscCall(DMDAVecGetArray(dau, user->Uloc, &u));
 95:     PetscCall(DMDAVecGetArray(dak, Kloc, &k));
 96:     PetscCall(DMDAVecGetArray(dak, F, &fk));
 97:     PetscCall(FormFunctionLocal_K(user, &infok, u, k, fk));
 98:     PetscCall(DMDAVecRestoreArray(dak, F, &fk));
 99:     PetscCall(DMDAVecRestoreArray(dau, user->Uloc, &u));
100:     PetscCall(DMDAVecRestoreArray(dak, Kloc, &k));
101:     break;
102:   case 2:
103:     PetscCall(DMCompositeScatter(user->pack, X, Uloc, Kloc));
104:     PetscCall(DMDAVecGetArray(dau, Uloc, &u));
105:     PetscCall(DMDAVecGetArray(dak, Kloc, &k));
106:     PetscCall(DMCompositeGetAccess(user->pack, F, &Fu, &Fk));
107:     PetscCall(DMDAVecGetArray(dau, Fu, &fu));
108:     PetscCall(DMDAVecGetArray(dak, Fk, &fk));
109:     PetscCall(FormFunctionLocal_U(user, &infou, u, k, fu));
110:     PetscCall(FormFunctionLocal_K(user, &infok, u, k, fk));
111:     PetscCall(DMDAVecRestoreArray(dau, Fu, &fu));
112:     PetscCall(DMDAVecRestoreArray(dak, Fk, &fk));
113:     PetscCall(DMCompositeRestoreAccess(user->pack, F, &Fu, &Fk));
114:     PetscCall(DMDAVecRestoreArray(dau, Uloc, &u));
115:     PetscCall(DMDAVecRestoreArray(dak, Kloc, &k));
116:     break;
117:   }
118:   PetscCall(DMCompositeRestoreLocalVectors(user->pack, &Uloc, &Kloc));
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

122: static PetscErrorCode FormJacobianLocal_U(User user, DMDALocalInfo *info, const PetscScalar u[], const PetscScalar k[], Mat Buu)
123: {
124:   PetscReal hx = 1. / info->mx;
125:   PetscInt  i;

127:   PetscFunctionBeginUser;
128:   for (i = info->xs; i < info->xs + info->xm; i++) {
129:     PetscInt    row = i - info->gxs, cols[] = {row - 1, row, row + 1};
130:     PetscScalar val = 1. / hx;
131:     if (i == 0) {
132:       PetscCall(MatSetValuesLocal(Buu, 1, &row, 1, &row, &val, INSERT_VALUES));
133:     } else if (i == info->mx - 1) {
134:       PetscCall(MatSetValuesLocal(Buu, 1, &row, 1, &row, &val, INSERT_VALUES));
135:     } else {
136:       PetscScalar vals[] = {-k[i - 1] / hx, (k[i - 1] + k[i]) / hx, -k[i] / hx};
137:       PetscCall(MatSetValuesLocal(Buu, 1, &row, 3, cols, vals, INSERT_VALUES));
138:     }
139:   }
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: static PetscErrorCode FormJacobianLocal_K(User user, DMDALocalInfo *info, const PetscScalar u[], const PetscScalar k[], Mat Bkk)
144: {
145:   PetscReal hx = 1. / info->mx;
146:   PetscInt  i;

148:   PetscFunctionBeginUser;
149:   for (i = info->xs; i < info->xs + info->xm; i++) {
150:     PetscInt    row    = i - info->gxs;
151:     PetscScalar vals[] = {hx * (PetscExpScalar(k[i] - 1.) + (PetscScalar)1.)};
152:     PetscCall(MatSetValuesLocal(Bkk, 1, &row, 1, &row, vals, INSERT_VALUES));
153:   }
154:   PetscFunctionReturn(PETSC_SUCCESS);
155: }

157: static PetscErrorCode FormJacobianLocal_UK(User user, DMDALocalInfo *info, DMDALocalInfo *infok, const PetscScalar u[], const PetscScalar k[], Mat Buk)
158: {
159:   PetscReal   hx = 1. / info->mx;
160:   PetscInt    i;
161:   PetscInt    row, cols[2];
162:   PetscScalar vals[2];

164:   PetscFunctionBeginUser;
165:   if (!Buk) PetscFunctionReturn(PETSC_SUCCESS); /* Not assembling this block */
166:   for (i = info->xs; i < info->xs + info->xm; i++) {
167:     if (i == 0 || i == info->mx - 1) continue;
168:     row     = i - info->gxs;
169:     cols[0] = i - 1 - infok->gxs;
170:     vals[0] = (u[i] - u[i - 1]) / hx;
171:     cols[1] = i - infok->gxs;
172:     vals[1] = (u[i] - u[i + 1]) / hx;
173:     PetscCall(MatSetValuesLocal(Buk, 1, &row, 2, cols, vals, INSERT_VALUES));
174:   }
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: static PetscErrorCode FormJacobianLocal_KU(User user, DMDALocalInfo *info, DMDALocalInfo *infok, const PetscScalar u[], const PetscScalar k[], Mat Bku)
179: {
180:   PetscInt  i;
181:   PetscReal hx = 1. / (info->mx - 1);

183:   PetscFunctionBeginUser;
184:   if (!Bku) PetscFunctionReturn(PETSC_SUCCESS); /* Not assembling this block */
185:   for (i = infok->xs; i < infok->xs + infok->xm; i++) {
186:     PetscInt    row = i - infok->gxs, cols[2];
187:     PetscScalar vals[2];
188:     const PetscScalar ubar = 0.5 * (u[i] + u[i + 1]), ubar_L = 0.5, ubar_R = 0.5, gradu = (u[i + 1] - u[i]) / hx, gradu_L = -1. / hx, gradu_R = 1. / hx, g = 1. + PetscSqr(gradu), g_gradu = 2. * gradu, w = 1. / (1. + ubar) + 1. / g, w_ubar = -1. / PetscSqr(1. + ubar), w_gradu = -g_gradu / PetscSqr(g), iw = 1. / w, iw_ubar = -w_ubar * PetscSqr(iw), iw_gradu = -w_gradu * PetscSqr(iw);
189:     cols[0] = i - info->gxs;
190:     vals[0] = -hx * (iw_ubar * ubar_L + iw_gradu * gradu_L);
191:     cols[1] = i + 1 - info->gxs;
192:     vals[1] = -hx * (iw_ubar * ubar_R + iw_gradu * gradu_R);
193:     PetscCall(MatSetValuesLocal(Bku, 1, &row, 2, cols, vals, INSERT_VALUES));
194:   }
195:   PetscFunctionReturn(PETSC_SUCCESS);
196: }

198: static PetscErrorCode FormJacobian_All(SNES snes, Vec X, Mat J, Mat B, void *ctx)
199: {
200:   User          user = (User)ctx;
201:   DM            dau, dak;
202:   DMDALocalInfo infou, infok;
203:   PetscScalar  *u, *k;
204:   Vec           Uloc, Kloc;

206:   PetscFunctionBeginUser;
207:   PetscCall(DMCompositeGetEntries(user->pack, &dau, &dak));
208:   PetscCall(DMDAGetLocalInfo(dau, &infou));
209:   PetscCall(DMDAGetLocalInfo(dak, &infok));
210:   PetscCall(DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc));
211:   switch (user->ptype) {
212:   case 0:
213:     PetscCall(DMGlobalToLocalBegin(dau, X, INSERT_VALUES, Uloc));
214:     PetscCall(DMGlobalToLocalEnd(dau, X, INSERT_VALUES, Uloc));
215:     PetscCall(DMDAVecGetArray(dau, Uloc, &u));
216:     PetscCall(DMDAVecGetArray(dak, user->Kloc, &k));
217:     PetscCall(FormJacobianLocal_U(user, &infou, u, k, B));
218:     PetscCall(DMDAVecRestoreArray(dau, Uloc, &u));
219:     PetscCall(DMDAVecRestoreArray(dak, user->Kloc, &k));
220:     break;
221:   case 1:
222:     PetscCall(DMGlobalToLocalBegin(dak, X, INSERT_VALUES, Kloc));
223:     PetscCall(DMGlobalToLocalEnd(dak, X, INSERT_VALUES, Kloc));
224:     PetscCall(DMDAVecGetArray(dau, user->Uloc, &u));
225:     PetscCall(DMDAVecGetArray(dak, Kloc, &k));
226:     PetscCall(FormJacobianLocal_K(user, &infok, u, k, B));
227:     PetscCall(DMDAVecRestoreArray(dau, user->Uloc, &u));
228:     PetscCall(DMDAVecRestoreArray(dak, Kloc, &k));
229:     break;
230:   case 2: {
231:     Mat       Buu, Buk, Bku, Bkk;
232:     PetscBool nest;
233:     IS       *is;
234:     PetscCall(DMCompositeScatter(user->pack, X, Uloc, Kloc));
235:     PetscCall(DMDAVecGetArray(dau, Uloc, &u));
236:     PetscCall(DMDAVecGetArray(dak, Kloc, &k));
237:     PetscCall(DMCompositeGetLocalISs(user->pack, &is));
238:     PetscCall(MatGetLocalSubMatrix(B, is[0], is[0], &Buu));
239:     PetscCall(MatGetLocalSubMatrix(B, is[0], is[1], &Buk));
240:     PetscCall(MatGetLocalSubMatrix(B, is[1], is[0], &Bku));
241:     PetscCall(MatGetLocalSubMatrix(B, is[1], is[1], &Bkk));
242:     PetscCall(FormJacobianLocal_U(user, &infou, u, k, Buu));
243:     PetscCall(PetscObjectTypeCompare((PetscObject)B, MATNEST, &nest));
244:     if (!nest) {
245:       /*
246:          DMCreateMatrix_Composite()  for a nested matrix does not generate off-block matrices that one can call MatSetValuesLocal() on, it just creates dummy
247:          matrices with no entries; there cannot insert values into them. Commit b6480e041dd2293a65f96222772d68cdb4ed6306
248:          changed Mat_Nest() from returning NULL pointers for these submatrices to dummy matrices because PCFIELDSPLIT could not
249:          handle the returned null matrices.
250:       */
251:       PetscCall(FormJacobianLocal_UK(user, &infou, &infok, u, k, Buk));
252:       PetscCall(FormJacobianLocal_KU(user, &infou, &infok, u, k, Bku));
253:     }
254:     PetscCall(FormJacobianLocal_K(user, &infok, u, k, Bkk));
255:     PetscCall(MatRestoreLocalSubMatrix(B, is[0], is[0], &Buu));
256:     PetscCall(MatRestoreLocalSubMatrix(B, is[0], is[1], &Buk));
257:     PetscCall(MatRestoreLocalSubMatrix(B, is[1], is[0], &Bku));
258:     PetscCall(MatRestoreLocalSubMatrix(B, is[1], is[1], &Bkk));
259:     PetscCall(DMDAVecRestoreArray(dau, Uloc, &u));
260:     PetscCall(DMDAVecRestoreArray(dak, Kloc, &k));

262:     PetscCall(ISDestroy(&is[0]));
263:     PetscCall(ISDestroy(&is[1]));
264:     PetscCall(PetscFree(is));
265:   } break;
266:   }
267:   PetscCall(DMCompositeRestoreLocalVectors(user->pack, &Uloc, &Kloc));
268:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
269:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
270:   if (J != B) {
271:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
272:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
273:   }
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: static PetscErrorCode FormInitial_Coupled(User user, Vec X)
278: {
279:   DM            dau, dak;
280:   DMDALocalInfo infou, infok;
281:   Vec           Xu, Xk;
282:   PetscScalar  *u, *k, hx;
283:   PetscInt      i;

285:   PetscFunctionBeginUser;
286:   PetscCall(DMCompositeGetEntries(user->pack, &dau, &dak));
287:   PetscCall(DMCompositeGetAccess(user->pack, X, &Xu, &Xk));
288:   PetscCall(DMDAVecGetArray(dau, Xu, &u));
289:   PetscCall(DMDAVecGetArray(dak, Xk, &k));
290:   PetscCall(DMDAGetLocalInfo(dau, &infou));
291:   PetscCall(DMDAGetLocalInfo(dak, &infok));
292:   hx = 1. / (infok.mx);
293:   for (i = infou.xs; i < infou.xs + infou.xm; i++) u[i] = (PetscScalar)i * hx * (1. - (PetscScalar)i * hx);
294:   for (i = infok.xs; i < infok.xs + infok.xm; i++) k[i] = 1.0 + 0.5 * PetscSinScalar(2 * PETSC_PI * i * hx);
295:   PetscCall(DMDAVecRestoreArray(dau, Xu, &u));
296:   PetscCall(DMDAVecRestoreArray(dak, Xk, &k));
297:   PetscCall(DMCompositeRestoreAccess(user->pack, X, &Xu, &Xk));
298:   PetscCall(DMCompositeScatter(user->pack, X, user->Uloc, user->Kloc));
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: int main(int argc, char *argv[])
303: {
304:   DM              dau, dak, pack;
305:   const PetscInt *lxu;
306:   PetscInt       *lxk, m, sizes;
307:   User            user;
308:   SNES            snes;
309:   Vec             X, F, Xu, Xk, Fu, Fk;
310:   Mat             B;
311:   IS             *isg;
312:   PetscBool       pass_dm;

314:   PetscFunctionBeginUser;
315:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
316:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 10, 1, 1, NULL, &dau));
317:   PetscCall(DMSetOptionsPrefix(dau, "u_"));
318:   PetscCall(DMSetFromOptions(dau));
319:   PetscCall(DMSetUp(dau));
320:   PetscCall(DMDAGetOwnershipRanges(dau, &lxu, 0, 0));
321:   PetscCall(DMDAGetInfo(dau, 0, &m, 0, 0, &sizes, 0, 0, 0, 0, 0, 0, 0, 0));
322:   PetscCall(PetscMalloc1(sizes, &lxk));
323:   PetscCall(PetscArraycpy(lxk, lxu, sizes));
324:   lxk[0]--;
325:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m - 1, 1, 1, lxk, &dak));
326:   PetscCall(DMSetOptionsPrefix(dak, "k_"));
327:   PetscCall(DMSetFromOptions(dak));
328:   PetscCall(DMSetUp(dak));
329:   PetscCall(PetscFree(lxk));

331:   PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &pack));
332:   PetscCall(DMSetOptionsPrefix(pack, "pack_"));
333:   PetscCall(DMCompositeAddDM(pack, dau));
334:   PetscCall(DMCompositeAddDM(pack, dak));
335:   PetscCall(DMDASetFieldName(dau, 0, "u"));
336:   PetscCall(DMDASetFieldName(dak, 0, "k"));
337:   PetscCall(DMSetFromOptions(pack));

339:   PetscCall(DMCreateGlobalVector(pack, &X));
340:   PetscCall(VecDuplicate(X, &F));

342:   PetscCall(PetscNew(&user));

344:   user->pack = pack;

346:   PetscCall(DMCompositeGetGlobalISs(pack, &isg));
347:   PetscCall(DMCompositeGetLocalVectors(pack, &user->Uloc, &user->Kloc));
348:   PetscCall(DMCompositeScatter(pack, X, user->Uloc, user->Kloc));

350:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Coupled problem options", "SNES");
351:   {
352:     user->ptype = 0;
353:     pass_dm     = PETSC_TRUE;

355:     PetscCall(PetscOptionsInt("-problem_type", "0: solve for u only, 1: solve for k only, 2: solve for both", 0, user->ptype, &user->ptype, NULL));
356:     PetscCall(PetscOptionsBool("-pass_dm", "Pass the packed DM to SNES to use when determining splits and forward into splits", 0, pass_dm, &pass_dm, NULL));
357:   }
358:   PetscOptionsEnd();

360:   PetscCall(FormInitial_Coupled(user, X));

362:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
363:   switch (user->ptype) {
364:   case 0:
365:     PetscCall(DMCompositeGetAccess(pack, X, &Xu, PETSC_NULLPTR));
366:     PetscCall(DMCompositeGetAccess(pack, F, &Fu, PETSC_NULLPTR));
367:     PetscCall(DMCreateMatrix(dau, &B));
368:     PetscCall(SNESSetFunction(snes, Fu, FormFunction_All, user));
369:     PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
370:     PetscCall(SNESSetFromOptions(snes));
371:     PetscCall(SNESSetDM(snes, dau));
372:     PetscCall(SNESSolve(snes, NULL, Xu));
373:     PetscCall(DMCompositeRestoreAccess(pack, X, &Xu, PETSC_NULLPTR));
374:     PetscCall(DMCompositeRestoreAccess(pack, F, &Fu, PETSC_NULLPTR));
375:     break;
376:   case 1:
377:     PetscCall(DMCompositeGetAccess(pack, X, PETSC_NULLPTR, &Xk));
378:     PetscCall(DMCompositeGetAccess(pack, F, PETSC_NULLPTR, &Fk));
379:     PetscCall(DMCreateMatrix(dak, &B));
380:     PetscCall(SNESSetFunction(snes, Fk, FormFunction_All, user));
381:     PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
382:     PetscCall(SNESSetFromOptions(snes));
383:     PetscCall(SNESSetDM(snes, dak));
384:     PetscCall(SNESSolve(snes, NULL, Xk));
385:     PetscCall(DMCompositeRestoreAccess(pack, X, PETSC_NULLPTR, &Xk));
386:     PetscCall(DMCompositeRestoreAccess(pack, F, PETSC_NULLPTR, &Fk));
387:     break;
388:   case 2:
389:     PetscCall(DMCreateMatrix(pack, &B));
390:     /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */
391:     PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
392:     PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
393:     PetscCall(SNESSetFunction(snes, F, FormFunction_All, user));
394:     PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
395:     PetscCall(SNESSetFromOptions(snes));
396:     if (!pass_dm) { /* Manually provide index sets and names for the splits */
397:       KSP ksp;
398:       PC  pc;
399:       PetscCall(SNESGetKSP(snes, &ksp));
400:       PetscCall(KSPGetPC(ksp, &pc));
401:       PetscCall(PCFieldSplitSetIS(pc, "u", isg[0]));
402:       PetscCall(PCFieldSplitSetIS(pc, "k", isg[1]));
403:     } else {
404:       /* The same names come from the options prefix for dau and dak. This option can support geometric multigrid inside
405:        * of splits, but it requires using a DM (perhaps your own implementation). */
406:       PetscCall(SNESSetDM(snes, pack));
407:     }
408:     PetscCall(SNESSolve(snes, NULL, X));
409:     break;
410:   }
411:   PetscCall(VecViewFromOptions(X, NULL, "-view_sol"));

413:   if (0) {
414:     PetscInt  col      = 0;
415:     PetscBool mult_dup = PETSC_FALSE, view_dup = PETSC_FALSE;
416:     Mat       D;
417:     Vec       Y;

419:     PetscCall(PetscOptionsGetInt(NULL, 0, "-col", &col, 0));
420:     PetscCall(PetscOptionsGetBool(NULL, 0, "-mult_dup", &mult_dup, 0));
421:     PetscCall(PetscOptionsGetBool(NULL, 0, "-view_dup", &view_dup, 0));

423:     PetscCall(VecDuplicate(X, &Y));
424:     /* PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); */
425:     /* PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); */
426:     PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &D));
427:     PetscCall(VecZeroEntries(X));
428:     PetscCall(VecSetValue(X, col, 1.0, INSERT_VALUES));
429:     PetscCall(VecAssemblyBegin(X));
430:     PetscCall(VecAssemblyEnd(X));
431:     PetscCall(MatMult(mult_dup ? D : B, X, Y));
432:     PetscCall(MatView(view_dup ? D : B, PETSC_VIEWER_STDOUT_WORLD));
433:     /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
434:     PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
435:     PetscCall(MatDestroy(&D));
436:     PetscCall(VecDestroy(&Y));
437:   }

439:   PetscCall(DMCompositeRestoreLocalVectors(pack, &user->Uloc, &user->Kloc));
440:   PetscCall(PetscFree(user));

442:   PetscCall(ISDestroy(&isg[0]));
443:   PetscCall(ISDestroy(&isg[1]));
444:   PetscCall(PetscFree(isg));
445:   PetscCall(VecDestroy(&X));
446:   PetscCall(VecDestroy(&F));
447:   PetscCall(MatDestroy(&B));
448:   PetscCall(DMDestroy(&dau));
449:   PetscCall(DMDestroy(&dak));
450:   PetscCall(DMDestroy(&pack));
451:   PetscCall(SNESDestroy(&snes));
452:   PetscCall(PetscFinalize());
453:   return 0;
454: }

456: /*TEST

458:    test:
459:       suffix: 0
460:       nsize: 3
461:       args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 0

463:    test:
464:       suffix: 1
465:       nsize: 3
466:       args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 1

468:    test:
469:       suffix: 2
470:       nsize: 3
471:       args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 2

473:    testset:
474:       suffix: 3
475:       nsize: 3
476:       output_file: output/ex28_3.out
477:       args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -ksp_monitor_short -problem_type 2 -snes_mf_operator -pc_type fieldsplit -pc_fieldsplit_dm_splits -pc_fieldsplit_type additive -fieldsplit_u_ksp_type gmres -fieldsplit_k_pc_type jacobi

479:       test:
480:         suffix: std
481:         args: -pack_dm_mat_type {{aij nest}}

483:       test:
484:         suffix: cuda
485:         requires: cuda
486:         args: -pack_dm_mat_type aijcusparse -pack_dm_vec_type cuda

488:       test:
489:         suffix: hip
490:         requires: hip
491:         args: -pack_dm_mat_type aijhipsparse -pack_dm_vec_type hip

493:       test:
494:         suffix: kok
495:         requires: kokkos_kernels
496:         args: -pack_dm_mat_type aijkokkos -pack_dm_vec_type kokkos

498:    test:
499:       requires: mumps
500:       suffix: 3_nest_lu
501:       nsize: {{1 3}}
502:       output_file: output/ex28_3.out
503:       args: -pack_dm_mat_type nest -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -ksp_monitor_short -problem_type 2 -snes_mf_operator -pc_type lu -pc_factor_mat_solver_type mumps

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:
516:       suffix: cuda
517:       nsize: 1
518:       requires: cuda
519:       args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 2 -pack_dm_mat_type aijcusparse -pack_dm_vec_type cuda

521:    test:
522:       suffix: viennacl
523:       nsize: 1
524:       requires: viennacl
525:       args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 2 -pack_dm_mat_type aijviennacl -pack_dm_vec_type viennacl

527: TEST*/