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(0); 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 60 ubar = 0.5*(u[i+1]+u[i]), 61 gradu = (u[i+1]-u[i])/hx, 62 g = 1. + gradu*gradu, 63 w = 1./(1.+ubar) + 1./g; 64 f[i] = hx*(PetscExpScalar(k[i]-1.0) + k[i] - 1./w); 65 } 66 PetscFunctionReturn(0); 67 } 68 69 static PetscErrorCode FormFunction_All(SNES snes,Vec X,Vec F,void *ctx) 70 { 71 User user = (User)ctx; 72 DM dau,dak; 73 DMDALocalInfo infou,infok; 74 PetscScalar *u,*k; 75 PetscScalar *fu,*fk; 76 Vec Uloc,Kloc,Fu,Fk; 77 78 PetscFunctionBeginUser; 79 CHKERRQ(DMCompositeGetEntries(user->pack,&dau,&dak)); 80 CHKERRQ(DMDAGetLocalInfo(dau,&infou)); 81 CHKERRQ(DMDAGetLocalInfo(dak,&infok)); 82 CHKERRQ(DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc)); 83 switch (user->ptype) { 84 case 0: 85 CHKERRQ(DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc)); 86 CHKERRQ(DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc)); 87 CHKERRQ(DMDAVecGetArray(dau,Uloc,&u)); 88 CHKERRQ(DMDAVecGetArray(dak,user->Kloc,&k)); 89 CHKERRQ(DMDAVecGetArray(dau,F,&fu)); 90 CHKERRQ(FormFunctionLocal_U(user,&infou,u,k,fu)); 91 CHKERRQ(DMDAVecRestoreArray(dau,F,&fu)); 92 CHKERRQ(DMDAVecRestoreArray(dau,Uloc,&u)); 93 CHKERRQ(DMDAVecRestoreArray(dak,user->Kloc,&k)); 94 break; 95 case 1: 96 CHKERRQ(DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc)); 97 CHKERRQ(DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc)); 98 CHKERRQ(DMDAVecGetArray(dau,user->Uloc,&u)); 99 CHKERRQ(DMDAVecGetArray(dak,Kloc,&k)); 100 CHKERRQ(DMDAVecGetArray(dak,F,&fk)); 101 CHKERRQ(FormFunctionLocal_K(user,&infok,u,k,fk)); 102 CHKERRQ(DMDAVecRestoreArray(dak,F,&fk)); 103 CHKERRQ(DMDAVecRestoreArray(dau,user->Uloc,&u)); 104 CHKERRQ(DMDAVecRestoreArray(dak,Kloc,&k)); 105 break; 106 case 2: 107 CHKERRQ(DMCompositeScatter(user->pack,X,Uloc,Kloc)); 108 CHKERRQ(DMDAVecGetArray(dau,Uloc,&u)); 109 CHKERRQ(DMDAVecGetArray(dak,Kloc,&k)); 110 CHKERRQ(DMCompositeGetAccess(user->pack,F,&Fu,&Fk)); 111 CHKERRQ(DMDAVecGetArray(dau,Fu,&fu)); 112 CHKERRQ(DMDAVecGetArray(dak,Fk,&fk)); 113 CHKERRQ(FormFunctionLocal_U(user,&infou,u,k,fu)); 114 CHKERRQ(FormFunctionLocal_K(user,&infok,u,k,fk)); 115 CHKERRQ(DMDAVecRestoreArray(dau,Fu,&fu)); 116 CHKERRQ(DMDAVecRestoreArray(dak,Fk,&fk)); 117 CHKERRQ(DMCompositeRestoreAccess(user->pack,F,&Fu,&Fk)); 118 CHKERRQ(DMDAVecRestoreArray(dau,Uloc,&u)); 119 CHKERRQ(DMDAVecRestoreArray(dak,Kloc,&k)); 120 break; 121 } 122 CHKERRQ(DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc)); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode FormJacobianLocal_U(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Buu) 127 { 128 PetscReal hx = 1./info->mx; 129 PetscInt i; 130 131 PetscFunctionBeginUser; 132 for (i=info->xs; i<info->xs+info->xm; i++) { 133 PetscInt row = i-info->gxs,cols[] = {row-1,row,row+1}; 134 PetscScalar val = 1./hx; 135 if (i == 0) { 136 CHKERRQ(MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES)); 137 } else if (i == info->mx-1) { 138 CHKERRQ(MatSetValuesLocal(Buu,1,&row,1,&row,&val,INSERT_VALUES)); 139 } else { 140 PetscScalar vals[] = {-k[i-1]/hx,(k[i-1]+k[i])/hx,-k[i]/hx}; 141 CHKERRQ(MatSetValuesLocal(Buu,1,&row,3,cols,vals,INSERT_VALUES)); 142 } 143 } 144 PetscFunctionReturn(0); 145 } 146 147 static PetscErrorCode FormJacobianLocal_K(User user,DMDALocalInfo *info,const PetscScalar u[],const PetscScalar k[],Mat Bkk) 148 { 149 PetscReal hx = 1./info->mx; 150 PetscInt i; 151 152 PetscFunctionBeginUser; 153 for (i=info->xs; i<info->xs+info->xm; i++) { 154 PetscInt row = i-info->gxs; 155 PetscScalar vals[] = {hx*(PetscExpScalar(k[i]-1.)+ (PetscScalar)1.)}; 156 CHKERRQ(MatSetValuesLocal(Bkk,1,&row,1,&row,vals,INSERT_VALUES)); 157 } 158 PetscFunctionReturn(0); 159 } 160 161 static PetscErrorCode FormJacobianLocal_UK(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Buk) 162 { 163 PetscReal hx = 1./info->mx; 164 PetscInt i; 165 PetscInt row,cols[2]; 166 PetscScalar vals[2]; 167 168 PetscFunctionBeginUser; 169 if (!Buk) PetscFunctionReturn(0); /* Not assembling this block */ 170 for (i=info->xs; i<info->xs+info->xm; i++) { 171 if (i == 0 || i == info->mx-1) continue; 172 row = i-info->gxs; 173 cols[0] = i-1-infok->gxs; vals[0] = (u[i]-u[i-1])/hx; 174 cols[1] = i-infok->gxs; vals[1] = (u[i]-u[i+1])/hx; 175 CHKERRQ(MatSetValuesLocal(Buk,1,&row,2,cols,vals,INSERT_VALUES)); 176 } 177 PetscFunctionReturn(0); 178 } 179 180 static PetscErrorCode FormJacobianLocal_KU(User user,DMDALocalInfo *info,DMDALocalInfo *infok,const PetscScalar u[],const PetscScalar k[],Mat Bku) 181 { 182 PetscInt i; 183 PetscReal hx = 1./(info->mx-1); 184 185 PetscFunctionBeginUser; 186 if (!Bku) PetscFunctionReturn(0); /* Not assembling this block */ 187 for (i=infok->xs; i<infok->xs+infok->xm; i++) { 188 PetscInt row = i-infok->gxs,cols[2]; 189 PetscScalar vals[2]; 190 const PetscScalar 191 ubar = 0.5*(u[i]+u[i+1]), 192 ubar_L = 0.5, 193 ubar_R = 0.5, 194 gradu = (u[i+1]-u[i])/hx, 195 gradu_L = -1./hx, 196 gradu_R = 1./hx, 197 g = 1. + PetscSqr(gradu), 198 g_gradu = 2.*gradu, 199 w = 1./(1.+ubar) + 1./g, 200 w_ubar = -1./PetscSqr(1.+ubar), 201 w_gradu = -g_gradu/PetscSqr(g), 202 iw = 1./w, 203 iw_ubar = -w_ubar * PetscSqr(iw), 204 iw_gradu = -w_gradu * PetscSqr(iw); 205 cols[0] = i-info->gxs; vals[0] = -hx*(iw_ubar*ubar_L + iw_gradu*gradu_L); 206 cols[1] = i+1-info->gxs; vals[1] = -hx*(iw_ubar*ubar_R + iw_gradu*gradu_R); 207 CHKERRQ(MatSetValuesLocal(Bku,1,&row,2,cols,vals,INSERT_VALUES)); 208 } 209 PetscFunctionReturn(0); 210 } 211 212 static PetscErrorCode FormJacobian_All(SNES snes,Vec X,Mat J,Mat B,void *ctx) 213 { 214 User user = (User)ctx; 215 DM dau,dak; 216 DMDALocalInfo infou,infok; 217 PetscScalar *u,*k; 218 Vec Uloc,Kloc; 219 220 PetscFunctionBeginUser; 221 CHKERRQ(DMCompositeGetEntries(user->pack,&dau,&dak)); 222 CHKERRQ(DMDAGetLocalInfo(dau,&infou)); 223 CHKERRQ(DMDAGetLocalInfo(dak,&infok)); 224 CHKERRQ(DMCompositeGetLocalVectors(user->pack,&Uloc,&Kloc)); 225 switch (user->ptype) { 226 case 0: 227 CHKERRQ(DMGlobalToLocalBegin(dau,X,INSERT_VALUES,Uloc)); 228 CHKERRQ(DMGlobalToLocalEnd (dau,X,INSERT_VALUES,Uloc)); 229 CHKERRQ(DMDAVecGetArray(dau,Uloc,&u)); 230 CHKERRQ(DMDAVecGetArray(dak,user->Kloc,&k)); 231 CHKERRQ(FormJacobianLocal_U(user,&infou,u,k,B)); 232 CHKERRQ(DMDAVecRestoreArray(dau,Uloc,&u)); 233 CHKERRQ(DMDAVecRestoreArray(dak,user->Kloc,&k)); 234 break; 235 case 1: 236 CHKERRQ(DMGlobalToLocalBegin(dak,X,INSERT_VALUES,Kloc)); 237 CHKERRQ(DMGlobalToLocalEnd (dak,X,INSERT_VALUES,Kloc)); 238 CHKERRQ(DMDAVecGetArray(dau,user->Uloc,&u)); 239 CHKERRQ(DMDAVecGetArray(dak,Kloc,&k)); 240 CHKERRQ(FormJacobianLocal_K(user,&infok,u,k,B)); 241 CHKERRQ(DMDAVecRestoreArray(dau,user->Uloc,&u)); 242 CHKERRQ(DMDAVecRestoreArray(dak,Kloc,&k)); 243 break; 244 case 2: { 245 Mat Buu,Buk,Bku,Bkk; 246 PetscBool nest; 247 IS *is; 248 CHKERRQ(DMCompositeScatter(user->pack,X,Uloc,Kloc)); 249 CHKERRQ(DMDAVecGetArray(dau,Uloc,&u)); 250 CHKERRQ(DMDAVecGetArray(dak,Kloc,&k)); 251 CHKERRQ(DMCompositeGetLocalISs(user->pack,&is)); 252 CHKERRQ(MatGetLocalSubMatrix(B,is[0],is[0],&Buu)); 253 CHKERRQ(MatGetLocalSubMatrix(B,is[0],is[1],&Buk)); 254 CHKERRQ(MatGetLocalSubMatrix(B,is[1],is[0],&Bku)); 255 CHKERRQ(MatGetLocalSubMatrix(B,is[1],is[1],&Bkk)); 256 CHKERRQ(FormJacobianLocal_U(user,&infou,u,k,Buu)); 257 CHKERRQ(PetscObjectTypeCompare((PetscObject)B,MATNEST,&nest)); 258 if (!nest) { 259 /* 260 DMCreateMatrix_Composite() for a nested matrix does not generate off-block matrices that one can call MatSetValuesLocal() on, it just creates dummy 261 matrices with no entries; there cannot insert values into them. Commit b6480e041dd2293a65f96222772d68cdb4ed6306 262 changed Mat_Nest() from returning NULL pointers for these submatrices to dummy matrices because PCFIELDSPLIT could not 263 handle the returned null matrices. 264 */ 265 CHKERRQ(FormJacobianLocal_UK(user,&infou,&infok,u,k,Buk)); 266 CHKERRQ(FormJacobianLocal_KU(user,&infou,&infok,u,k,Bku)); 267 } 268 CHKERRQ(FormJacobianLocal_K(user,&infok,u,k,Bkk)); 269 CHKERRQ(MatRestoreLocalSubMatrix(B,is[0],is[0],&Buu)); 270 CHKERRQ(MatRestoreLocalSubMatrix(B,is[0],is[1],&Buk)); 271 CHKERRQ(MatRestoreLocalSubMatrix(B,is[1],is[0],&Bku)); 272 CHKERRQ(MatRestoreLocalSubMatrix(B,is[1],is[1],&Bkk)); 273 CHKERRQ(DMDAVecRestoreArray(dau,Uloc,&u)); 274 CHKERRQ(DMDAVecRestoreArray(dak,Kloc,&k)); 275 276 CHKERRQ(ISDestroy(&is[0])); 277 CHKERRQ(ISDestroy(&is[1])); 278 CHKERRQ(PetscFree(is)); 279 } break; 280 } 281 CHKERRQ(DMCompositeRestoreLocalVectors(user->pack,&Uloc,&Kloc)); 282 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 283 CHKERRQ(MatAssemblyEnd (B,MAT_FINAL_ASSEMBLY)); 284 if (J != B) { 285 CHKERRQ(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 286 CHKERRQ(MatAssemblyEnd (J,MAT_FINAL_ASSEMBLY)); 287 } 288 PetscFunctionReturn(0); 289 } 290 291 static PetscErrorCode FormInitial_Coupled(User user,Vec X) 292 { 293 DM dau,dak; 294 DMDALocalInfo infou,infok; 295 Vec Xu,Xk; 296 PetscScalar *u,*k,hx; 297 PetscInt i; 298 299 PetscFunctionBeginUser; 300 CHKERRQ(DMCompositeGetEntries(user->pack,&dau,&dak)); 301 CHKERRQ(DMCompositeGetAccess(user->pack,X,&Xu,&Xk)); 302 CHKERRQ(DMDAVecGetArray(dau,Xu,&u)); 303 CHKERRQ(DMDAVecGetArray(dak,Xk,&k)); 304 CHKERRQ(DMDAGetLocalInfo(dau,&infou)); 305 CHKERRQ(DMDAGetLocalInfo(dak,&infok)); 306 hx = 1./(infok.mx); 307 for (i=infou.xs; i<infou.xs+infou.xm; i++) u[i] = (PetscScalar)i*hx * (1.-(PetscScalar)i*hx); 308 for (i=infok.xs; i<infok.xs+infok.xm; i++) k[i] = 1.0 + 0.5*PetscSinScalar(2*PETSC_PI*i*hx); 309 CHKERRQ(DMDAVecRestoreArray(dau,Xu,&u)); 310 CHKERRQ(DMDAVecRestoreArray(dak,Xk,&k)); 311 CHKERRQ(DMCompositeRestoreAccess(user->pack,X,&Xu,&Xk)); 312 CHKERRQ(DMCompositeScatter(user->pack,X,user->Uloc,user->Kloc)); 313 PetscFunctionReturn(0); 314 } 315 316 int main(int argc, char *argv[]) 317 { 318 PetscErrorCode ierr; 319 DM dau,dak,pack; 320 const PetscInt *lxu; 321 PetscInt *lxk,m,sizes; 322 User user; 323 SNES snes; 324 Vec X,F,Xu,Xk,Fu,Fk; 325 Mat B; 326 IS *isg; 327 PetscBool pass_dm; 328 329 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 330 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,10,1,1,NULL,&dau)); 331 CHKERRQ(DMSetOptionsPrefix(dau,"u_")); 332 CHKERRQ(DMSetFromOptions(dau)); 333 CHKERRQ(DMSetUp(dau)); 334 CHKERRQ(DMDAGetOwnershipRanges(dau,&lxu,0,0)); 335 CHKERRQ(DMDAGetInfo(dau,0, &m,0,0, &sizes,0,0, 0,0,0,0,0,0)); 336 CHKERRQ(PetscMalloc1(sizes,&lxk)); 337 CHKERRQ(PetscArraycpy(lxk,lxu,sizes)); 338 lxk[0]--; 339 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m-1,1,1,lxk,&dak)); 340 CHKERRQ(DMSetOptionsPrefix(dak,"k_")); 341 CHKERRQ(DMSetFromOptions(dak)); 342 CHKERRQ(DMSetUp(dak)); 343 CHKERRQ(PetscFree(lxk)); 344 345 CHKERRQ(DMCompositeCreate(PETSC_COMM_WORLD,&pack)); 346 CHKERRQ(DMSetOptionsPrefix(pack,"pack_")); 347 CHKERRQ(DMCompositeAddDM(pack,dau)); 348 CHKERRQ(DMCompositeAddDM(pack,dak)); 349 CHKERRQ(DMDASetFieldName(dau,0,"u")); 350 CHKERRQ(DMDASetFieldName(dak,0,"k")); 351 CHKERRQ(DMSetFromOptions(pack)); 352 353 CHKERRQ(DMCreateGlobalVector(pack,&X)); 354 CHKERRQ(VecDuplicate(X,&F)); 355 356 CHKERRQ(PetscNew(&user)); 357 358 user->pack = pack; 359 360 CHKERRQ(DMCompositeGetGlobalISs(pack,&isg)); 361 CHKERRQ(DMCompositeGetLocalVectors(pack,&user->Uloc,&user->Kloc)); 362 CHKERRQ(DMCompositeScatter(pack,X,user->Uloc,user->Kloc)); 363 364 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Coupled problem options","SNES");CHKERRQ(ierr); 365 { 366 user->ptype = 0; pass_dm = PETSC_TRUE; 367 368 CHKERRQ(PetscOptionsInt("-problem_type","0: solve for u only, 1: solve for k only, 2: solve for both",0,user->ptype,&user->ptype,NULL)); 369 CHKERRQ(PetscOptionsBool("-pass_dm","Pass the packed DM to SNES to use when determining splits and forward into splits",0,pass_dm,&pass_dm,NULL)); 370 } 371 ierr = PetscOptionsEnd();CHKERRQ(ierr); 372 373 CHKERRQ(FormInitial_Coupled(user,X)); 374 375 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 376 switch (user->ptype) { 377 case 0: 378 CHKERRQ(DMCompositeGetAccess(pack,X,&Xu,0)); 379 CHKERRQ(DMCompositeGetAccess(pack,F,&Fu,0)); 380 CHKERRQ(DMCreateMatrix(dau,&B)); 381 CHKERRQ(SNESSetFunction(snes,Fu,FormFunction_All,user)); 382 CHKERRQ(SNESSetJacobian(snes,B,B,FormJacobian_All,user)); 383 CHKERRQ(SNESSetFromOptions(snes)); 384 CHKERRQ(SNESSetDM(snes,dau)); 385 CHKERRQ(SNESSolve(snes,NULL,Xu)); 386 CHKERRQ(DMCompositeRestoreAccess(pack,X,&Xu,0)); 387 CHKERRQ(DMCompositeRestoreAccess(pack,F,&Fu,0)); 388 break; 389 case 1: 390 CHKERRQ(DMCompositeGetAccess(pack,X,0,&Xk)); 391 CHKERRQ(DMCompositeGetAccess(pack,F,0,&Fk)); 392 CHKERRQ(DMCreateMatrix(dak,&B)); 393 CHKERRQ(SNESSetFunction(snes,Fk,FormFunction_All,user)); 394 CHKERRQ(SNESSetJacobian(snes,B,B,FormJacobian_All,user)); 395 CHKERRQ(SNESSetFromOptions(snes)); 396 CHKERRQ(SNESSetDM(snes,dak)); 397 CHKERRQ(SNESSolve(snes,NULL,Xk)); 398 CHKERRQ(DMCompositeRestoreAccess(pack,X,0,&Xk)); 399 CHKERRQ(DMCompositeRestoreAccess(pack,F,0,&Fk)); 400 break; 401 case 2: 402 CHKERRQ(DMCreateMatrix(pack,&B)); 403 /* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */ 404 CHKERRQ(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE)); 405 CHKERRQ(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 406 CHKERRQ(SNESSetFunction(snes,F,FormFunction_All,user)); 407 CHKERRQ(SNESSetJacobian(snes,B,B,FormJacobian_All,user)); 408 CHKERRQ(SNESSetFromOptions(snes)); 409 if (!pass_dm) { /* Manually provide index sets and names for the splits */ 410 KSP ksp; 411 PC pc; 412 CHKERRQ(SNESGetKSP(snes,&ksp)); 413 CHKERRQ(KSPGetPC(ksp,&pc)); 414 CHKERRQ(PCFieldSplitSetIS(pc,"u",isg[0])); 415 CHKERRQ(PCFieldSplitSetIS(pc,"k",isg[1])); 416 } else { 417 /* The same names come from the options prefix for dau and dak. This option can support geometric multigrid inside 418 * of splits, but it requires using a DM (perhaps your own implementation). */ 419 CHKERRQ(SNESSetDM(snes,pack)); 420 } 421 CHKERRQ(SNESSolve(snes,NULL,X)); 422 break; 423 } 424 CHKERRQ(VecViewFromOptions(X,NULL,"-view_sol")); 425 426 if (0) { 427 PetscInt col = 0; 428 PetscBool mult_dup = PETSC_FALSE,view_dup = PETSC_FALSE; 429 Mat D; 430 Vec Y; 431 432 CHKERRQ(PetscOptionsGetInt(NULL,0,"-col",&col,0)); 433 CHKERRQ(PetscOptionsGetBool(NULL,0,"-mult_dup",&mult_dup,0)); 434 CHKERRQ(PetscOptionsGetBool(NULL,0,"-view_dup",&view_dup,0)); 435 436 CHKERRQ(VecDuplicate(X,&Y)); 437 /* CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); */ 438 /* CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); */ 439 CHKERRQ(MatConvert(B,MATAIJ,MAT_INITIAL_MATRIX,&D)); 440 CHKERRQ(VecZeroEntries(X)); 441 CHKERRQ(VecSetValue(X,col,1.0,INSERT_VALUES)); 442 CHKERRQ(VecAssemblyBegin(X)); 443 CHKERRQ(VecAssemblyEnd(X)); 444 CHKERRQ(MatMult(mult_dup ? D : B,X,Y)); 445 CHKERRQ(MatView(view_dup ? D : B,PETSC_VIEWER_STDOUT_WORLD)); 446 /* CHKERRQ(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */ 447 CHKERRQ(VecView(Y,PETSC_VIEWER_STDOUT_WORLD)); 448 CHKERRQ(MatDestroy(&D)); 449 CHKERRQ(VecDestroy(&Y)); 450 } 451 452 CHKERRQ(DMCompositeRestoreLocalVectors(pack,&user->Uloc,&user->Kloc)); 453 CHKERRQ(PetscFree(user)); 454 455 CHKERRQ(ISDestroy(&isg[0])); 456 CHKERRQ(ISDestroy(&isg[1])); 457 CHKERRQ(PetscFree(isg)); 458 CHKERRQ(VecDestroy(&X)); 459 CHKERRQ(VecDestroy(&F)); 460 CHKERRQ(MatDestroy(&B)); 461 CHKERRQ(DMDestroy(&dau)); 462 CHKERRQ(DMDestroy(&dak)); 463 CHKERRQ(DMDestroy(&pack)); 464 CHKERRQ(SNESDestroy(&snes)); 465 ierr = PetscFinalize(); 466 return ierr; 467 } 468 469 /*TEST 470 471 test: 472 suffix: 0 473 nsize: 3 474 args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 0 475 476 test: 477 suffix: 1 478 nsize: 3 479 args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 1 480 481 test: 482 suffix: 2 483 nsize: 3 484 args: -u_da_grid_x 20 -snes_converged_reason -snes_monitor_short -problem_type 2 485 486 test: 487 suffix: 3 488 nsize: 3 489 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 490 491 test: 492 suffix: 4 493 nsize: 6 494 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 495 requires: !single 496 497 test: 498 suffix: glvis_composite_da_1d 499 args: -u_da_grid_x 20 -problem_type 0 -snes_monitor_solution glvis: 500 501 test: 502 suffix: cuda 503 nsize: 1 504 requires: cuda 505 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 506 507 test: 508 suffix: viennacl 509 nsize: 1 510 requires: viennacl 511 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 512 513 TEST*/ 514