1 static const char help[] = "1D multiphysics prototype with analytic Jacobians to solve individual problems and a coupled problem.\n\n";
2
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 */
25
26 #include <petscsnes.h>
27 #include <petscdm.h>
28 #include <petscdmda.h>
29 #include <petscdmcomposite.h>
30
31 typedef struct _UserCtx *User;
32 struct _UserCtx {
33 PetscInt ptype;
34 DM pack;
35 Vec Uloc, Kloc;
36 };
37
FormFunctionLocal_U(User user,DMDALocalInfo * info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])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;
42
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 }
51
FormFunctionLocal_K(User user,DMDALocalInfo * info,const PetscScalar u[],const PetscScalar k[],PetscScalar f[])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;
56
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 }
64
FormFunction_All(SNES snes,Vec X,Vec F,PetscCtx ctx)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;
73
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 }
121
FormJacobianLocal_U(User user,DMDALocalInfo * info,const PetscScalar u[],const PetscScalar k[],Mat Buu)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;
126
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 }
142
FormJacobianLocal_K(User user,DMDALocalInfo * info,const PetscScalar u[],const PetscScalar k[],Mat Bkk)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;
147
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 }
156
FormJacobianLocal_UK(User user,DMDALocalInfo * info,DMDALocalInfo * infok,const PetscScalar u[],const PetscScalar k[],Mat Buk)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];
163
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 }
177
FormJacobianLocal_KU(User user,DMDALocalInfo * info,DMDALocalInfo * infok,const PetscScalar u[],const PetscScalar k[],Mat Bku)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);
182
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 }
197
FormJacobian_All(SNES snes,Vec X,Mat J,Mat B,PetscCtx ctx)198 static PetscErrorCode FormJacobian_All(SNES snes, Vec X, Mat J, Mat B, PetscCtx ctx)
199 {
200 User user = (User)ctx;
201 DM dau, dak;
202 DMDALocalInfo infou, infok;
203 PetscScalar *u, *k;
204 Vec Uloc, Kloc;
205
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));
261
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 }
276
FormInitial_Coupled(User user,Vec X)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;
284
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 }
301
main(int argc,char * argv[])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;
313 VecType vtype;
314
315 PetscFunctionBeginUser;
316 PetscCall(PetscInitialize(&argc, &argv, 0, help));
317 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 10, 1, 1, NULL, &dau));
318 PetscCall(DMSetOptionsPrefix(dau, "u_"));
319 PetscCall(DMSetFromOptions(dau));
320 PetscCall(DMSetUp(dau));
321 PetscCall(DMDAGetOwnershipRanges(dau, &lxu, 0, 0));
322 PetscCall(DMDAGetInfo(dau, 0, &m, 0, 0, &sizes, 0, 0, 0, 0, 0, 0, 0, 0));
323 PetscCall(PetscMalloc1(sizes, &lxk));
324 PetscCall(PetscArraycpy(lxk, lxu, sizes));
325 lxk[0]--;
326 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m - 1, 1, 1, lxk, &dak));
327 PetscCall(DMSetOptionsPrefix(dak, "k_"));
328 PetscCall(DMSetFromOptions(dak));
329 PetscCall(DMSetUp(dak));
330 PetscCall(PetscFree(lxk));
331
332 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &pack));
333 PetscCall(DMSetOptionsPrefix(pack, "pack_"));
334 PetscCall(DMCompositeAddDM(pack, dau));
335 PetscCall(DMCompositeAddDM(pack, dak));
336 PetscCall(DMDASetFieldName(dau, 0, "u"));
337 PetscCall(DMDASetFieldName(dak, 0, "k"));
338 PetscCall(DMSetFromOptions(pack));
339
340 /* if we call -pack_dm_vec_type xxx it does not get propagated to the subDMs */
341 PetscCall(DMGetVecType(pack, &vtype));
342 PetscCall(DMSetVecType(dau, vtype));
343 PetscCall(DMSetVecType(dak, vtype));
344
345 PetscCall(DMCreateGlobalVector(pack, &X));
346 PetscCall(VecDuplicate(X, &F));
347
348 PetscCall(PetscNew(&user));
349
350 user->pack = pack;
351
352 PetscCall(DMCompositeGetGlobalISs(pack, &isg));
353 PetscCall(DMCompositeGetLocalVectors(pack, &user->Uloc, &user->Kloc));
354 PetscCall(DMCompositeScatter(pack, X, user->Uloc, user->Kloc));
355
356 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Coupled problem options", "SNES");
357 {
358 user->ptype = 0;
359 pass_dm = PETSC_TRUE;
360
361 PetscCall(PetscOptionsInt("-problem_type", "0: solve for u only, 1: solve for k only, 2: solve for both", 0, user->ptype, &user->ptype, NULL));
362 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));
363 }
364 PetscOptionsEnd();
365
366 PetscCall(FormInitial_Coupled(user, X));
367
368 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
369 switch (user->ptype) {
370 case 0:
371 PetscCall(DMCompositeGetAccess(pack, X, &Xu, PETSC_NULLPTR));
372 PetscCall(DMCompositeGetAccess(pack, F, &Fu, PETSC_NULLPTR));
373 PetscCall(DMCreateMatrix(dau, &B));
374 PetscCall(SNESSetFunction(snes, Fu, FormFunction_All, user));
375 PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
376 PetscCall(SNESSetFromOptions(snes));
377 PetscCall(SNESSetDM(snes, dau));
378 PetscCall(SNESSolve(snes, NULL, Xu));
379 PetscCall(DMCompositeRestoreAccess(pack, X, &Xu, PETSC_NULLPTR));
380 PetscCall(DMCompositeRestoreAccess(pack, F, &Fu, PETSC_NULLPTR));
381 break;
382 case 1:
383 PetscCall(DMCompositeGetAccess(pack, X, PETSC_NULLPTR, &Xk));
384 PetscCall(DMCompositeGetAccess(pack, F, PETSC_NULLPTR, &Fk));
385 PetscCall(DMCreateMatrix(dak, &B));
386 PetscCall(SNESSetFunction(snes, Fk, FormFunction_All, user));
387 PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
388 PetscCall(SNESSetFromOptions(snes));
389 PetscCall(SNESSetDM(snes, dak));
390 PetscCall(SNESSolve(snes, NULL, Xk));
391 PetscCall(DMCompositeRestoreAccess(pack, X, PETSC_NULLPTR, &Xk));
392 PetscCall(DMCompositeRestoreAccess(pack, F, PETSC_NULLPTR, &Fk));
393 break;
394 case 2:
395 PetscCall(DMCreateMatrix(pack, &B));
396 /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */
397 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
398 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
399 PetscCall(SNESSetFunction(snes, F, FormFunction_All, user));
400 PetscCall(SNESSetJacobian(snes, B, B, FormJacobian_All, user));
401 PetscCall(SNESSetFromOptions(snes));
402 if (!pass_dm) { /* Manually provide index sets and names for the splits */
403 KSP ksp;
404 PC pc;
405 PetscCall(SNESGetKSP(snes, &ksp));
406 PetscCall(KSPGetPC(ksp, &pc));
407 PetscCall(PCFieldSplitSetIS(pc, "u", isg[0]));
408 PetscCall(PCFieldSplitSetIS(pc, "k", isg[1]));
409 } else {
410 /* The same names come from the options prefix for dau and dak. This option can support geometric multigrid inside
411 * of splits, but it requires using a DM (perhaps your own implementation). */
412 PetscCall(SNESSetDM(snes, pack));
413 }
414 PetscCall(SNESSolve(snes, NULL, X));
415 break;
416 }
417 PetscCall(VecViewFromOptions(X, NULL, "-view_sol"));
418
419 if (0) {
420 PetscInt col = 0;
421 PetscBool mult_dup = PETSC_FALSE, view_dup = PETSC_FALSE;
422 Mat D;
423 Vec Y;
424
425 PetscCall(PetscOptionsGetInt(NULL, 0, "-col", &col, 0));
426 PetscCall(PetscOptionsGetBool(NULL, 0, "-mult_dup", &mult_dup, 0));
427 PetscCall(PetscOptionsGetBool(NULL, 0, "-view_dup", &view_dup, 0));
428
429 PetscCall(VecDuplicate(X, &Y));
430 /* PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); */
431 /* PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); */
432 PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &D));
433 PetscCall(VecZeroEntries(X));
434 PetscCall(VecSetValue(X, col, 1.0, INSERT_VALUES));
435 PetscCall(VecAssemblyBegin(X));
436 PetscCall(VecAssemblyEnd(X));
437 PetscCall(MatMult(mult_dup ? D : B, X, Y));
438 PetscCall(MatView(view_dup ? D : B, PETSC_VIEWER_STDOUT_WORLD));
439 /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
440 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
441 PetscCall(MatDestroy(&D));
442 PetscCall(VecDestroy(&Y));
443 }
444
445 PetscCall(DMCompositeRestoreLocalVectors(pack, &user->Uloc, &user->Kloc));
446 PetscCall(PetscFree(user));
447
448 PetscCall(ISDestroy(&isg[0]));
449 PetscCall(ISDestroy(&isg[1]));
450 PetscCall(PetscFree(isg));
451 PetscCall(VecDestroy(&X));
452 PetscCall(VecDestroy(&F));
453 PetscCall(MatDestroy(&B));
454 PetscCall(DMDestroy(&dau));
455 PetscCall(DMDestroy(&dak));
456 PetscCall(DMDestroy(&pack));
457 PetscCall(SNESDestroy(&snes));
458 PetscCall(PetscFinalize());
459 return 0;
460 }
461
462 /*TEST
463
464 test:
465 suffix: 0
466 nsize: 3
467 args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 0
468
469 test:
470 suffix: 1
471 nsize: 3
472 args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 1
473
474 test:
475 suffix: 2
476 nsize: 3
477 args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 2
478
479 testset:
480 suffix: 3
481 nsize: 3
482 output_file: output/ex28_3.out
483 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
484
485 test:
486 suffix: std
487 args: -pack_dm_mat_type {{aij nest}}
488
489 test:
490 suffix: cuda
491 requires: cuda
492 args: -pack_dm_mat_type aijcusparse -pack_dm_vec_type cuda
493
494 test:
495 suffix: hip
496 requires: hip
497 args: -pack_dm_mat_type aijhipsparse -pack_dm_vec_type hip
498
499 test:
500 suffix: kok
501 requires: kokkos_kernels
502 args: -pack_dm_mat_type aijkokkos -pack_dm_vec_type kokkos
503
504 test:
505 requires: mumps
506 suffix: 3_nest_lu
507 nsize: {{1 3}}
508 output_file: output/ex28_3.out
509 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
510
511 test:
512 suffix: 4
513 nsize: 6
514 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
515 requires: !single
516
517 test:
518 suffix: glvis_composite_da_1d
519 args: -u_da_grid_x 20 -problem_type 0 -snes_monitor_solution glvis:
520
521 test:
522 suffix: cuda
523 nsize: 1
524 requires: cuda
525 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
526
527 test:
528 suffix: viennacl
529 nsize: 1
530 requires: viennacl
531 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
532
533 TEST*/
534