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, PetscCtx 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 row, cols[2];
161: PetscScalar vals[2];
163: PetscFunctionBeginUser;
164: if (!Buk) PetscFunctionReturn(PETSC_SUCCESS); /* Not assembling this block */
165: for (PetscInt i = info->xs; i < info->xs + info->xm; i++) {
166: if (i == 0 || i == info->mx - 1) continue;
167: row = i - info->gxs;
168: cols[0] = i - 1 - infok->gxs;
169: vals[0] = (u[i] - u[i - 1]) / hx;
170: cols[1] = i - infok->gxs;
171: vals[1] = (u[i] - u[i + 1]) / hx;
172: PetscCall(MatSetValuesLocal(Buk, 1, &row, 2, cols, vals, INSERT_VALUES));
173: }
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode FormJacobianLocal_KU(User user, DMDALocalInfo *info, DMDALocalInfo *infok, const PetscScalar u[], const PetscScalar k[], Mat Bku)
178: {
179: PetscReal hx = 1. / (info->mx - 1);
181: PetscFunctionBeginUser;
182: if (!Bku) PetscFunctionReturn(PETSC_SUCCESS); /* Not assembling this block */
183: for (PetscInt i = infok->xs; i < infok->xs + infok->xm; i++) {
184: PetscInt row = i - infok->gxs, cols[2];
185: PetscScalar vals[2];
186: 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);
187: cols[0] = i - info->gxs;
188: vals[0] = -hx * (iw_ubar * ubar_L + iw_gradu * gradu_L);
189: cols[1] = i + 1 - info->gxs;
190: vals[1] = -hx * (iw_ubar * ubar_R + iw_gradu * gradu_R);
191: PetscCall(MatSetValuesLocal(Bku, 1, &row, 2, cols, vals, INSERT_VALUES));
192: }
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: static PetscErrorCode FormJacobian_All(SNES snes, Vec X, Mat J, Mat B, PetscCtx ctx)
197: {
198: User user = (User)ctx;
199: DM dau, dak;
200: DMDALocalInfo infou, infok;
201: PetscScalar *u, *k;
202: Vec Uloc, Kloc;
204: PetscFunctionBeginUser;
205: PetscCall(DMCompositeGetEntries(user->pack, &dau, &dak));
206: PetscCall(DMDAGetLocalInfo(dau, &infou));
207: PetscCall(DMDAGetLocalInfo(dak, &infok));
208: PetscCall(DMCompositeGetLocalVectors(user->pack, &Uloc, &Kloc));
209: switch (user->ptype) {
210: case 0:
211: PetscCall(DMGlobalToLocalBegin(dau, X, INSERT_VALUES, Uloc));
212: PetscCall(DMGlobalToLocalEnd(dau, X, INSERT_VALUES, Uloc));
213: PetscCall(DMDAVecGetArray(dau, Uloc, &u));
214: PetscCall(DMDAVecGetArray(dak, user->Kloc, &k));
215: PetscCall(FormJacobianLocal_U(user, &infou, u, k, B));
216: PetscCall(DMDAVecRestoreArray(dau, Uloc, &u));
217: PetscCall(DMDAVecRestoreArray(dak, user->Kloc, &k));
218: break;
219: case 1:
220: PetscCall(DMGlobalToLocalBegin(dak, X, INSERT_VALUES, Kloc));
221: PetscCall(DMGlobalToLocalEnd(dak, X, INSERT_VALUES, Kloc));
222: PetscCall(DMDAVecGetArray(dau, user->Uloc, &u));
223: PetscCall(DMDAVecGetArray(dak, Kloc, &k));
224: PetscCall(FormJacobianLocal_K(user, &infok, u, k, B));
225: PetscCall(DMDAVecRestoreArray(dau, user->Uloc, &u));
226: PetscCall(DMDAVecRestoreArray(dak, Kloc, &k));
227: break;
228: case 2: {
229: Mat Buu, Buk, Bku, Bkk;
230: PetscBool nest;
231: IS *is;
232: PetscCall(DMCompositeScatter(user->pack, X, Uloc, Kloc));
233: PetscCall(DMDAVecGetArray(dau, Uloc, &u));
234: PetscCall(DMDAVecGetArray(dak, Kloc, &k));
235: PetscCall(DMCompositeGetLocalISs(user->pack, &is));
236: PetscCall(MatGetLocalSubMatrix(B, is[0], is[0], &Buu));
237: PetscCall(MatGetLocalSubMatrix(B, is[0], is[1], &Buk));
238: PetscCall(MatGetLocalSubMatrix(B, is[1], is[0], &Bku));
239: PetscCall(MatGetLocalSubMatrix(B, is[1], is[1], &Bkk));
240: PetscCall(FormJacobianLocal_U(user, &infou, u, k, Buu));
241: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATNEST, &nest));
242: if (!nest) {
243: /*
244: DMCreateMatrix_Composite() for a nested matrix does not generate off-block matrices that one can call MatSetValuesLocal() on, it just creates dummy
245: matrices with no entries; there cannot insert values into them. Commit b6480e041dd2293a65f96222772d68cdb4ed6306
246: changed Mat_Nest() from returning NULL pointers for these submatrices to dummy matrices because PCFIELDSPLIT could not
247: handle the returned null matrices.
248: */
249: PetscCall(FormJacobianLocal_UK(user, &infou, &infok, u, k, Buk));
250: PetscCall(FormJacobianLocal_KU(user, &infou, &infok, u, k, Bku));
251: }
252: PetscCall(FormJacobianLocal_K(user, &infok, u, k, Bkk));
253: PetscCall(MatRestoreLocalSubMatrix(B, is[0], is[0], &Buu));
254: PetscCall(MatRestoreLocalSubMatrix(B, is[0], is[1], &Buk));
255: PetscCall(MatRestoreLocalSubMatrix(B, is[1], is[0], &Bku));
256: PetscCall(MatRestoreLocalSubMatrix(B, is[1], is[1], &Bkk));
257: PetscCall(DMDAVecRestoreArray(dau, Uloc, &u));
258: PetscCall(DMDAVecRestoreArray(dak, Kloc, &k));
260: PetscCall(ISDestroy(&is[0]));
261: PetscCall(ISDestroy(&is[1]));
262: PetscCall(PetscFree(is));
263: } break;
264: }
265: PetscCall(DMCompositeRestoreLocalVectors(user->pack, &Uloc, &Kloc));
266: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
267: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
268: if (J != B) {
269: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
270: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: static PetscErrorCode FormInitial_Coupled(User user, Vec X)
276: {
277: DM dau, dak;
278: DMDALocalInfo infou, infok;
279: Vec Xu, Xk;
280: PetscScalar *u, *k, hx;
282: PetscFunctionBeginUser;
283: PetscCall(DMCompositeGetEntries(user->pack, &dau, &dak));
284: PetscCall(DMCompositeGetAccess(user->pack, X, &Xu, &Xk));
285: PetscCall(DMDAVecGetArray(dau, Xu, &u));
286: PetscCall(DMDAVecGetArray(dak, Xk, &k));
287: PetscCall(DMDAGetLocalInfo(dau, &infou));
288: PetscCall(DMDAGetLocalInfo(dak, &infok));
289: hx = 1. / (infok.mx);
290: for (PetscInt i = infou.xs; i < infou.xs + infou.xm; i++) u[i] = (PetscScalar)i * hx * (1. - (PetscScalar)i * hx);
291: for (PetscInt i = infok.xs; i < infok.xs + infok.xm; i++) k[i] = 1.0 + 0.5 * PetscSinScalar(2 * PETSC_PI * i * hx);
292: PetscCall(DMDAVecRestoreArray(dau, Xu, &u));
293: PetscCall(DMDAVecRestoreArray(dak, Xk, &k));
294: PetscCall(DMCompositeRestoreAccess(user->pack, X, &Xu, &Xk));
295: PetscCall(DMCompositeScatter(user->pack, X, user->Uloc, user->Kloc));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: int main(int argc, char *argv[])
300: {
301: DM dau, dak, pack;
302: const PetscInt *lxu;
303: PetscInt *lxk, m, sizes;
304: User user;
305: SNES snes;
306: Vec X, F, Xu, Xk, Fu, Fk;
307: Mat B;
308: IS *isg;
309: PetscBool pass_dm;
310: VecType vtype;
312: PetscFunctionBeginUser;
313: PetscCall(PetscInitialize(&argc, &argv, 0, help));
314: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 10, 1, 1, NULL, &dau));
315: PetscCall(DMSetOptionsPrefix(dau, "u_"));
316: PetscCall(DMSetFromOptions(dau));
317: PetscCall(DMSetUp(dau));
318: PetscCall(DMDAGetOwnershipRanges(dau, &lxu, 0, 0));
319: PetscCall(DMDAGetInfo(dau, 0, &m, 0, 0, &sizes, 0, 0, 0, 0, 0, 0, 0, 0));
320: PetscCall(PetscMalloc1(sizes, &lxk));
321: PetscCall(PetscArraycpy(lxk, lxu, sizes));
322: lxk[0]--;
323: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m - 1, 1, 1, lxk, &dak));
324: PetscCall(DMSetOptionsPrefix(dak, "k_"));
325: PetscCall(DMSetFromOptions(dak));
326: PetscCall(DMSetUp(dak));
327: PetscCall(PetscFree(lxk));
329: PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &pack));
330: PetscCall(DMSetOptionsPrefix(pack, "pack_"));
331: PetscCall(DMCompositeAddDM(pack, dau));
332: PetscCall(DMCompositeAddDM(pack, dak));
333: PetscCall(DMDASetFieldName(dau, 0, "u"));
334: PetscCall(DMDASetFieldName(dak, 0, "k"));
335: PetscCall(DMSetFromOptions(pack));
337: /* if we call -pack_dm_vec_type xxx it does not get propagated to the subDMs */
338: PetscCall(DMGetVecType(pack, &vtype));
339: PetscCall(DMSetVecType(dau, vtype));
340: PetscCall(DMSetVecType(dak, vtype));
342: PetscCall(DMCreateGlobalVector(pack, &X));
343: PetscCall(VecDuplicate(X, &F));
345: PetscCall(PetscNew(&user));
347: user->pack = pack;
349: PetscCall(DMCompositeGetGlobalISs(pack, &isg));
350: PetscCall(DMCompositeGetLocalVectors(pack, &user->Uloc, &user->Kloc));
351: PetscCall(DMCompositeScatter(pack, X, user->Uloc, user->Kloc));
353: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Coupled problem options", "SNES");
354: {
355: user->ptype = 0;
356: pass_dm = PETSC_TRUE;
358: PetscCall(PetscOptionsInt("-problem_type", "0: solve for u only, 1: solve for k only, 2: solve for both", 0, user->ptype, &user->ptype, NULL));
359: 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));
360: }
361: PetscOptionsEnd();
363: PetscCall(FormInitial_Coupled(user, X));
365: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
366: switch (user->ptype) {
367: case 0:
368: PetscCall(DMCompositeGetAccess(pack, X, &Xu, PETSC_NULLPTR));
369: PetscCall(DMCompositeGetAccess(pack, F, &Fu, PETSC_NULLPTR));
370: PetscCall(DMCreateMatrix(dau, &B));
371: PetscCall(SNESSetFunction(snes, Fu, FormFunction_All, user));
372: PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
373: PetscCall(SNESSetFromOptions(snes));
374: PetscCall(SNESSetDM(snes, dau));
375: PetscCall(SNESSolve(snes, NULL, Xu));
376: PetscCall(DMCompositeRestoreAccess(pack, X, &Xu, PETSC_NULLPTR));
377: PetscCall(DMCompositeRestoreAccess(pack, F, &Fu, PETSC_NULLPTR));
378: break;
379: case 1:
380: PetscCall(DMCompositeGetAccess(pack, X, PETSC_NULLPTR, &Xk));
381: PetscCall(DMCompositeGetAccess(pack, F, PETSC_NULLPTR, &Fk));
382: PetscCall(DMCreateMatrix(dak, &B));
383: PetscCall(SNESSetFunction(snes, Fk, FormFunction_All, user));
384: PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
385: PetscCall(SNESSetFromOptions(snes));
386: PetscCall(SNESSetDM(snes, dak));
387: PetscCall(SNESSolve(snes, NULL, Xk));
388: PetscCall(DMCompositeRestoreAccess(pack, X, PETSC_NULLPTR, &Xk));
389: PetscCall(DMCompositeRestoreAccess(pack, F, PETSC_NULLPTR, &Fk));
390: break;
391: case 2:
392: PetscCall(DMCreateMatrix(pack, &B));
393: /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */
394: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
395: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
396: PetscCall(SNESSetFunction(snes, F, FormFunction_All, user));
397: PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
398: PetscCall(SNESSetFromOptions(snes));
399: if (!pass_dm) { /* Manually provide index sets and names for the splits */
400: KSP ksp;
401: PC pc;
402: PetscCall(SNESGetKSP(snes, &ksp));
403: PetscCall(KSPGetPC(ksp, &pc));
404: PetscCall(PCFieldSplitSetIS(pc, "u", isg[0]));
405: PetscCall(PCFieldSplitSetIS(pc, "k", isg[1]));
406: } else {
407: /* The same names come from the options prefix for dau and dak. This option can support geometric multigrid inside
408: * of splits, but it requires using a DM (perhaps your own implementation). */
409: PetscCall(SNESSetDM(snes, pack));
410: }
411: PetscCall(SNESSolve(snes, NULL, X));
412: break;
413: }
414: PetscCall(VecViewFromOptions(X, NULL, "-view_sol"));
416: if (0) {
417: PetscInt col = 0;
418: PetscBool mult_dup = PETSC_FALSE, view_dup = PETSC_FALSE;
419: Mat D;
420: Vec Y;
422: PetscCall(PetscOptionsGetInt(NULL, 0, "-col", &col, 0));
423: PetscCall(PetscOptionsGetBool(NULL, 0, "-mult_dup", &mult_dup, 0));
424: PetscCall(PetscOptionsGetBool(NULL, 0, "-view_dup", &view_dup, 0));
426: PetscCall(VecDuplicate(X, &Y));
427: /* PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); */
428: /* PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); */
429: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &D));
430: PetscCall(VecZeroEntries(X));
431: PetscCall(VecSetValue(X, col, 1.0, INSERT_VALUES));
432: PetscCall(VecAssemblyBegin(X));
433: PetscCall(VecAssemblyEnd(X));
434: PetscCall(MatMult(mult_dup ? D : B, X, Y));
435: PetscCall(MatView(view_dup ? D : B, PETSC_VIEWER_STDOUT_WORLD));
436: /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
437: PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
438: PetscCall(MatDestroy(&D));
439: PetscCall(VecDestroy(&Y));
440: }
442: PetscCall(DMCompositeRestoreLocalVectors(pack, &user->Uloc, &user->Kloc));
443: PetscCall(PetscFree(user));
445: PetscCall(ISDestroy(&isg[0]));
446: PetscCall(ISDestroy(&isg[1]));
447: PetscCall(PetscFree(isg));
448: PetscCall(VecDestroy(&X));
449: PetscCall(VecDestroy(&F));
450: PetscCall(MatDestroy(&B));
451: PetscCall(DMDestroy(&dau));
452: PetscCall(DMDestroy(&dak));
453: PetscCall(DMDestroy(&pack));
454: PetscCall(SNESDestroy(&snes));
455: PetscCall(PetscFinalize());
456: return 0;
457: }
459: /*TEST
461: test:
462: suffix: 0
463: nsize: 3
464: args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 0
466: test:
467: suffix: 1
468: nsize: 3
469: args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 1
471: test:
472: suffix: 2
473: nsize: 3
474: args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 2
476: testset:
477: suffix: 3
478: nsize: 3
479: output_file: output/ex28_3.out
480: 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
482: test:
483: suffix: std
484: args: -pack_dm_mat_type {{aij nest}}
486: test:
487: suffix: cuda
488: requires: cuda
489: args: -pack_dm_mat_type aijcusparse -pack_dm_vec_type cuda
491: test:
492: suffix: hip
493: requires: hip
494: args: -pack_dm_mat_type aijhipsparse -pack_dm_vec_type hip
496: test:
497: suffix: kok
498: requires: kokkos_kernels
499: args: -pack_dm_mat_type aijkokkos -pack_dm_vec_type kokkos
501: test:
502: requires: mumps
503: suffix: 3_nest_lu
504: nsize: {{1 3}}
505: output_file: output/ex28_3.out
506: 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
508: test:
509: suffix: 4
510: nsize: 6
511: 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
512: requires: !single
514: test:
515: suffix: glvis_composite_da_1d
516: args: -u_da_grid_x 20 -problem_type 0 -snes_monitor_solution glvis:
518: test:
519: suffix: cuda
520: nsize: 1
521: requires: cuda
522: 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
524: test:
525: suffix: viennacl
526: nsize: 1
527: requires: viennacl
528: 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
530: TEST*/