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