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