1! 2! Description: Solves a nonlinear system in parallel with SNES. 3! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular 4! domain, using distributed arrays (DMDAs) to partition the parallel grid. 5! The command line options include: 6! -par <parameter>, where <parameter> indicates the nonlinearity of the problem 7! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 8! 9 10! 11! -------------------------------------------------------------------------- 12! 13! Solid Fuel Ignition (SFI) problem. This problem is modeled by 14! the partial differential equation 15! 16! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 17! 18! with boundary conditions 19! 20! u = 0 for x = 0, x = 1, y = 0, y = 1. 21! 22! A finite difference approximation with the usual 5-point stencil 23! is used to discretize the boundary value problem to obtain a nonlinear 24! system of equations. 25! 26! The uniprocessor version of this code is snes/tutorials/ex4f.F 27! 28! -------------------------------------------------------------------------- 29! The following define must be used before including any PETSc include files 30! into a module or interface. This is because they can't handle declarations 31! in them 32! 33 34 module ex5f90module 35 use petscsnes 36 use petscdmda 37#include <petsc/finclude/petscsnes.h> 38 type userctx 39 PetscInt xs,xe,xm,gxs,gxe,gxm 40 PetscInt ys,ye,ym,gys,gye,gym 41 PetscInt mx,my 42 PetscMPIInt rank 43 PetscReal lambda 44 end type userctx 45 46 contains 47! --------------------------------------------------------------------- 48! 49! FormFunction - Evaluates nonlinear function, F(x). 50! 51! Input Parameters: 52! snes - the SNES context 53! X - input vector 54! dummy - optional user-defined context, as set by SNESSetFunction() 55! (not used here) 56! 57! Output Parameter: 58! F - function vector 59! 60! Notes: 61! This routine serves as a wrapper for the lower-level routine 62! "FormFunctionLocal", where the actual computations are 63! done using the standard Fortran style of treating the local 64! vector data as a multidimensional array over the local mesh. 65! This routine merely handles ghost point scatters and accesses 66! the local vector data via VecGetArrayF90() and VecRestoreArrayF90(). 67! 68 subroutine FormFunction(snes,X,F,user,ierr) 69 implicit none 70 71! Input/output variables: 72 SNES snes 73 Vec X,F 74 PetscErrorCode ierr 75 type (userctx) user 76 DM da 77 78! Declarations for use with local arrays: 79 PetscScalar,pointer :: lx_v(:),lf_v(:) 80 Vec localX 81 82! Scatter ghost points to local vector, using the 2-step process 83! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 84! By placing code between these two statements, computations can 85! be done while messages are in transition. 86 PetscCall(SNESGetDM(snes,da,ierr)) 87 PetscCall(DMGetLocalVector(da,localX,ierr)) 88 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)) 89 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)) 90 91! Get a pointer to vector data. 92! - For default PETSc vectors, VecGetArrayF90() returns a pointer to 93! the data array. Otherwise, the routine is implementation dependent. 94! - You MUST call VecRestoreArrayF90() when you no longer need access to 95! the array. 96! - Note that the interface to VecGetArrayF90() differs from VecGetArray(). 97 98 PetscCall(VecGetArrayF90(localX,lx_v,ierr)) 99 PetscCall(VecGetArrayF90(F,lf_v,ierr)) 100 101! Compute function over the locally owned part of the grid 102 PetscCall(FormFunctionLocal(lx_v,lf_v,user,ierr)) 103 104! Restore vectors 105 PetscCall(VecRestoreArrayF90(localX,lx_v,ierr)) 106 PetscCall(VecRestoreArrayF90(F,lf_v,ierr)) 107 108! Insert values into global vector 109 110 PetscCall(DMRestoreLocalVector(da,localX,ierr)) 111 PetscCall(PetscLogFlops(11.0d0*user%ym*user%xm,ierr)) 112 113! PetscCallA(VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)) 114! PetscCallA(VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr)) 115 end subroutine formfunction 116 end module ex5f90module 117 118 module ex5f90moduleinterfaces 119 use ex5f90module 120 121 Interface SNESSetApplicationContext 122 Subroutine SNESSetApplicationContext(snes,ctx,ierr) 123 use ex5f90module 124 SNES snes 125 type(userctx) ctx 126 PetscErrorCode ierr 127 End Subroutine 128 End Interface SNESSetApplicationContext 129 130 Interface SNESGetApplicationContext 131 Subroutine SNESGetApplicationContext(snes,ctx,ierr) 132 use ex5f90module 133 SNES snes 134 type(userctx), pointer :: ctx 135 PetscErrorCode ierr 136 End Subroutine 137 End Interface SNESGetApplicationContext 138 end module ex5f90moduleinterfaces 139 140 program main 141 use ex5f90module 142 use ex5f90moduleinterfaces 143 implicit none 144! 145 146! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147! Variable declarations 148! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 149! 150! Variables: 151! snes - nonlinear solver 152! x, r - solution, residual vectors 153! J - Jacobian matrix 154! its - iterations for convergence 155! Nx, Ny - number of preocessors in x- and y- directions 156! matrix_free - flag - 1 indicates matrix-free version 157! 158 SNES snes 159 Vec x,r 160 Mat J 161 PetscErrorCode ierr 162 PetscInt its 163 PetscBool flg,matrix_free 164 PetscInt ione,nfour 165 PetscReal lambda_max,lambda_min 166 type (userctx) user 167 DM da 168 169! Note: Any user-defined Fortran routines (such as FormJacobian) 170! MUST be declared as external. 171 external FormInitialGuess,FormJacobian 172 173! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 174! Initialize program 175! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 176 PetscCallA(PetscInitialize(ierr)) 177 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,user%rank,ierr)) 178 179! Initialize problem parameters 180 lambda_max = 6.81 181 lambda_min = 0.0 182 user%lambda = 6.0 183 ione = 1 184 nfour = 4 185 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',user%lambda,flg,ierr)) 186 PetscCheckA(user%lambda .lt. lambda_max .and. user%lambda .gt. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range') 187 188! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 189! Create nonlinear solver context 190! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 191 PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr)) 192 193! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194! Create vector data structures; set function evaluation routine 195! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196 197! Create distributed array (DMDA) to manage parallel grid and vectors 198 199! This really needs only the star-type stencil, but we use the box 200! stencil temporarily. 201 PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)) 202 PetscCallA(DMSetFromOptions(da,ierr)) 203 PetscCallA(DMSetUp(da,ierr)) 204 205 PetscCallA(DMDAGetInfo(da,PETSC_NULL_INTEGER,user%mx,user%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 206 207! 208! Visualize the distribution of the array across the processors 209! 210! PetscCallA(DMView(da,PETSC_VIEWER_DRAW_WORLD,ierr)) 211 212! Extract global and local vectors from DMDA; then duplicate for remaining 213! vectors that are the same types 214 PetscCallA(DMCreateGlobalVector(da,x,ierr)) 215 PetscCallA(VecDuplicate(x,r,ierr)) 216 217! Get local grid boundaries (for 2-dimensional DMDA) 218 PetscCallA(DMDAGetCorners(da,user%xs,user%ys,PETSC_NULL_INTEGER,user%xm,user%ym,PETSC_NULL_INTEGER,ierr)) 219 PetscCallA(DMDAGetGhostCorners(da,user%gxs,user%gys,PETSC_NULL_INTEGER,user%gxm,user%gym,PETSC_NULL_INTEGER,ierr)) 220 221! Here we shift the starting indices up by one so that we can easily 222! use the Fortran convention of 1-based indices (rather 0-based indices). 223 user%xs = user%xs+1 224 user%ys = user%ys+1 225 user%gxs = user%gxs+1 226 user%gys = user%gys+1 227 228 user%ye = user%ys+user%ym-1 229 user%xe = user%xs+user%xm-1 230 user%gye = user%gys+user%gym-1 231 user%gxe = user%gxs+user%gxm-1 232 233 PetscCallA(SNESSetApplicationContext(snes,user,ierr)) 234 235! Set function evaluation routine and vector 236 PetscCallA(SNESSetFunction(snes,r,FormFunction,user,ierr)) 237 238! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 239! Create matrix data structure; set Jacobian evaluation routine 240! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 241 242! Set Jacobian matrix data structure and default Jacobian evaluation 243! routine. User can override with: 244! -snes_fd : default finite differencing approximation of Jacobian 245! -snes_mf : matrix-free Newton-Krylov method with no preconditioning 246! (unless user explicitly sets preconditioner) 247! -snes_mf_operator : form preconditioning matrix as set by the user, 248! but use matrix-free approx for Jacobian-vector 249! products within Newton-Krylov method 250! 251! Note: For the parallel case, vectors and matrices MUST be partitioned 252! accordingly. When using distributed arrays (DMDAs) to create vectors, 253! the DMDAs determine the problem partitioning. We must explicitly 254! specify the local matrix dimensions upon its creation for compatibility 255! with the vector distribution. Thus, the generic MatCreate() routine 256! is NOT sufficient when working with distributed arrays. 257! 258! Note: Here we only approximately preallocate storage space for the 259! Jacobian. See the users manual for a discussion of better techniques 260! for preallocating matrix memory. 261 262 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)) 263 if (.not. matrix_free) then 264 PetscCallA(DMSetMatType(da,MATAIJ,ierr)) 265 PetscCallA(DMCreateMatrix(da,J,ierr)) 266 PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,user,ierr)) 267 endif 268 269! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 270! Customize nonlinear solver; set runtime options 271! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 272! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 273 PetscCallA(SNESSetDM(snes,da,ierr)) 274 PetscCallA(SNESSetFromOptions(snes,ierr)) 275 276! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 277! Evaluate initial guess; then solve nonlinear system. 278! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 279! Note: The user should initialize the vector, x, with the initial guess 280! for the nonlinear solver prior to calling SNESSolve(). In particular, 281! to employ an initial guess of zero, the user should explicitly set 282! this vector to zero by calling VecSet(). 283 284 PetscCallA(FormInitialGuess(snes,x,ierr)) 285 PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr)) 286 PetscCallA(SNESGetIterationNumber(snes,its,ierr)) 287 if (user%rank .eq. 0) then 288 write(6,100) its 289 endif 290 100 format('Number of SNES iterations = ',i5) 291 292! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 293! Free work space. All PETSc objects should be destroyed when they 294! are no longer needed. 295! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 296 if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr)) 297 PetscCallA(VecDestroy(x,ierr)) 298 PetscCallA(VecDestroy(r,ierr)) 299 PetscCallA(SNESDestroy(snes,ierr)) 300 PetscCallA(DMDestroy(da,ierr)) 301 302 PetscCallA(PetscFinalize(ierr)) 303 end 304 305! --------------------------------------------------------------------- 306! 307! FormInitialGuess - Forms initial approximation. 308! 309! Input Parameters: 310! X - vector 311! 312! Output Parameter: 313! X - vector 314! 315! Notes: 316! This routine serves as a wrapper for the lower-level routine 317! "InitialGuessLocal", where the actual computations are 318! done using the standard Fortran style of treating the local 319! vector data as a multidimensional array over the local mesh. 320! This routine merely handles ghost point scatters and accesses 321! the local vector data via VecGetArrayF90() and VecRestoreArrayF90(). 322! 323 subroutine FormInitialGuess(snes,X,ierr) 324 use ex5f90module 325 use ex5f90moduleinterfaces 326 implicit none 327 328! Input/output variables: 329 SNES snes 330 type(userctx), pointer:: puser 331 Vec X 332 PetscErrorCode ierr 333 DM da 334 335! Declarations for use with local arrays: 336 PetscScalar,pointer :: lx_v(:) 337 338 ierr = 0 339 PetscCallA(SNESGetDM(snes,da,ierr)) 340 PetscCallA(SNESGetApplicationContext(snes,puser,ierr)) 341! Get a pointer to vector data. 342! - For default PETSc vectors, VecGetArrayF90() returns a pointer to 343! the data array. Otherwise, the routine is implementation dependent. 344! - You MUST call VecRestoreArrayF90() when you no longer need access to 345! the array. 346! - Note that the interface to VecGetArrayF90() differs from VecGetArray(). 347 348 PetscCallA(VecGetArrayF90(X,lx_v,ierr)) 349 350! Compute initial guess over the locally owned part of the grid 351 PetscCallA(InitialGuessLocal(puser,lx_v,ierr)) 352 353! Restore vector 354 PetscCallA(VecRestoreArrayF90(X,lx_v,ierr)) 355 356! Insert values into global vector 357 358 end 359 360! --------------------------------------------------------------------- 361! 362! InitialGuessLocal - Computes initial approximation, called by 363! the higher level routine FormInitialGuess(). 364! 365! Input Parameter: 366! x - local vector data 367! 368! Output Parameters: 369! x - local vector data 370! ierr - error code 371! 372! Notes: 373! This routine uses standard Fortran-style computations over a 2-dim array. 374! 375 subroutine InitialGuessLocal(user,x,ierr) 376 use ex5f90module 377 implicit none 378 379! Input/output variables: 380 type (userctx) user 381 PetscScalar x(user%xs:user%xe,user%ys:user%ye) 382 PetscErrorCode ierr 383 384! Local variables: 385 PetscInt i,j 386 PetscReal temp1,temp,hx,hy 387 PetscReal one 388 389! Set parameters 390 391 ierr = 0 392 one = 1.0 393 hx = one/(user%mx-1) 394 hy = one/(user%my-1) 395 temp1 = user%lambda/(user%lambda + one) 396 397 do 20 j=user%ys,user%ye 398 temp = min(j-1,user%my-j)*hy 399 do 10 i=user%xs,user%xe 400 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then 401 x(i,j) = 0.0 402 else 403 x(i,j) = temp1 * sqrt(min(hx*min(i-1,user%mx-i),temp)) 404 endif 405 10 continue 406 20 continue 407 408 end 409 410! --------------------------------------------------------------------- 411! 412! FormFunctionLocal - Computes nonlinear function, called by 413! the higher level routine FormFunction(). 414! 415! Input Parameter: 416! x - local vector data 417! 418! Output Parameters: 419! f - local vector data, f(x) 420! ierr - error code 421! 422! Notes: 423! This routine uses standard Fortran-style computations over a 2-dim array. 424! 425 subroutine FormFunctionLocal(x,f,user,ierr) 426 use ex5f90module 427 428 implicit none 429 430! Input/output variables: 431 type (userctx) user 432 PetscScalar x(user%gxs:user%gxe,user%gys:user%gye) 433 PetscScalar f(user%xs:user%xe,user%ys:user%ye) 434 PetscErrorCode ierr 435 436! Local variables: 437 PetscScalar two,one,hx,hy,hxdhy,hydhx,sc 438 PetscScalar u,uxx,uyy 439 PetscInt i,j 440 441 one = 1.0 442 two = 2.0 443 hx = one/(user%mx-1) 444 hy = one/(user%my-1) 445 sc = hx*hy*user%lambda 446 hxdhy = hx/hy 447 hydhx = hy/hx 448 449! Compute function over the locally owned part of the grid 450 451 do 20 j=user%ys,user%ye 452 do 10 i=user%xs,user%xe 453 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then 454 f(i,j) = x(i,j) 455 else 456 u = x(i,j) 457 uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j)) 458 uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1)) 459 f(i,j) = uxx + uyy - sc*exp(u) 460 endif 461 10 continue 462 20 continue 463 464 end 465 466! --------------------------------------------------------------------- 467! 468! FormJacobian - Evaluates Jacobian matrix. 469! 470! Input Parameters: 471! snes - the SNES context 472! x - input vector 473! dummy - optional user-defined context, as set by SNESSetJacobian() 474! (not used here) 475! 476! Output Parameters: 477! jac - Jacobian matrix 478! jac_prec - optionally different preconditioning matrix (not used here) 479! flag - flag indicating matrix structure 480! 481! Notes: 482! This routine serves as a wrapper for the lower-level routine 483! "FormJacobianLocal", where the actual computations are 484! done using the standard Fortran style of treating the local 485! vector data as a multidimensional array over the local mesh. 486! This routine merely accesses the local vector data via 487! VecGetArrayF90() and VecRestoreArrayF90(). 488! 489! Notes: 490! Due to grid point reordering with DMDAs, we must always work 491! with the local grid points, and then transform them to the new 492! global numbering with the "ltog" mapping 493! We cannot work directly with the global numbers for the original 494! uniprocessor grid! 495! 496! Two methods are available for imposing this transformation 497! when setting matrix entries: 498! (A) MatSetValuesLocal(), using the local ordering (including 499! ghost points!) 500! - Set matrix entries using the local ordering 501! by calling MatSetValuesLocal() 502! (B) MatSetValues(), using the global ordering 503 504! - Set matrix entries using the global ordering by calling 505! MatSetValues() 506! Option (A) seems cleaner/easier in many cases, and is the procedure 507! used in this example. 508! 509 subroutine FormJacobian(snes,X,jac,jac_prec,user,ierr) 510 use ex5f90module 511 implicit none 512 513! Input/output variables: 514 SNES snes 515 Vec X 516 Mat jac,jac_prec 517 type(userctx) user 518 PetscErrorCode ierr 519 DM da 520 521! Declarations for use with local arrays: 522 PetscScalar,pointer :: lx_v(:) 523 Vec localX 524 525! Scatter ghost points to local vector, using the 2-step process 526! DMGlobalToLocalBegin(), DMGlobalToLocalEnd() 527! Computations can be done while messages are in transition, 528! by placing code between these two statements. 529 530 PetscCallA(SNESGetDM(snes,da,ierr)) 531 PetscCallA(DMGetLocalVector(da,localX,ierr)) 532 PetscCallA(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)) 533 PetscCallA(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)) 534 535! Get a pointer to vector data 536 PetscCallA(VecGetArrayF90(localX,lx_v,ierr)) 537 538! Compute entries for the locally owned part of the Jacobian preconditioner. 539 PetscCallA(FormJacobianLocal(lx_v,jac_prec,user,ierr)) 540 541! Assemble matrix, using the 2-step process: 542! MatAssemblyBegin(), MatAssemblyEnd() 543! Computations can be done while messages are in transition, 544! by placing code between these two statements. 545 546 PetscCallA(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 547 if (jac .ne. jac_prec) then 548 PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 549 endif 550 PetscCallA(VecRestoreArrayF90(localX,lx_v,ierr)) 551 PetscCallA(DMRestoreLocalVector(da,localX,ierr)) 552 PetscCallA(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 553 if (jac .ne. jac_prec) then 554 PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 555 endif 556 557! Tell the matrix we will never add a new nonzero location to the 558! matrix. If we do it will generate an error. 559 560 PetscCallA(MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 561 562 end 563 564! --------------------------------------------------------------------- 565! 566! FormJacobianLocal - Computes Jacobian preconditioner matrix, 567! called by the higher level routine FormJacobian(). 568! 569! Input Parameters: 570! x - local vector data 571! 572! Output Parameters: 573! jac_prec - Jacobian preconditioner matrix 574! ierr - error code 575! 576! Notes: 577! This routine uses standard Fortran-style computations over a 2-dim array. 578! 579! Notes: 580! Due to grid point reordering with DMDAs, we must always work 581! with the local grid points, and then transform them to the new 582! global numbering with the "ltog" mapping 583! We cannot work directly with the global numbers for the original 584! uniprocessor grid! 585! 586! Two methods are available for imposing this transformation 587! when setting matrix entries: 588! (A) MatSetValuesLocal(), using the local ordering (including 589! ghost points!) 590! - Set matrix entries using the local ordering 591! by calling MatSetValuesLocal() 592! (B) MatSetValues(), using the global ordering 593! - Then apply this map explicitly yourself 594! - Set matrix entries using the global ordering by calling 595! MatSetValues() 596! Option (A) seems cleaner/easier in many cases, and is the procedure 597! used in this example. 598! 599 subroutine FormJacobianLocal(x,jac_prec,user,ierr) 600 use ex5f90module 601 implicit none 602 603! Input/output variables: 604 type (userctx) user 605 PetscScalar x(user%gxs:user%gxe,user%gys:user%gye) 606 Mat jac_prec 607 PetscErrorCode ierr 608 609! Local variables: 610 PetscInt row,col(5),i,j 611 PetscInt ione,ifive 612 PetscScalar two,one,hx,hy,hxdhy 613 PetscScalar hydhx,sc,v(5) 614 615! Set parameters 616 ione = 1 617 ifive = 5 618 one = 1.0 619 two = 2.0 620 hx = one/(user%mx-1) 621 hy = one/(user%my-1) 622 sc = hx*hy 623 hxdhy = hx/hy 624 hydhx = hy/hx 625 626! Compute entries for the locally owned part of the Jacobian. 627! - Currently, all PETSc parallel matrix formats are partitioned by 628! contiguous chunks of rows across the processors. 629! - Each processor needs to insert only elements that it owns 630! locally (but any non-local elements will be sent to the 631! appropriate processor during matrix assembly). 632! - Here, we set all entries for a particular row at once. 633! - We can set matrix entries either using either 634! MatSetValuesLocal() or MatSetValues(), as discussed above. 635! - Note that MatSetValues() uses 0-based row and column numbers 636! in Fortran as well as in C. 637 638 do 20 j=user%ys,user%ye 639 row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1 640 do 10 i=user%xs,user%xe 641 row = row + 1 642! boundary points 643 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then 644 col(1) = row 645 v(1) = one 646 PetscCallA(MatSetValuesLocal(jac_prec,ione,row,ione,col,v,INSERT_VALUES,ierr)) 647! interior grid points 648 else 649 v(1) = -hxdhy 650 v(2) = -hydhx 651 v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i,j)) 652 v(4) = -hydhx 653 v(5) = -hxdhy 654 col(1) = row - user%gxm 655 col(2) = row - 1 656 col(3) = row 657 col(4) = row + 1 658 col(5) = row + user%gxm 659 PetscCallA(MatSetValuesLocal(jac_prec,ione,row,ifive,col,v,INSERT_VALUES,ierr)) 660 endif 661 10 continue 662 20 continue 663 664 end 665 666! 667!/*TEST 668! 669! test: 670! nsize: 4 671! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 672! requires: !single 673! 674! test: 675! suffix: 2 676! nsize: 4 677! args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 678! requires: !single 679! 680! test: 681! suffix: 3 682! nsize: 3 683! args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 684! requires: !single 685! 686! test: 687! suffix: 4 688! nsize: 3 689! args: -snes_mf_operator -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 690! requires: !single 691! 692! test: 693! suffix: 5 694! requires: !single 695! 696!TEST*/ 697