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