1 2 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMView_DA_2d" 6 PetscErrorCode DMView_DA_2d(DM da,PetscViewer viewer) 7 { 8 PetscErrorCode ierr; 9 PetscMPIInt rank; 10 PetscBool iascii,isdraw,isbinary; 11 DM_DA *dd = (DM_DA*)da->data; 12 #if defined(PETSC_HAVE_MATLAB_ENGINE) 13 PetscBool ismatlab; 14 #endif 15 16 PetscFunctionBegin; 17 ierr = MPI_Comm_rank(((PetscObject)da)->comm,&rank);CHKERRQ(ierr); 18 19 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 20 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 21 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 22 #if defined(PETSC_HAVE_MATLAB_ENGINE) 23 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);CHKERRQ(ierr); 24 #endif 25 if (iascii) { 26 PetscViewerFormat format; 27 28 ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr); 29 if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) { 30 DMDALocalInfo info; 31 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 32 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 33 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D m %D n %D w %D s %D\n",rank,dd->M,dd->N,dd->m,dd->n,dd->w,dd->s);CHKERRQ(ierr); 34 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\n",info.xs,info.xs+info.xm,info.ys,info.ys+info.ym);CHKERRQ(ierr); 35 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 36 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 37 } else { 38 ierr = DMView_DA_VTK(da,viewer);CHKERRQ(ierr); 39 } 40 } else if (isdraw) { 41 PetscDraw draw; 42 double ymin = -1*dd->s-1,ymax = dd->N+dd->s; 43 double xmin = -1*dd->s-1,xmax = dd->M+dd->s; 44 double x,y; 45 PetscInt base,*idx; 46 char node[10]; 47 PetscBool isnull; 48 49 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 50 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 51 if (!dd->coordinates) { 52 ierr = PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);CHKERRQ(ierr); 53 } 54 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 55 56 /* first processor draw all node lines */ 57 if (!rank) { 58 ymin = 0.0; ymax = dd->N - 1; 59 for (xmin=0; xmin<dd->M; xmin++) { 60 ierr = PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);CHKERRQ(ierr); 61 } 62 xmin = 0.0; xmax = dd->M - 1; 63 for (ymin=0; ymin<dd->N; ymin++) { 64 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);CHKERRQ(ierr); 65 } 66 } 67 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 68 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 69 70 /* draw my box */ 71 ymin = dd->ys; ymax = dd->ye - 1; xmin = dd->xs/dd->w; 72 xmax =(dd->xe-1)/dd->w; 73 ierr = PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);CHKERRQ(ierr); 74 ierr = PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 75 ierr = PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 76 ierr = PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);CHKERRQ(ierr); 77 78 /* put in numbers */ 79 base = (dd->base)/dd->w; 80 for (y=ymin; y<=ymax; y++) { 81 for (x=xmin; x<=xmax; x++) { 82 sprintf(node,"%d",(int)base++); 83 ierr = PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);CHKERRQ(ierr); 84 } 85 } 86 87 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 88 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 89 /* overlay ghost numbers, useful for error checking */ 90 /* put in numbers */ 91 92 base = 0; idx = dd->idx; 93 ymin = dd->Ys; ymax = dd->Ye; xmin = dd->Xs; xmax = dd->Xe; 94 for (y=ymin; y<ymax; y++) { 95 for (x=xmin; x<xmax; x++) { 96 if ((base % dd->w) == 0) { 97 sprintf(node,"%d",(int)(idx[base]/dd->w)); 98 ierr = PetscDrawString(draw,x/dd->w,y,PETSC_DRAW_BLUE,node);CHKERRQ(ierr); 99 } 100 base++; 101 } 102 } 103 ierr = PetscDrawSynchronizedFlush(draw);CHKERRQ(ierr); 104 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 105 } else if (isbinary){ 106 ierr = DMView_DA_Binary(da,viewer);CHKERRQ(ierr); 107 #if defined(PETSC_HAVE_MATLAB_ENGINE) 108 } else if (ismatlab) { 109 ierr = DMView_DA_Matlab(da,viewer);CHKERRQ(ierr); 110 #endif 111 } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for DMDA 1d",((PetscObject)viewer)->type_name); 112 PetscFunctionReturn(0); 113 } 114 115 /* 116 M is number of grid points 117 m is number of processors 118 119 */ 120 #undef __FUNCT__ 121 #define __FUNCT__ "DMDASplitComm2d" 122 PetscErrorCode DMDASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm) 123 { 124 PetscErrorCode ierr; 125 PetscInt m,n = 0,x = 0,y = 0; 126 PetscMPIInt size,csize,rank; 127 128 PetscFunctionBegin; 129 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 130 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 131 132 csize = 4*size; 133 do { 134 if (csize % 4) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot split communicator of size %d tried %d %D %D",size,csize,x,y); 135 csize = csize/4; 136 137 m = (PetscInt)(0.5 + sqrt(((double)M)*((double)csize)/((double)N))); 138 if (!m) m = 1; 139 while (m > 0) { 140 n = csize/m; 141 if (m*n == csize) break; 142 m--; 143 } 144 if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;} 145 146 x = M/m + ((M % m) > ((csize-1) % m)); 147 y = (N + (csize-1)/m)/n; 148 } while ((x < 4 || y < 4) && csize > 1); 149 if (size != csize) { 150 MPI_Group entire_group,sub_group; 151 PetscMPIInt i,*groupies; 152 153 ierr = MPI_Comm_group(comm,&entire_group);CHKERRQ(ierr); 154 ierr = PetscMalloc(csize*sizeof(PetscInt),&groupies);CHKERRQ(ierr); 155 for (i=0; i<csize; i++) { 156 groupies[i] = (rank/csize)*csize + i; 157 } 158 ierr = MPI_Group_incl(entire_group,csize,groupies,&sub_group);CHKERRQ(ierr); 159 ierr = PetscFree(groupies);CHKERRQ(ierr); 160 ierr = MPI_Comm_create(comm,sub_group,outcomm);CHKERRQ(ierr); 161 ierr = MPI_Group_free(&entire_group);CHKERRQ(ierr); 162 ierr = MPI_Group_free(&sub_group);CHKERRQ(ierr); 163 ierr = PetscInfo1(0,"DMDASplitComm2d:Creating redundant coarse problems of size %d\n",csize);CHKERRQ(ierr); 164 } else { 165 *outcomm = comm; 166 } 167 PetscFunctionReturn(0); 168 } 169 170 #undef __FUNCT__ 171 #define __FUNCT__ "DMDAFunction" 172 static PetscErrorCode DMDAFunction(DM dm,Vec x,Vec F) 173 { 174 PetscErrorCode ierr; 175 Vec localX; 176 177 PetscFunctionBegin; 178 ierr = DMGetLocalVector(dm,&localX);CHKERRQ(ierr); 179 ierr = DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);CHKERRQ(ierr); 180 ierr = DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);CHKERRQ(ierr); 181 ierr = DMDAComputeFunction1(dm,localX,F,dm->ctx);CHKERRQ(ierr); 182 ierr = DMRestoreLocalVector(dm,&localX);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "DMDASetLocalFunction" 188 /*@C 189 DMDASetLocalFunction - Caches in a DM a local function. 190 191 Logically Collective on DMDA 192 193 Input Parameter: 194 + da - initial distributed array 195 - lf - the local function 196 197 Level: intermediate 198 199 Notes: 200 If you used SNESSetFunction(snes,r,SNESDMDAComputeFunction,ctx); then the context passed to your function is the ctx set here. 201 202 If you use SNESSetDM() and did not use SNESSetFunction() then the context passed to your function is the context set with DMSetApplicationContext() 203 204 Developer Notes: It is possibly confusing which context is passed to the user function, it would be nice to unify them somehow. 205 206 .keywords: distributed array, refine 207 208 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunctioni() 209 @*/ 210 PetscErrorCode DMDASetLocalFunction(DM da,DMDALocalFunction1 lf) 211 { 212 PetscErrorCode ierr; 213 DM_DA *dd = (DM_DA*)da->data; 214 215 PetscFunctionBegin; 216 PetscValidHeaderSpecific(da,DM_CLASSID,1); 217 ierr = DMSetFunction(da,DMDAFunction);CHKERRQ(ierr); 218 dd->lf = lf; 219 PetscFunctionReturn(0); 220 } 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "DMDASetLocalFunctioni" 224 /*@C 225 DMDASetLocalFunctioni - Caches in a DM a local function that evaluates a single component 226 227 Logically Collective on DMDA 228 229 Input Parameter: 230 + da - initial distributed array 231 - lfi - the local function 232 233 Level: intermediate 234 235 .keywords: distributed array, refine 236 237 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction() 238 @*/ 239 PetscErrorCode DMDASetLocalFunctioni(DM da,PetscErrorCode (*lfi)(DMDALocalInfo*,MatStencil*,void*,PetscScalar*,void*)) 240 { 241 DM_DA *dd = (DM_DA*)da->data; 242 PetscFunctionBegin; 243 PetscValidHeaderSpecific(da,DM_CLASSID,1); 244 dd->lfi = lfi; 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "DMDASetLocalFunctionib" 250 /*@C 251 DMDASetLocalFunctionib - Caches in a DM a block local function that evaluates a single component 252 253 Logically Collective on DMDA 254 255 Input Parameter: 256 + da - initial distributed array 257 - lfi - the local function 258 259 Level: intermediate 260 261 .keywords: distributed array, refine 262 263 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction() 264 @*/ 265 PetscErrorCode DMDASetLocalFunctionib(DM da,PetscErrorCode (*lfi)(DMDALocalInfo*,MatStencil*,void*,PetscScalar*,void*)) 266 { 267 DM_DA *dd = (DM_DA*)da->data; 268 PetscFunctionBegin; 269 PetscValidHeaderSpecific(da,DM_CLASSID,1); 270 dd->lfib = lfi; 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "DMDASetLocalAdicFunction_Private" 276 PetscErrorCode DMDASetLocalAdicFunction_Private(DM da,DMDALocalFunction1 ad_lf) 277 { 278 DM_DA *dd = (DM_DA*)da->data; 279 PetscFunctionBegin; 280 PetscValidHeaderSpecific(da,DM_CLASSID,1); 281 dd->adic_lf = ad_lf; 282 PetscFunctionReturn(0); 283 } 284 285 /*MC 286 DMDASetLocalAdicFunctioni - Caches in a DM a local functioni computed by ADIC/ADIFOR 287 288 Synopsis: 289 PetscErrorCode DMDASetLocalAdicFunctioni(DM da,PetscInt (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*) 290 291 Logically Collective on DMDA 292 293 Input Parameter: 294 + da - initial distributed array 295 - ad_lfi - the local function as computed by ADIC/ADIFOR 296 297 Level: intermediate 298 299 .keywords: distributed array, refine 300 301 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(), 302 DMDASetLocalJacobian(), DMDASetLocalFunctioni() 303 M*/ 304 305 #undef __FUNCT__ 306 #define __FUNCT__ "DMDASetLocalAdicFunctioni_Private" 307 PetscErrorCode DMDASetLocalAdicFunctioni_Private(DM da,PetscErrorCode (*ad_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*)) 308 { 309 DM_DA *dd = (DM_DA*)da->data; 310 PetscFunctionBegin; 311 PetscValidHeaderSpecific(da,DM_CLASSID,1); 312 dd->adic_lfi = ad_lfi; 313 PetscFunctionReturn(0); 314 } 315 316 /*MC 317 DMDASetLocalAdicMFFunctioni - Caches in a DM a local functioni computed by ADIC/ADIFOR 318 319 Synopsis: 320 PetscErrorCode DMDASetLocalAdicFunctioni(DM da,int (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*) 321 322 Logically Collective on DMDA 323 324 Input Parameter: 325 + da - initial distributed array 326 - admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR 327 328 Level: intermediate 329 330 .keywords: distributed array, refine 331 332 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(), 333 DMDASetLocalJacobian(), DMDASetLocalFunctioni() 334 M*/ 335 336 #undef __FUNCT__ 337 #define __FUNCT__ "DMDASetLocalAdicMFFunctioni_Private" 338 PetscErrorCode DMDASetLocalAdicMFFunctioni_Private(DM da,PetscErrorCode (*admf_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*)) 339 { 340 DM_DA *dd = (DM_DA*)da->data; 341 PetscFunctionBegin; 342 PetscValidHeaderSpecific(da,DM_CLASSID,1); 343 dd->adicmf_lfi = admf_lfi; 344 PetscFunctionReturn(0); 345 } 346 347 /*MC 348 DMDASetLocalAdicFunctionib - Caches in a DM a block local functioni computed by ADIC/ADIFOR 349 350 Synopsis: 351 PetscErrorCode DMDASetLocalAdicFunctionib(DM da,PetscInt (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*) 352 353 Logically Collective on DMDA 354 355 Input Parameter: 356 + da - initial distributed array 357 - ad_lfi - the local function as computed by ADIC/ADIFOR 358 359 Level: intermediate 360 361 .keywords: distributed array, refine 362 363 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(), 364 DMDASetLocalJacobian(), DMDASetLocalFunctionib() 365 M*/ 366 367 #undef __FUNCT__ 368 #define __FUNCT__ "DMDASetLocalAdicFunctionib_Private" 369 PetscErrorCode DMDASetLocalAdicFunctionib_Private(DM da,PetscErrorCode (*ad_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*)) 370 { 371 DM_DA *dd = (DM_DA*)da->data; 372 PetscFunctionBegin; 373 PetscValidHeaderSpecific(da,DM_CLASSID,1); 374 dd->adic_lfib = ad_lfi; 375 PetscFunctionReturn(0); 376 } 377 378 /*MC 379 DMDASetLocalAdicMFFunctionib - Caches in a DM a block local functioni computed by ADIC/ADIFOR 380 381 Synopsis: 382 PetscErrorCode DMDASetLocalAdicFunctionib(DM da,int (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*) 383 384 Logically Collective on DMDA 385 386 Input Parameter: 387 + da - initial distributed array 388 - admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR 389 390 Level: intermediate 391 392 .keywords: distributed array, refine 393 394 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(), 395 DMDASetLocalJacobian(), DMDASetLocalFunctionib() 396 M*/ 397 398 #undef __FUNCT__ 399 #define __FUNCT__ "DMDASetLocalAdicMFFunctionib_Private" 400 PetscErrorCode DMDASetLocalAdicMFFunctionib_Private(DM da,PetscErrorCode (*admf_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*)) 401 { 402 DM_DA *dd = (DM_DA*)da->data; 403 PetscFunctionBegin; 404 PetscValidHeaderSpecific(da,DM_CLASSID,1); 405 dd->adicmf_lfib = admf_lfi; 406 PetscFunctionReturn(0); 407 } 408 409 /*MC 410 DMDASetLocalAdicMFFunction - Caches in a DM a local function computed by ADIC/ADIFOR 411 412 Synopsis: 413 PetscErrorCode DMDASetLocalAdicMFFunction(DM da,DMDALocalFunction1 ad_lf) 414 415 Logically Collective on DMDA 416 417 Input Parameter: 418 + da - initial distributed array 419 - ad_lf - the local function as computed by ADIC/ADIFOR 420 421 Level: intermediate 422 423 .keywords: distributed array, refine 424 425 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(), 426 DMDASetLocalJacobian() 427 M*/ 428 429 #undef __FUNCT__ 430 #define __FUNCT__ "DMDASetLocalAdicMFFunction_Private" 431 PetscErrorCode DMDASetLocalAdicMFFunction_Private(DM da,DMDALocalFunction1 ad_lf) 432 { 433 DM_DA *dd = (DM_DA*)da->data; 434 PetscFunctionBegin; 435 PetscValidHeaderSpecific(da,DM_CLASSID,1); 436 dd->adicmf_lf = ad_lf; 437 PetscFunctionReturn(0); 438 } 439 440 #undef __FUNCT__ 441 #define __FUNCT__ "DMDAJacobianDefaultLocal" 442 PetscErrorCode DMDAJacobianLocal(DM dm,Vec x,Mat A,Mat B, MatStructure *str) 443 { 444 PetscErrorCode ierr; 445 Vec localX; 446 447 PetscFunctionBegin; 448 ierr = DMGetLocalVector(dm,&localX);CHKERRQ(ierr); 449 ierr = DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);CHKERRQ(ierr); 450 ierr = DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);CHKERRQ(ierr); 451 ierr = MatFDColoringApply(B,dm->fd,localX,str,dm);CHKERRQ(ierr); 452 ierr = DMRestoreLocalVector(dm,&localX);CHKERRQ(ierr); 453 /* Assemble true Jacobian; if it is different */ 454 if (A != B) { 455 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 456 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 457 } 458 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 459 *str = SAME_NONZERO_PATTERN; 460 PetscFunctionReturn(0); 461 } 462 463 464 #undef __FUNCT__ 465 #define __FUNCT__ "DMDAJacobian" 466 static PetscErrorCode DMDAJacobian(DM dm,Vec x,Mat A,Mat B, MatStructure *str) 467 { 468 PetscErrorCode ierr; 469 Vec localX; 470 471 PetscFunctionBegin; 472 ierr = DMGetLocalVector(dm,&localX);CHKERRQ(ierr); 473 ierr = DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);CHKERRQ(ierr); 474 ierr = DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);CHKERRQ(ierr); 475 ierr = DMDAComputeJacobian1(dm,localX,B,dm->ctx);CHKERRQ(ierr); 476 ierr = DMRestoreLocalVector(dm,&localX);CHKERRQ(ierr); 477 /* Assemble true Jacobian; if it is different */ 478 if (A != B) { 479 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 480 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 481 } 482 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 483 *str = SAME_NONZERO_PATTERN; 484 PetscFunctionReturn(0); 485 } 486 487 /*@C 488 DMDASetLocalJacobian - Caches in a DM a local Jacobian computation function 489 490 Logically Collective on DMDA 491 492 493 Input Parameter: 494 + da - initial distributed array 495 - lj - the local Jacobian 496 497 Level: intermediate 498 499 Notes: The routine SNESDMDAComputeFunction() uses this the cached function to evaluate the user provided function. 500 501 .keywords: distributed array, refine 502 503 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction() 504 @*/ 505 #undef __FUNCT__ 506 #define __FUNCT__ "DMDASetLocalJacobian" 507 PetscErrorCode DMDASetLocalJacobian(DM da,DMDALocalFunction1 lj) 508 { 509 PetscErrorCode ierr; 510 DM_DA *dd = (DM_DA*)da->data; 511 512 PetscFunctionBegin; 513 PetscValidHeaderSpecific(da,DM_CLASSID,1); 514 ierr = DMSetJacobian(da,DMDAJacobian);CHKERRQ(ierr); 515 dd->lj = lj; 516 PetscFunctionReturn(0); 517 } 518 519 #undef __FUNCT__ 520 #define __FUNCT__ "DMDAGetLocalFunction" 521 /*@C 522 DMDAGetLocalFunction - Gets from a DM a local function and its ADIC/ADIFOR Jacobian 523 524 Note Collective 525 526 Input Parameter: 527 . da - initial distributed array 528 529 Output Parameter: 530 . lf - the local function 531 532 Level: intermediate 533 534 .keywords: distributed array, refine 535 536 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalJacobian(), DMDASetLocalFunction() 537 @*/ 538 PetscErrorCode DMDAGetLocalFunction(DM da,DMDALocalFunction1 *lf) 539 { 540 DM_DA *dd = (DM_DA*)da->data; 541 PetscFunctionBegin; 542 PetscValidHeaderSpecific(da,DM_CLASSID,1); 543 if (lf) *lf = dd->lf; 544 PetscFunctionReturn(0); 545 } 546 547 #undef __FUNCT__ 548 #define __FUNCT__ "DMDAGetLocalJacobian" 549 /*@C 550 DMDAGetLocalJacobian - Gets from a DM a local jacobian 551 552 Not Collective 553 554 Input Parameter: 555 . da - initial distributed array 556 557 Output Parameter: 558 . lj - the local jacobian 559 560 Level: intermediate 561 562 .keywords: distributed array, refine 563 564 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalJacobian() 565 @*/ 566 PetscErrorCode DMDAGetLocalJacobian(DM da,DMDALocalFunction1 *lj) 567 { 568 DM_DA *dd = (DM_DA*)da->data; 569 PetscFunctionBegin; 570 PetscValidHeaderSpecific(da,DM_CLASSID,1); 571 if (lj) *lj = dd->lj; 572 PetscFunctionReturn(0); 573 } 574 575 #undef __FUNCT__ 576 #define __FUNCT__ "DMDAComputeFunction" 577 /*@ 578 DMDAComputeFunction - Evaluates a user provided function on each processor that 579 share a DMDA 580 581 Input Parameters: 582 + da - the DM that defines the grid 583 . vu - input vector 584 . vfu - output vector 585 - w - any user data 586 587 Notes: Does NOT do ghost updates on vu upon entry 588 589 This should eventually replace DMDAComputeFunction1 590 591 Level: advanced 592 593 .seealso: DMDAComputeJacobian1WithAdic() 594 595 @*/ 596 PetscErrorCode DMDAComputeFunction(DM da,PetscErrorCode (*lf)(void),Vec vu,Vec vfu,void *w) 597 { 598 PetscErrorCode ierr; 599 void *u,*fu; 600 DMDALocalInfo info; 601 PetscErrorCode (*f)(DMDALocalInfo*,void*,void*,void*) = (PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))lf; 602 603 PetscFunctionBegin; 604 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 605 ierr = DMDAVecGetArray(da,vu,&u);CHKERRQ(ierr); 606 ierr = DMDAVecGetArray(da,vfu,&fu);CHKERRQ(ierr); 607 608 ierr = (*f)(&info,u,fu,w);CHKERRQ(ierr); 609 610 ierr = DMDAVecRestoreArray(da,vu,&u);CHKERRQ(ierr); 611 ierr = DMDAVecRestoreArray(da,vfu,&fu);CHKERRQ(ierr); 612 PetscFunctionReturn(0); 613 } 614 615 #undef __FUNCT__ 616 #define __FUNCT__ "DMDAComputeFunctionLocal" 617 /*@C 618 DMDAComputeFunctionLocal - This is a universal function evaluation routine for 619 a local DM function. 620 621 Collective on DMDA 622 623 Input Parameters: 624 + da - the DM context 625 . func - The local function 626 . X - input vector 627 . F - function vector 628 - ctx - A user context 629 630 Level: intermediate 631 632 .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(), 633 SNESSetFunction(), SNESSetJacobian() 634 635 @*/ 636 PetscErrorCode DMDAComputeFunctionLocal(DM da, DMDALocalFunction1 func, Vec X, Vec F, void *ctx) 637 { 638 Vec localX; 639 DMDALocalInfo info; 640 void *u; 641 void *fu; 642 PetscErrorCode ierr; 643 644 PetscFunctionBegin; 645 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 646 /* 647 Scatter ghost points to local vector, using the 2-step process 648 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 649 */ 650 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 651 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 652 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 653 ierr = DMDAVecGetArray(da,localX,&u);CHKERRQ(ierr); 654 ierr = DMDAVecGetArray(da,F,&fu);CHKERRQ(ierr); 655 ierr = (*func)(&info,u,fu,ctx);CHKERRQ(ierr); 656 ierr = DMDAVecRestoreArray(da,localX,&u);CHKERRQ(ierr); 657 ierr = DMDAVecRestoreArray(da,F,&fu);CHKERRQ(ierr); 658 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 662 #undef __FUNCT__ 663 #define __FUNCT__ "DMDAComputeFunctionLocalGhost" 664 /*@C 665 DMDAComputeFunctionLocalGhost - This is a universal function evaluation routine for 666 a local DM function, but the ghost values of the output are communicated and added. 667 668 Collective on DMDA 669 670 Input Parameters: 671 + da - the DM context 672 . func - The local function 673 . X - input vector 674 . F - function vector 675 - ctx - A user context 676 677 Level: intermediate 678 679 .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(), 680 SNESSetFunction(), SNESSetJacobian() 681 682 @*/ 683 PetscErrorCode DMDAComputeFunctionLocalGhost(DM da, DMDALocalFunction1 func, Vec X, Vec F, void *ctx) 684 { 685 Vec localX, localF; 686 DMDALocalInfo info; 687 void *u; 688 void *fu; 689 PetscErrorCode ierr; 690 691 PetscFunctionBegin; 692 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 693 ierr = DMGetLocalVector(da,&localF);CHKERRQ(ierr); 694 /* 695 Scatter ghost points to local vector, using the 2-step process 696 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 697 */ 698 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 699 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 700 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 701 ierr = VecSet(localF, 0.0);CHKERRQ(ierr); 702 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 703 ierr = DMDAVecGetArray(da,localX,&u);CHKERRQ(ierr); 704 ierr = DMDAVecGetArray(da,localF,&fu);CHKERRQ(ierr); 705 ierr = (*func)(&info,u,fu,ctx);CHKERRQ(ierr); 706 ierr = DMLocalToGlobalBegin(da,localF,ADD_VALUES,F);CHKERRQ(ierr); 707 ierr = DMLocalToGlobalEnd(da,localF,ADD_VALUES,F);CHKERRQ(ierr); 708 ierr = DMDAVecRestoreArray(da,localX,&u);CHKERRQ(ierr); 709 ierr = DMDAVecRestoreArray(da,localF,&fu);CHKERRQ(ierr); 710 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 711 ierr = DMRestoreLocalVector(da,&localF);CHKERRQ(ierr); 712 PetscFunctionReturn(0); 713 } 714 715 #undef __FUNCT__ 716 #define __FUNCT__ "DMDAComputeFunction1" 717 /*@ 718 DMDAComputeFunction1 - Evaluates a user provided function on each processor that 719 share a DMDA 720 721 Input Parameters: 722 + da - the DM that defines the grid 723 . vu - input vector 724 . vfu - output vector 725 - w - any user data 726 727 Notes: Does NOT do ghost updates on vu upon entry 728 729 Level: advanced 730 731 .seealso: DMDAComputeJacobian1WithAdic() 732 733 @*/ 734 PetscErrorCode DMDAComputeFunction1(DM da,Vec vu,Vec vfu,void *w) 735 { 736 PetscErrorCode ierr; 737 void *u,*fu; 738 DMDALocalInfo info; 739 DM_DA *dd = (DM_DA*)da->data; 740 741 PetscFunctionBegin; 742 if (!dd->lf) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_NULL,"DMDASetLocalFunction() never called"); 743 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 744 ierr = DMDAVecGetArray(da,vu,&u);CHKERRQ(ierr); 745 ierr = DMDAVecGetArray(da,vfu,&fu);CHKERRQ(ierr); 746 747 CHKMEMQ; 748 ierr = (*dd->lf)(&info,u,fu,w);CHKERRQ(ierr); 749 CHKMEMQ; 750 751 ierr = DMDAVecRestoreArray(da,vu,&u);CHKERRQ(ierr); 752 ierr = DMDAVecRestoreArray(da,vfu,&fu);CHKERRQ(ierr); 753 PetscFunctionReturn(0); 754 } 755 756 #undef __FUNCT__ 757 #define __FUNCT__ "DMDAComputeFunctioniTest1" 758 PetscErrorCode DMDAComputeFunctioniTest1(DM da,void *w) 759 { 760 Vec vu,fu,fui; 761 PetscErrorCode ierr; 762 PetscInt i,n; 763 PetscScalar *ui; 764 PetscRandom rnd; 765 PetscReal norm; 766 767 PetscFunctionBegin; 768 ierr = DMGetLocalVector(da,&vu);CHKERRQ(ierr); 769 ierr = PetscRandomCreate(PETSC_COMM_SELF,&rnd);CHKERRQ(ierr); 770 ierr = PetscRandomSetFromOptions(rnd);CHKERRQ(ierr); 771 ierr = VecSetRandom(vu,rnd);CHKERRQ(ierr); 772 ierr = PetscRandomDestroy(&rnd);CHKERRQ(ierr); 773 774 ierr = DMGetGlobalVector(da,&fu);CHKERRQ(ierr); 775 ierr = DMGetGlobalVector(da,&fui);CHKERRQ(ierr); 776 777 ierr = DMDAComputeFunction1(da,vu,fu,w);CHKERRQ(ierr); 778 779 ierr = VecGetArray(fui,&ui);CHKERRQ(ierr); 780 ierr = VecGetLocalSize(fui,&n);CHKERRQ(ierr); 781 for (i=0; i<n; i++) { 782 ierr = DMDAComputeFunctioni1(da,i,vu,ui+i,w);CHKERRQ(ierr); 783 } 784 ierr = VecRestoreArray(fui,&ui);CHKERRQ(ierr); 785 786 ierr = VecAXPY(fui,-1.0,fu);CHKERRQ(ierr); 787 ierr = VecNorm(fui,NORM_2,&norm);CHKERRQ(ierr); 788 ierr = PetscPrintf(((PetscObject)da)->comm,"Norm of difference in vectors %G\n",norm);CHKERRQ(ierr); 789 ierr = VecView(fu,0);CHKERRQ(ierr); 790 ierr = VecView(fui,0);CHKERRQ(ierr); 791 792 ierr = DMRestoreLocalVector(da,&vu);CHKERRQ(ierr); 793 ierr = DMRestoreGlobalVector(da,&fu);CHKERRQ(ierr); 794 ierr = DMRestoreGlobalVector(da,&fui);CHKERRQ(ierr); 795 PetscFunctionReturn(0); 796 } 797 798 #undef __FUNCT__ 799 #define __FUNCT__ "DMDAComputeFunctioni1" 800 /*@ 801 DMDAComputeFunctioni1 - Evaluates a user provided point-wise function 802 803 Input Parameters: 804 + da - the DM that defines the grid 805 . i - the component of the function we wish to compute (must be local) 806 . vu - input vector 807 . vfu - output value 808 - w - any user data 809 810 Notes: Does NOT do ghost updates on vu upon entry 811 812 Level: advanced 813 814 .seealso: DMDAComputeJacobian1WithAdic() 815 816 @*/ 817 PetscErrorCode DMDAComputeFunctioni1(DM da,PetscInt i,Vec vu,PetscScalar *vfu,void *w) 818 { 819 PetscErrorCode ierr; 820 void *u; 821 DMDALocalInfo info; 822 MatStencil stencil; 823 DM_DA *dd = (DM_DA*)da->data; 824 825 PetscFunctionBegin; 826 827 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 828 ierr = DMDAVecGetArray(da,vu,&u);CHKERRQ(ierr); 829 830 /* figure out stencil value from i */ 831 stencil.c = i % info.dof; 832 stencil.i = (i % (info.xm*info.dof))/info.dof; 833 stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof); 834 stencil.k = i/(info.xm*info.ym*info.dof); 835 836 ierr = (*dd->lfi)(&info,&stencil,u,vfu,w);CHKERRQ(ierr); 837 838 ierr = DMDAVecRestoreArray(da,vu,&u);CHKERRQ(ierr); 839 PetscFunctionReturn(0); 840 } 841 842 #undef __FUNCT__ 843 #define __FUNCT__ "DMDAComputeFunctionib1" 844 /*@ 845 DMDAComputeFunctionib1 - Evaluates a user provided point-block function 846 847 Input Parameters: 848 + da - the DM that defines the grid 849 . i - the component of the function we wish to compute (must be local) 850 . vu - input vector 851 . vfu - output value 852 - w - any user data 853 854 Notes: Does NOT do ghost updates on vu upon entry 855 856 Level: advanced 857 858 .seealso: DMDAComputeJacobian1WithAdic() 859 860 @*/ 861 PetscErrorCode DMDAComputeFunctionib1(DM da,PetscInt i,Vec vu,PetscScalar *vfu,void *w) 862 { 863 PetscErrorCode ierr; 864 void *u; 865 DMDALocalInfo info; 866 MatStencil stencil; 867 DM_DA *dd = (DM_DA*)da->data; 868 869 PetscFunctionBegin; 870 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 871 ierr = DMDAVecGetArray(da,vu,&u);CHKERRQ(ierr); 872 873 /* figure out stencil value from i */ 874 stencil.c = i % info.dof; 875 if (stencil.c) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Point-block functions can only be called for the entire block"); 876 stencil.i = (i % (info.xm*info.dof))/info.dof; 877 stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof); 878 stencil.k = i/(info.xm*info.ym*info.dof); 879 880 ierr = (*dd->lfib)(&info,&stencil,u,vfu,w);CHKERRQ(ierr); 881 882 ierr = DMDAVecRestoreArray(da,vu,&u);CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 #if defined(new) 887 #undef __FUNCT__ 888 #define __FUNCT__ "DMDAGetDiagonal_MFFD" 889 /* 890 DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local 891 function lives on a DMDA 892 893 y ~= (F(u + ha) - F(u))/h, 894 where F = nonlinear function, as set by SNESSetFunction() 895 u = current iterate 896 h = difference interval 897 */ 898 PetscErrorCode DMDAGetDiagonal_MFFD(DM da,Vec U,Vec a) 899 { 900 PetscScalar h,*aa,*ww,v; 901 PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON; 902 PetscErrorCode ierr; 903 PetscInt gI,nI; 904 MatStencil stencil; 905 DMDALocalInfo info; 906 907 PetscFunctionBegin; 908 ierr = (*ctx->func)(0,U,a,ctx->funcctx);CHKERRQ(ierr); 909 ierr = (*ctx->funcisetbase)(U,ctx->funcctx);CHKERRQ(ierr); 910 911 ierr = VecGetArray(U,&ww);CHKERRQ(ierr); 912 ierr = VecGetArray(a,&aa);CHKERRQ(ierr); 913 914 nI = 0; 915 h = ww[gI]; 916 if (h == 0.0) h = 1.0; 917 #if !defined(PETSC_USE_COMPLEX) 918 if (h < umin && h >= 0.0) h = umin; 919 else if (h < 0.0 && h > -umin) h = -umin; 920 #else 921 if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin; 922 else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin; 923 #endif 924 h *= epsilon; 925 926 ww[gI] += h; 927 ierr = (*ctx->funci)(i,w,&v,ctx->funcctx);CHKERRQ(ierr); 928 aa[nI] = (v - aa[nI])/h; 929 ww[gI] -= h; 930 nI++; 931 } 932 ierr = VecRestoreArray(U,&ww);CHKERRQ(ierr); 933 ierr = VecRestoreArray(a,&aa);CHKERRQ(ierr); 934 PetscFunctionReturn(0); 935 } 936 #endif 937 938 #if defined(PETSC_HAVE_ADIC) 939 EXTERN_C_BEGIN 940 #include <adic/ad_utils.h> 941 EXTERN_C_END 942 943 #undef __FUNCT__ 944 #define __FUNCT__ "DMDAComputeJacobian1WithAdic" 945 /*@C 946 DMDAComputeJacobian1WithAdic - Evaluates a adiC provided Jacobian function on each processor that 947 share a DMDA 948 949 Input Parameters: 950 + da - the DM that defines the grid 951 . vu - input vector (ghosted) 952 . J - output matrix 953 - w - any user data 954 955 Level: advanced 956 957 Notes: Does NOT do ghost updates on vu upon entry 958 959 .seealso: DMDAComputeFunction1() 960 961 @*/ 962 PetscErrorCode DMDAComputeJacobian1WithAdic(DM da,Vec vu,Mat J,void *w) 963 { 964 PetscErrorCode ierr; 965 PetscInt gtdof,tdof; 966 PetscScalar *ustart; 967 DMDALocalInfo info; 968 void *ad_u,*ad_f,*ad_ustart,*ad_fstart; 969 ISColoring iscoloring; 970 971 PetscFunctionBegin; 972 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 973 974 PetscADResetIndep(); 975 976 /* get space for derivative objects. */ 977 ierr = DMDAGetAdicArray(da,PETSC_TRUE,&ad_u,&ad_ustart,>dof);CHKERRQ(ierr); 978 ierr = DMDAGetAdicArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);CHKERRQ(ierr); 979 ierr = VecGetArray(vu,&ustart);CHKERRQ(ierr); 980 ierr = DMCreateColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);CHKERRQ(ierr); 981 982 PetscADSetValueAndColor(ad_ustart,gtdof,iscoloring->colors,ustart); 983 984 ierr = VecRestoreArray(vu,&ustart);CHKERRQ(ierr); 985 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 986 ierr = PetscADIncrementTotalGradSize(iscoloring->n);CHKERRQ(ierr); 987 PetscADSetIndepDone(); 988 989 ierr = PetscLogEventBegin(DMDA_LocalADFunction,0,0,0,0);CHKERRQ(ierr); 990 ierr = (*dd->adic_lf)(&info,ad_u,ad_f,w);CHKERRQ(ierr); 991 ierr = PetscLogEventEnd(DMDA_LocalADFunction,0,0,0,0);CHKERRQ(ierr); 992 993 /* stick the values into the matrix */ 994 ierr = MatSetValuesAdic(J,(PetscScalar**)ad_fstart);CHKERRQ(ierr); 995 996 /* return space for derivative objects. */ 997 ierr = DMDARestoreAdicArray(da,PETSC_TRUE,&ad_u,&ad_ustart,>dof);CHKERRQ(ierr); 998 ierr = DMDARestoreAdicArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);CHKERRQ(ierr); 999 PetscFunctionReturn(0); 1000 } 1001 1002 #undef __FUNCT__ 1003 #define __FUNCT__ "DMDAMultiplyByJacobian1WithAdic" 1004 /*@C 1005 DMDAMultiplyByJacobian1WithAdic - Applies an ADIC-provided Jacobian function to a vector on 1006 each processor that shares a DMDA. 1007 1008 Input Parameters: 1009 + da - the DM that defines the grid 1010 . vu - Jacobian is computed at this point (ghosted) 1011 . v - product is done on this vector (ghosted) 1012 . fu - output vector = J(vu)*v (not ghosted) 1013 - w - any user data 1014 1015 Notes: 1016 This routine does NOT do ghost updates on vu upon entry. 1017 1018 Level: advanced 1019 1020 .seealso: DMDAComputeFunction1() 1021 1022 @*/ 1023 PetscErrorCode DMDAMultiplyByJacobian1WithAdic(DM da,Vec vu,Vec v,Vec f,void *w) 1024 { 1025 PetscErrorCode ierr; 1026 PetscInt i,gtdof,tdof; 1027 PetscScalar *avu,*av,*af,*ad_vustart,*ad_fstart; 1028 DMDALocalInfo info; 1029 void *ad_vu,*ad_f; 1030 1031 PetscFunctionBegin; 1032 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 1033 1034 /* get space for derivative objects. */ 1035 ierr = DMDAGetAdicMFArray(da,PETSC_TRUE,&ad_vu,&ad_vustart,>dof);CHKERRQ(ierr); 1036 ierr = DMDAGetAdicMFArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);CHKERRQ(ierr); 1037 1038 /* copy input vector into derivative object */ 1039 ierr = VecGetArray(vu,&avu);CHKERRQ(ierr); 1040 ierr = VecGetArray(v,&av);CHKERRQ(ierr); 1041 for (i=0; i<gtdof; i++) { 1042 ad_vustart[2*i] = avu[i]; 1043 ad_vustart[2*i+1] = av[i]; 1044 } 1045 ierr = VecRestoreArray(vu,&avu);CHKERRQ(ierr); 1046 ierr = VecRestoreArray(v,&av);CHKERRQ(ierr); 1047 1048 PetscADResetIndep(); 1049 ierr = PetscADIncrementTotalGradSize(1);CHKERRQ(ierr); 1050 PetscADSetIndepDone(); 1051 1052 ierr = (*dd->adicmf_lf)(&info,ad_vu,ad_f,w);CHKERRQ(ierr); 1053 1054 /* stick the values into the vector */ 1055 ierr = VecGetArray(f,&af);CHKERRQ(ierr); 1056 for (i=0; i<tdof; i++) { 1057 af[i] = ad_fstart[2*i+1]; 1058 } 1059 ierr = VecRestoreArray(f,&af);CHKERRQ(ierr); 1060 1061 /* return space for derivative objects. */ 1062 ierr = DMDARestoreAdicMFArray(da,PETSC_TRUE,&ad_vu,&ad_vustart,>dof);CHKERRQ(ierr); 1063 ierr = DMDARestoreAdicMFArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);CHKERRQ(ierr); 1064 PetscFunctionReturn(0); 1065 } 1066 #endif 1067 1068 #undef __FUNCT__ 1069 #define __FUNCT__ "DMDAComputeJacobian1" 1070 /*@ 1071 DMDAComputeJacobian1 - Evaluates a local Jacobian function on each processor that 1072 share a DMDA 1073 1074 Input Parameters: 1075 + da - the DM that defines the grid 1076 . vu - input vector (ghosted) 1077 . J - output matrix 1078 - w - any user data 1079 1080 Notes: Does NOT do ghost updates on vu upon entry 1081 1082 Level: advanced 1083 1084 .seealso: DMDAComputeFunction1() 1085 1086 @*/ 1087 PetscErrorCode DMDAComputeJacobian1(DM da,Vec vu,Mat J,void *w) 1088 { 1089 PetscErrorCode ierr; 1090 void *u; 1091 DMDALocalInfo info; 1092 DM_DA *dd = (DM_DA*)da->data; 1093 1094 PetscFunctionBegin; 1095 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 1096 ierr = DMDAVecGetArray(da,vu,&u);CHKERRQ(ierr); 1097 ierr = (*dd->lj)(&info,u,J,w);CHKERRQ(ierr); 1098 ierr = DMDAVecRestoreArray(da,vu,&u);CHKERRQ(ierr); 1099 PetscFunctionReturn(0); 1100 } 1101 1102 1103 #undef __FUNCT__ 1104 #define __FUNCT__ "DMDAComputeJacobian1WithAdifor" 1105 /* 1106 DMDAComputeJacobian1WithAdifor - Evaluates a ADIFOR provided Jacobian local function on each processor that 1107 share a DMDA 1108 1109 Input Parameters: 1110 + da - the DM that defines the grid 1111 . vu - input vector (ghosted) 1112 . J - output matrix 1113 - w - any user data 1114 1115 Notes: Does NOT do ghost updates on vu upon entry 1116 1117 .seealso: DMDAComputeFunction1() 1118 1119 */ 1120 PetscErrorCode DMDAComputeJacobian1WithAdifor(DM da,Vec vu,Mat J,void *w) 1121 { 1122 PetscErrorCode ierr; 1123 PetscInt i,Nc,N; 1124 ISColoringValue *color; 1125 DMDALocalInfo info; 1126 PetscScalar *u,*g_u,*g_f,*f = 0,*p_u; 1127 ISColoring iscoloring; 1128 DM_DA *dd = (DM_DA*)da->data; 1129 void (*lf)(PetscInt*,DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*) = 1130 (void (*)(PetscInt*,DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*))*dd->adifor_lf; 1131 1132 PetscFunctionBegin; 1133 ierr = DMCreateColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);CHKERRQ(ierr); 1134 Nc = iscoloring->n; 1135 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 1136 N = info.gxm*info.gym*info.gzm*info.dof; 1137 1138 /* get space for derivative objects. */ 1139 ierr = PetscMalloc(Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar),&g_u);CHKERRQ(ierr); 1140 ierr = PetscMemzero(g_u,Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar));CHKERRQ(ierr); 1141 p_u = g_u; 1142 color = iscoloring->colors; 1143 for (i=0; i<N; i++) { 1144 p_u[*color++] = 1.0; 1145 p_u += Nc; 1146 } 1147 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 1148 ierr = PetscMalloc2(Nc*info.xm*info.ym*info.zm*info.dof,PetscScalar,&g_f,info.xm*info.ym*info.zm*info.dof,PetscScalar,&f);CHKERRQ(ierr); 1149 1150 /* Seed the input array g_u with coloring information */ 1151 1152 ierr = VecGetArray(vu,&u);CHKERRQ(ierr); 1153 (lf)(&Nc,&info,u,g_u,&Nc,f,g_f,&Nc,w,&ierr);CHKERRQ(ierr); 1154 ierr = VecRestoreArray(vu,&u);CHKERRQ(ierr); 1155 1156 /* stick the values into the matrix */ 1157 /* PetscScalarView(Nc*info.xm*info.ym,g_f,0); */ 1158 ierr = MatSetValuesAdifor(J,Nc,g_f);CHKERRQ(ierr); 1159 1160 /* return space for derivative objects. */ 1161 ierr = PetscFree(g_u);CHKERRQ(ierr); 1162 ierr = PetscFree2(g_f,f);CHKERRQ(ierr); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "DMDAFormJacobianLocal" 1168 /*@C 1169 DMDAFormjacobianLocal - This is a universal Jacobian evaluation routine for 1170 a local DM function. 1171 1172 Collective on DMDA 1173 1174 Input Parameters: 1175 + da - the DM context 1176 . func - The local function 1177 . X - input vector 1178 . J - Jacobian matrix 1179 - ctx - A user context 1180 1181 Level: intermediate 1182 1183 .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(), 1184 SNESSetFunction(), SNESSetJacobian() 1185 1186 @*/ 1187 PetscErrorCode DMDAFormJacobianLocal(DM da, DMDALocalFunction1 func, Vec X, Mat J, void *ctx) 1188 { 1189 Vec localX; 1190 DMDALocalInfo info; 1191 void *u; 1192 PetscErrorCode ierr; 1193 1194 PetscFunctionBegin; 1195 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 1196 /* 1197 Scatter ghost points to local vector, using the 2-step process 1198 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 1199 */ 1200 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 1201 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 1202 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 1203 ierr = DMDAVecGetArray(da,localX,&u);CHKERRQ(ierr); 1204 ierr = (*func)(&info,u,J,ctx);CHKERRQ(ierr); 1205 ierr = DMDAVecRestoreArray(da,localX,&u);CHKERRQ(ierr); 1206 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 1207 PetscFunctionReturn(0); 1208 } 1209 1210 #undef __FUNCT__ 1211 #define __FUNCT__ "DMDAMultiplyByJacobian1WithAD" 1212 /*@C 1213 DMDAMultiplyByJacobian1WithAD - Applies a Jacobian function supplied by ADIFOR or ADIC 1214 to a vector on each processor that shares a DMDA. 1215 1216 Input Parameters: 1217 + da - the DM that defines the grid 1218 . vu - Jacobian is computed at this point (ghosted) 1219 . v - product is done on this vector (ghosted) 1220 . fu - output vector = J(vu)*v (not ghosted) 1221 - w - any user data 1222 1223 Notes: 1224 This routine does NOT do ghost updates on vu and v upon entry. 1225 1226 Automatically calls DMDAMultiplyByJacobian1WithAdifor() or DMDAMultiplyByJacobian1WithAdic() 1227 depending on whether DMDASetLocalAdicMFFunction() or DMDASetLocalAdiforMFFunction() was called. 1228 1229 Level: advanced 1230 1231 .seealso: DMDAComputeFunction1(), DMDAMultiplyByJacobian1WithAdifor(), DMDAMultiplyByJacobian1WithAdic() 1232 1233 @*/ 1234 PetscErrorCode DMDAMultiplyByJacobian1WithAD(DM da,Vec u,Vec v,Vec f,void *w) 1235 { 1236 PetscErrorCode ierr; 1237 DM_DA *dd = (DM_DA*)da->data; 1238 1239 PetscFunctionBegin; 1240 if (dd->adicmf_lf) { 1241 #if defined(PETSC_HAVE_ADIC) 1242 ierr = DMDAMultiplyByJacobian1WithAdic(da,u,v,f,w);CHKERRQ(ierr); 1243 #else 1244 SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP_SYS,"Requires ADIC to be installed and cannot use complex numbers"); 1245 #endif 1246 } else if (dd->adiformf_lf) { 1247 ierr = DMDAMultiplyByJacobian1WithAdifor(da,u,v,f,w);CHKERRQ(ierr); 1248 } else SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ORDER,"Must call DMDASetLocalAdiforMFFunction() or DMDASetLocalAdicMFFunction() before using"); 1249 PetscFunctionReturn(0); 1250 } 1251 1252 1253 #undef __FUNCT__ 1254 #define __FUNCT__ "DMDAMultiplyByJacobian1WithAdifor" 1255 /*@C 1256 DMDAMultiplyByJacobian1WithAdifor - Applies a ADIFOR provided Jacobian function on each processor that 1257 share a DM to a vector 1258 1259 Input Parameters: 1260 + da - the DM that defines the grid 1261 . vu - Jacobian is computed at this point (ghosted) 1262 . v - product is done on this vector (ghosted) 1263 . fu - output vector = J(vu)*v (not ghosted) 1264 - w - any user data 1265 1266 Notes: Does NOT do ghost updates on vu and v upon entry 1267 1268 Level: advanced 1269 1270 .seealso: DMDAComputeFunction1() 1271 1272 @*/ 1273 PetscErrorCode DMDAMultiplyByJacobian1WithAdifor(DM da,Vec u,Vec v,Vec f,void *w) 1274 { 1275 PetscErrorCode ierr; 1276 PetscScalar *au,*av,*af,*awork; 1277 Vec work; 1278 DMDALocalInfo info; 1279 DM_DA *dd = (DM_DA*)da->data; 1280 void (*lf)(DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*) = 1281 (void (*)(DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*))*dd->adiformf_lf; 1282 1283 PetscFunctionBegin; 1284 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 1285 1286 ierr = DMGetGlobalVector(da,&work);CHKERRQ(ierr); 1287 ierr = VecGetArray(u,&au);CHKERRQ(ierr); 1288 ierr = VecGetArray(v,&av);CHKERRQ(ierr); 1289 ierr = VecGetArray(f,&af);CHKERRQ(ierr); 1290 ierr = VecGetArray(work,&awork);CHKERRQ(ierr); 1291 (lf)(&info,au,av,awork,af,w,&ierr);CHKERRQ(ierr); 1292 ierr = VecRestoreArray(u,&au);CHKERRQ(ierr); 1293 ierr = VecRestoreArray(v,&av);CHKERRQ(ierr); 1294 ierr = VecRestoreArray(f,&af);CHKERRQ(ierr); 1295 ierr = VecRestoreArray(work,&awork);CHKERRQ(ierr); 1296 ierr = DMRestoreGlobalVector(da,&work);CHKERRQ(ierr); 1297 1298 PetscFunctionReturn(0); 1299 } 1300 1301 #undef __FUNCT__ 1302 #define __FUNCT__ "DMSetUp_DA_2D" 1303 PetscErrorCode DMSetUp_DA_2D(DM da) 1304 { 1305 DM_DA *dd = (DM_DA*)da->data; 1306 const PetscInt M = dd->M; 1307 const PetscInt N = dd->N; 1308 PetscInt m = dd->m; 1309 PetscInt n = dd->n; 1310 const PetscInt dof = dd->w; 1311 const PetscInt s = dd->s; 1312 const DMDABoundaryType bx = dd->bx; 1313 const DMDABoundaryType by = dd->by; 1314 const DMDAStencilType stencil_type = dd->stencil_type; 1315 PetscInt *lx = dd->lx; 1316 PetscInt *ly = dd->ly; 1317 MPI_Comm comm; 1318 PetscMPIInt rank,size; 1319 PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end,IXs,IXe,IYs,IYe; 1320 PetscInt up,down,left,right,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn,*idx_cpy; 1321 const PetscInt *idx_full; 1322 PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count; 1323 PetscInt s_x,s_y; /* s proportionalized to w */ 1324 PetscInt sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0; 1325 Vec local,global; 1326 VecScatter ltog,gtol; 1327 IS to,from,ltogis; 1328 PetscErrorCode ierr; 1329 1330 PetscFunctionBegin; 1331 if (dof < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof); 1332 if (s < 0) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s); 1333 ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 1334 #if !defined(PETSC_USE_64BIT_INDICES) 1335 if (((Petsc64bitInt) M)*((Petsc64bitInt) N)*((Petsc64bitInt) dof) > (Petsc64bitInt) PETSC_MPI_INT_MAX) SETERRQ3(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D (dof) is too large for 32 bit indices",M,N,dof); 1336 #endif 1337 1338 if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof); 1339 if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s); 1340 1341 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1342 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1343 1344 dd->dim = 2; 1345 ierr = PetscMalloc(dof*sizeof(char*),&dd->fieldname);CHKERRQ(ierr); 1346 ierr = PetscMemzero(dd->fieldname,dof*sizeof(char*));CHKERRQ(ierr); 1347 1348 if (m != PETSC_DECIDE) { 1349 if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m); 1350 else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size); 1351 } 1352 if (n != PETSC_DECIDE) { 1353 if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n); 1354 else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size); 1355 } 1356 1357 if (m == PETSC_DECIDE || n == PETSC_DECIDE) { 1358 if (n != PETSC_DECIDE) { 1359 m = size/n; 1360 } else if (m != PETSC_DECIDE) { 1361 n = size/m; 1362 } else { 1363 /* try for squarish distribution */ 1364 m = (PetscInt)(0.5 + sqrt(((double)M)*((double)size)/((double)N))); 1365 if (!m) m = 1; 1366 while (m > 0) { 1367 n = size/m; 1368 if (m*n == size) break; 1369 m--; 1370 } 1371 if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;} 1372 } 1373 if (m*n != size) SETERRQ(comm,PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n "); 1374 } else if (m*n != size) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition"); 1375 1376 if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m); 1377 if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n); 1378 1379 /* 1380 Determine locally owned region 1381 xs is the first local node number, x is the number of local nodes 1382 */ 1383 if (!lx) { 1384 ierr = PetscMalloc(m*sizeof(PetscInt), &dd->lx);CHKERRQ(ierr); 1385 lx = dd->lx; 1386 for (i=0; i<m; i++) { 1387 lx[i] = M/m + ((M % m) > i); 1388 } 1389 } 1390 x = lx[rank % m]; 1391 xs = 0; 1392 for (i=0; i<(rank % m); i++) { 1393 xs += lx[i]; 1394 } 1395 #if defined(PETSC_USE_DEBUG) 1396 left = xs; 1397 for (i=(rank % m); i<m; i++) { 1398 left += lx[i]; 1399 } 1400 if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M); 1401 #endif 1402 1403 /* 1404 Determine locally owned region 1405 ys is the first local node number, y is the number of local nodes 1406 */ 1407 if (!ly) { 1408 ierr = PetscMalloc(n*sizeof(PetscInt), &dd->ly);CHKERRQ(ierr); 1409 ly = dd->ly; 1410 for (i=0; i<n; i++) { 1411 ly[i] = N/n + ((N % n) > i); 1412 } 1413 } 1414 y = ly[rank/m]; 1415 ys = 0; 1416 for (i=0; i<(rank/m); i++) { 1417 ys += ly[i]; 1418 } 1419 #if defined(PETSC_USE_DEBUG) 1420 left = ys; 1421 for (i=(rank/m); i<n; i++) { 1422 left += ly[i]; 1423 } 1424 if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N); 1425 #endif 1426 1427 /* 1428 check if the scatter requires more than one process neighbor or wraps around 1429 the domain more than once 1430 */ 1431 if ((x < s) && ((m > 1) || (bx == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s); 1432 if ((y < s) && ((n > 1) || (by == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s); 1433 xe = xs + x; 1434 ye = ys + y; 1435 1436 /* determine ghost region (Xs) and region scattered into (IXs) */ 1437 /* Assume No Periodicity */ 1438 if (xs-s > 0) { Xs = xs - s; IXs = xs - s; } else { Xs = 0; IXs = 0; } 1439 if (xe+s <= M) { Xe = xe + s; IXe = xe + s; } else { Xe = M; IXe = M; } 1440 if (ys-s > 0) { Ys = ys - s; IYs = ys - s; } else { Ys = 0; IYs = 0; } 1441 if (ye+s <= N) { Ye = ye + s; IYe = ye + s; } else { Ye = N; IYe = N; } 1442 1443 /* fix for periodicity/ghosted */ 1444 if (bx) { Xs = xs - s; Xe = xe + s; } 1445 if (bx == DMDA_BOUNDARY_PERIODIC) { IXs = xs - s; IXe = xe + s; } 1446 if (by) { Ys = ys - s; Ye = ye + s; } 1447 if (by == DMDA_BOUNDARY_PERIODIC) { IYs = ys - s; IYe = ye + s; } 1448 1449 /* Resize all X parameters to reflect w */ 1450 s_x = s; 1451 s_y = s; 1452 1453 /* determine starting point of each processor */ 1454 nn = x*y; 1455 ierr = PetscMalloc2(size+1,PetscInt,&bases,size,PetscInt,&ldims);CHKERRQ(ierr); 1456 ierr = MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);CHKERRQ(ierr); 1457 bases[0] = 0; 1458 for (i=1; i<=size; i++) { 1459 bases[i] = ldims[i-1]; 1460 } 1461 for (i=1; i<=size; i++) { 1462 bases[i] += bases[i-1]; 1463 } 1464 base = bases[rank]*dof; 1465 1466 /* allocate the base parallel and sequential vectors */ 1467 dd->Nlocal = x*y*dof; 1468 ierr = VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);CHKERRQ(ierr); 1469 dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof; 1470 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);CHKERRQ(ierr); 1471 1472 /* generate appropriate vector scatters */ 1473 /* local to global inserts non-ghost point region into global */ 1474 ierr = VecGetOwnershipRange(global,&start,&end);CHKERRQ(ierr); 1475 ierr = ISCreateStride(comm,x*y*dof,start,1,&to);CHKERRQ(ierr); 1476 1477 count = x*y; 1478 ierr = PetscMalloc(x*y*sizeof(PetscInt),&idx);CHKERRQ(ierr); 1479 left = xs - Xs; right = left + x; 1480 down = ys - Ys; up = down + y; 1481 count = 0; 1482 for (i=down; i<up; i++) { 1483 for (j=left; j<right; j++) { 1484 idx[count++] = i*(Xe-Xs) + j; 1485 } 1486 } 1487 1488 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);CHKERRQ(ierr); 1489 ierr = VecScatterCreate(local,from,global,to,<og);CHKERRQ(ierr); 1490 ierr = PetscLogObjectParent(dd,ltog);CHKERRQ(ierr); 1491 ierr = ISDestroy(&from);CHKERRQ(ierr); 1492 ierr = ISDestroy(&to);CHKERRQ(ierr); 1493 1494 /* global to local must include ghost points within the domain, 1495 but not ghost points outside the domain that aren't periodic */ 1496 if (stencil_type == DMDA_STENCIL_BOX) { 1497 count = (IXe-IXs)*(IYe-IYs); 1498 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 1499 1500 left = IXs - Xs; right = left + (IXe-IXs); 1501 down = IYs - Ys; up = down + (IYe-IYs); 1502 count = 0; 1503 for (i=down; i<up; i++) { 1504 for (j=left; j<right; j++) { 1505 idx[count++] = j + i*(Xe-Xs); 1506 } 1507 } 1508 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 1509 1510 } else { 1511 /* must drop into cross shape region */ 1512 /* ---------| 1513 | top | 1514 |--- ---| up 1515 | middle | 1516 | | 1517 ---- ---- down 1518 | bottom | 1519 ----------- 1520 Xs xs xe Xe */ 1521 count = (ys-IYs)*x + y*(IXe-IXs) + (IYe-ye)*x; 1522 ierr = PetscMalloc(count*sizeof(PetscInt),&idx);CHKERRQ(ierr); 1523 1524 left = xs - Xs; right = left + x; 1525 down = ys - Ys; up = down + y; 1526 count = 0; 1527 /* bottom */ 1528 for (i=(IYs-Ys); i<down; i++) { 1529 for (j=left; j<right; j++) { 1530 idx[count++] = j + i*(Xe-Xs); 1531 } 1532 } 1533 /* middle */ 1534 for (i=down; i<up; i++) { 1535 for (j=(IXs-Xs); j<(IXe-Xs); j++) { 1536 idx[count++] = j + i*(Xe-Xs); 1537 } 1538 } 1539 /* top */ 1540 for (i=up; i<up+IYe-ye; i++) { 1541 for (j=left; j<right; j++) { 1542 idx[count++] = j + i*(Xe-Xs); 1543 } 1544 } 1545 ierr = ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);CHKERRQ(ierr); 1546 } 1547 1548 1549 /* determine who lies on each side of us stored in n6 n7 n8 1550 n3 n5 1551 n0 n1 n2 1552 */ 1553 1554 /* Assume the Non-Periodic Case */ 1555 n1 = rank - m; 1556 if (rank % m) { 1557 n0 = n1 - 1; 1558 } else { 1559 n0 = -1; 1560 } 1561 if ((rank+1) % m) { 1562 n2 = n1 + 1; 1563 n5 = rank + 1; 1564 n8 = rank + m + 1; if (n8 >= m*n) n8 = -1; 1565 } else { 1566 n2 = -1; n5 = -1; n8 = -1; 1567 } 1568 if (rank % m) { 1569 n3 = rank - 1; 1570 n6 = n3 + m; if (n6 >= m*n) n6 = -1; 1571 } else { 1572 n3 = -1; n6 = -1; 1573 } 1574 n7 = rank + m; if (n7 >= m*n) n7 = -1; 1575 1576 if (bx == DMDA_BOUNDARY_PERIODIC && by == DMDA_BOUNDARY_PERIODIC) { 1577 /* Modify for Periodic Cases */ 1578 /* Handle all four corners */ 1579 if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1; 1580 if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0; 1581 if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m; 1582 if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1; 1583 1584 /* Handle Top and Bottom Sides */ 1585 if (n1 < 0) n1 = rank + m * (n-1); 1586 if (n7 < 0) n7 = rank - m * (n-1); 1587 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 1588 if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1; 1589 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 1590 if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1; 1591 1592 /* Handle Left and Right Sides */ 1593 if (n3 < 0) n3 = rank + (m-1); 1594 if (n5 < 0) n5 = rank - (m-1); 1595 if ((n1 >= 0) && (n0 < 0)) n0 = rank-1; 1596 if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1; 1597 if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1; 1598 if ((n7 >= 0) && (n8 < 0)) n8 = rank+1; 1599 } else if (by == DMDA_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */ 1600 if (n1 < 0) n1 = rank + m * (n-1); 1601 if (n7 < 0) n7 = rank - m * (n-1); 1602 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 1603 if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1; 1604 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 1605 if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1; 1606 } else if (bx == DMDA_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */ 1607 if (n3 < 0) n3 = rank + (m-1); 1608 if (n5 < 0) n5 = rank - (m-1); 1609 if ((n1 >= 0) && (n0 < 0)) n0 = rank-1; 1610 if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1; 1611 if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1; 1612 if ((n7 >= 0) && (n8 < 0)) n8 = rank+1; 1613 } 1614 1615 ierr = PetscMalloc(9*sizeof(PetscInt),&dd->neighbors);CHKERRQ(ierr); 1616 dd->neighbors[0] = n0; 1617 dd->neighbors[1] = n1; 1618 dd->neighbors[2] = n2; 1619 dd->neighbors[3] = n3; 1620 dd->neighbors[4] = rank; 1621 dd->neighbors[5] = n5; 1622 dd->neighbors[6] = n6; 1623 dd->neighbors[7] = n7; 1624 dd->neighbors[8] = n8; 1625 1626 if (stencil_type == DMDA_STENCIL_STAR) { 1627 /* save corner processor numbers */ 1628 sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8; 1629 n0 = n2 = n6 = n8 = -1; 1630 } 1631 1632 ierr = PetscMalloc((Xe-Xs)*(Ye-Ys)*sizeof(PetscInt),&idx);CHKERRQ(ierr); 1633 ierr = PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*sizeof(PetscInt));CHKERRQ(ierr); 1634 1635 nn = 0; 1636 xbase = bases[rank]; 1637 for (i=1; i<=s_y; i++) { 1638 if (n0 >= 0) { /* left below */ 1639 x_t = lx[n0 % m]; 1640 y_t = ly[(n0/m)]; 1641 s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x; 1642 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1643 } 1644 if (n1 >= 0) { /* directly below */ 1645 x_t = x; 1646 y_t = ly[(n1/m)]; 1647 s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t; 1648 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1649 } 1650 if (n2 >= 0) { /* right below */ 1651 x_t = lx[n2 % m]; 1652 y_t = ly[(n2/m)]; 1653 s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t; 1654 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1655 } 1656 } 1657 1658 for (i=0; i<y; i++) { 1659 if (n3 >= 0) { /* directly left */ 1660 x_t = lx[n3 % m]; 1661 /* y_t = y; */ 1662 s_t = bases[n3] + (i+1)*x_t - s_x; 1663 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1664 } 1665 1666 for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */ 1667 1668 if (n5 >= 0) { /* directly right */ 1669 x_t = lx[n5 % m]; 1670 /* y_t = y; */ 1671 s_t = bases[n5] + (i)*x_t; 1672 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1673 } 1674 } 1675 1676 for (i=1; i<=s_y; i++) { 1677 if (n6 >= 0) { /* left above */ 1678 x_t = lx[n6 % m]; 1679 /* y_t = ly[(n6/m)]; */ 1680 s_t = bases[n6] + (i)*x_t - s_x; 1681 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1682 } 1683 if (n7 >= 0) { /* directly above */ 1684 x_t = x; 1685 /* y_t = ly[(n7/m)]; */ 1686 s_t = bases[n7] + (i-1)*x_t; 1687 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1688 } 1689 if (n8 >= 0) { /* right above */ 1690 x_t = lx[n8 % m]; 1691 /* y_t = ly[(n8/m)]; */ 1692 s_t = bases[n8] + (i-1)*x_t; 1693 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1694 } 1695 } 1696 1697 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 1698 ierr = VecScatterCreate(global,from,local,to,>ol);CHKERRQ(ierr); 1699 ierr = PetscLogObjectParent(da,gtol);CHKERRQ(ierr); 1700 ierr = ISDestroy(&to);CHKERRQ(ierr); 1701 ierr = ISDestroy(&from);CHKERRQ(ierr); 1702 1703 if (stencil_type == DMDA_STENCIL_STAR) { 1704 n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8; 1705 } 1706 1707 if ((stencil_type == DMDA_STENCIL_STAR) || 1708 (bx && bx != DMDA_BOUNDARY_PERIODIC) || 1709 (by && by != DMDA_BOUNDARY_PERIODIC)) { 1710 /* 1711 Recompute the local to global mappings, this time keeping the 1712 information about the cross corner processor numbers and any ghosted 1713 but not periodic indices. 1714 */ 1715 nn = 0; 1716 xbase = bases[rank]; 1717 for (i=1; i<=s_y; i++) { 1718 if (n0 >= 0) { /* left below */ 1719 x_t = lx[n0 % m]; 1720 y_t = ly[(n0/m)]; 1721 s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x; 1722 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1723 } else if (xs-Xs > 0 && ys-Ys > 0) { 1724 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1725 } 1726 if (n1 >= 0) { /* directly below */ 1727 x_t = x; 1728 y_t = ly[(n1/m)]; 1729 s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t; 1730 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1731 } else if (ys-Ys > 0) { 1732 for (j=0; j<x; j++) { idx[nn++] = -1;} 1733 } 1734 if (n2 >= 0) { /* right below */ 1735 x_t = lx[n2 % m]; 1736 y_t = ly[(n2/m)]; 1737 s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t; 1738 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1739 } else if (Xe-xe> 0 && ys-Ys > 0) { 1740 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1741 } 1742 } 1743 1744 for (i=0; i<y; i++) { 1745 if (n3 >= 0) { /* directly left */ 1746 x_t = lx[n3 % m]; 1747 /* y_t = y; */ 1748 s_t = bases[n3] + (i+1)*x_t - s_x; 1749 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1750 } else if (xs-Xs > 0) { 1751 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1752 } 1753 1754 for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */ 1755 1756 if (n5 >= 0) { /* directly right */ 1757 x_t = lx[n5 % m]; 1758 /* y_t = y; */ 1759 s_t = bases[n5] + (i)*x_t; 1760 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1761 } else if (Xe-xe > 0) { 1762 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1763 } 1764 } 1765 1766 for (i=1; i<=s_y; i++) { 1767 if (n6 >= 0) { /* left above */ 1768 x_t = lx[n6 % m]; 1769 /* y_t = ly[(n6/m)]; */ 1770 s_t = bases[n6] + (i)*x_t - s_x; 1771 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1772 } else if (xs-Xs > 0 && Ye-ye > 0) { 1773 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1774 } 1775 if (n7 >= 0) { /* directly above */ 1776 x_t = x; 1777 /* y_t = ly[(n7/m)]; */ 1778 s_t = bases[n7] + (i-1)*x_t; 1779 for (j=0; j<x_t; j++) { idx[nn++] = s_t++;} 1780 } else if (Ye-ye > 0) { 1781 for (j=0; j<x; j++) { idx[nn++] = -1;} 1782 } 1783 if (n8 >= 0) { /* right above */ 1784 x_t = lx[n8 % m]; 1785 /* y_t = ly[(n8/m)]; */ 1786 s_t = bases[n8] + (i-1)*x_t; 1787 for (j=0; j<s_x; j++) { idx[nn++] = s_t++;} 1788 } else if (Xe-xe > 0 && Ye-ye > 0) { 1789 for (j=0; j<s_x; j++) { idx[nn++] = -1;} 1790 } 1791 } 1792 } 1793 /* 1794 Set the local to global ordering in the global vector, this allows use 1795 of VecSetValuesLocal(). 1796 */ 1797 ierr = ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);CHKERRQ(ierr); 1798 ierr = PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);CHKERRQ(ierr); 1799 ierr = PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1800 ierr = ISGetIndices(ltogis, &idx_full); 1801 ierr = PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));CHKERRQ(ierr); 1802 ierr = ISRestoreIndices(ltogis, &idx_full); 1803 ierr = ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);CHKERRQ(ierr); 1804 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1805 ierr = ISDestroy(<ogis);CHKERRQ(ierr); 1806 ierr = ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);CHKERRQ(ierr); 1807 ierr = PetscLogObjectParent(da,da->ltogmap);CHKERRQ(ierr); 1808 1809 ierr = PetscFree2(bases,ldims);CHKERRQ(ierr); 1810 dd->m = m; dd->n = n; 1811 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 1812 dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1; 1813 dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1; 1814 1815 ierr = VecDestroy(&local);CHKERRQ(ierr); 1816 ierr = VecDestroy(&global);CHKERRQ(ierr); 1817 1818 dd->gtol = gtol; 1819 dd->ltog = ltog; 1820 dd->idx = idx_cpy; 1821 dd->Nl = nn*dof; 1822 dd->base = base; 1823 da->ops->view = DMView_DA_2d; 1824 dd->ltol = PETSC_NULL; 1825 dd->ao = PETSC_NULL; 1826 1827 PetscFunctionReturn(0); 1828 } 1829 1830 #undef __FUNCT__ 1831 #define __FUNCT__ "DMDACreate2d" 1832 /*@C 1833 DMDACreate2d - Creates an object that will manage the communication of two-dimensional 1834 regular array data that is distributed across some processors. 1835 1836 Collective on MPI_Comm 1837 1838 Input Parameters: 1839 + comm - MPI communicator 1840 . bx,by - type of ghost nodes the array have. 1841 Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC. 1842 . stencil_type - stencil type. Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR. 1843 . M,N - global dimension in each direction of the array (use -M and or -N to indicate that it may be set to a different value 1844 from the command line with -da_grid_x <M> -da_grid_y <N>) 1845 . m,n - corresponding number of processors in each dimension 1846 (or PETSC_DECIDE to have calculated) 1847 . dof - number of degrees of freedom per node 1848 . s - stencil width 1849 - lx, ly - arrays containing the number of nodes in each cell along 1850 the x and y coordinates, or PETSC_NULL. If non-null, these 1851 must be of length as m and n, and the corresponding 1852 m and n cannot be PETSC_DECIDE. The sum of the lx[] entries 1853 must be M, and the sum of the ly[] entries must be N. 1854 1855 Output Parameter: 1856 . da - the resulting distributed array object 1857 1858 Options Database Key: 1859 + -da_view - Calls DMView() at the conclusion of DMDACreate2d() 1860 . -da_grid_x <nx> - number of grid points in x direction, if M < 0 1861 . -da_grid_y <ny> - number of grid points in y direction, if N < 0 1862 . -da_processors_x <nx> - number of processors in x direction 1863 . -da_processors_y <ny> - number of processors in y direction 1864 . -da_refine_x <rx> - refinement ratio in x direction 1865 . -da_refine_y <ry> - refinement ratio in y direction 1866 - -da_refine <n> - refine the DMDA n times before creating, if M or N < 0 1867 1868 1869 Level: beginner 1870 1871 Notes: 1872 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 1873 standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 1874 the standard 9-pt stencil. 1875 1876 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 1877 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 1878 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 1879 1880 .keywords: distributed array, create, two-dimensional 1881 1882 .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(), 1883 DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(), 1884 DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges() 1885 1886 @*/ 1887 1888 PetscErrorCode DMDACreate2d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDAStencilType stencil_type, 1889 PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da) 1890 { 1891 PetscErrorCode ierr; 1892 1893 PetscFunctionBegin; 1894 ierr = DMDACreate(comm, da);CHKERRQ(ierr); 1895 ierr = DMDASetDim(*da, 2);CHKERRQ(ierr); 1896 ierr = DMDASetSizes(*da, M, N, 1);CHKERRQ(ierr); 1897 ierr = DMDASetNumProcs(*da, m, n, PETSC_DECIDE);CHKERRQ(ierr); 1898 ierr = DMDASetBoundaryType(*da, bx, by, DMDA_BOUNDARY_NONE);CHKERRQ(ierr); 1899 ierr = DMDASetDof(*da, dof);CHKERRQ(ierr); 1900 ierr = DMDASetStencilType(*da, stencil_type);CHKERRQ(ierr); 1901 ierr = DMDASetStencilWidth(*da, s);CHKERRQ(ierr); 1902 ierr = DMDASetOwnershipRanges(*da, lx, ly, PETSC_NULL);CHKERRQ(ierr); 1903 /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */ 1904 ierr = DMSetFromOptions(*da);CHKERRQ(ierr); 1905 ierr = DMSetUp(*da);CHKERRQ(ierr); 1906 ierr = DMView_DA_Private(*da);CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909