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