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! Solid Fuel Ignition (SFI) problem. This problem is modeled by 13! the partial differential equation 14! 15! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 16! 17! with boundary conditions 18! 19! u = 0 for x = 0, x = 1, y = 0, y = 1. 20! 21! A finite difference approximation with the usual 5-point stencil 22! is used to discretize the boundary value problem to obtain a nonlinear 23! system of equations. 24! 25! The uniprocessor version of this code is snes/tutorials/ex4f.F 26! 27! -------------------------------------------------------------------------- 28! The following define must be used before including any PETSc include files 29! into a module or interface. This is because they can't handle declarations 30! in them 31! 32#include <petsc/finclude/petscdmda.h> 33#include <petsc/finclude/petscsnes.h> 34#include <petsc/finclude/petscsys.h> 35#include <petsc/finclude/petscmat.h> 36module ex5f90tmodule 37 use petscdmda 38 type userctx 39 type(tDM) da 40 PetscInt xs, xe, xm, gxs, gxe, gxm 41 PetscInt ys, ye, ym, gys, gye, gym 42 PetscInt mx, my 43 PetscMPIInt rank 44 PetscReal lambda 45 end type userctx 46 47contains 48! --------------------------------------------------------------------- 49! 50! FormFunction - Evaluates nonlinear function, F(x). 51! 52! Input Parameters: 53! snes - the SNES context 54! X - input vector 55! dummy - optional user-defined context, as set by SNESSetFunction() 56! (not used here) 57! 58! Output Parameter: 59! F - function vector 60! 61! Notes: 62! This routine serves as a wrapper for the lower-level routine 63! "FormFunctionLocal", where the actual computations are 64! done using the standard Fortran style of treating the local 65! vector data as a multidimensional array over the local mesh. 66! This routine merely handles ghost point scatters and accesses 67! the local vector data via VecGetArray() and VecRestoreArray(). 68! 69 subroutine FormFunction(snesIn, X, F, user, ierr) 70 use petscsnes 71 use petscdmda 72 73! Input/output variables: 74 type(tSNES) snesIn 75 type(tVec) X, F 76 PetscErrorCode ierr 77 type(userctx) user 78 79! Declarations for use with local arrays: 80 PetscScalar, pointer :: lx_v(:), lf_v(:) 81 type(tVec) localX 82 83! Scatter ghost points to local vector, using the 2-step process 84! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 85! By placing code between these two statements, computations can 86! be done while messages are in transition. 87 PetscCall(DMGetLocalVector(user%da, localX, ierr)) 88 PetscCall(DMGlobalToLocalBegin(user%da, X, INSERT_VALUES, localX, ierr)) 89 PetscCall(DMGlobalToLocalEnd(user%da, X, INSERT_VALUES, localX, ierr)) 90 91! Get a pointer to vector data. 92! - VecGetArray90() returns a pointer to the data array. 93! - You MUST call VecRestoreArray() when you no longer need access to 94! the array. 95 96 PetscCall(VecGetArray(localX, lx_v, ierr)) 97 PetscCall(VecGetArray(F, lf_v, ierr)) 98 99! Compute function over the locally owned part of the grid 100 PetscCall(FormFunctionLocal(lx_v, lf_v, user, ierr)) 101 102! Restore vectors 103 PetscCall(VecRestoreArray(localX, lx_v, ierr)) 104 PetscCall(VecRestoreArray(F, lf_v, ierr)) 105 106! Insert values into global vector 107 108 PetscCall(DMRestoreLocalVector(user%da, localX, ierr)) 109 PetscCall(PetscLogFlops(11.0d0*user%ym*user%xm, ierr)) 110 111! PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)) 112! PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr)) 113 end subroutine formfunction 114end module ex5f90tmodule 115 116module f90moduleinterfacest 117 use ex5f90tmodule 118 119 Interface SNESSetApplicationContext 120 Subroutine SNESSetApplicationContext(snesIn, ctx, ierr) 121 use petscsnes 122 use ex5f90tmodule 123 type(tSNES) snesIn 124 type(userctx) ctx 125 PetscErrorCode ierr 126 End Subroutine 127 End Interface SNESSetApplicationContext 128 129 Interface SNESGetApplicationContext 130 Subroutine SNESGetApplicationContext(snesIn, ctx, ierr) 131 use petscsnes 132 use ex5f90tmodule 133 type(tSNES) snesIn 134 type(userctx), pointer :: ctx 135 PetscErrorCode ierr 136 End Subroutine 137 End Interface SNESGetApplicationContext 138end module f90moduleinterfacest 139 140program main 141 use petscdmda 142 use petscsnes 143 use ex5f90tmodule 144 use f90moduleinterfacest 145 implicit none 146! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147! Variable declarations 148! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 149! 150! Variables: 151! mysnes - 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 type(tSNES) mysnes 159 type(tVec) x, r 160 type(tMat) J 161 PetscErrorCode ierr 162 PetscInt its 163 PetscBool flg, matrix_free, set 164 PetscInt ione, nfour 165 PetscReal lambda_max, lambda_min 166 type(userctx) user 167 type(userctx), pointer:: puser 168 type(tPetscOptions) :: options 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 options%v = 0 182 lambda_max = 6.81 183 lambda_min = 0.0 184 user%lambda = 6.0 185 ione = 1 186 nfour = 4 187 PetscCallA(PetscOptionsGetReal(options, PETSC_NULL_CHARACTER, '-par', user%lambda, flg, ierr)) 188 PetscCheckA(user%lambda < lambda_max .and. user%lambda > lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range') 189 190! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 191! Create nonlinear solver context 192! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 193 PetscCallA(SNESCreate(PETSC_COMM_WORLD, mysnes, ierr)) 194 195! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196! Create vector data structures; set function evaluation routine 197! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 198 199! Create distributed array (DMDA) to manage parallel grid and vectors 200 201! This really needs only the star-type stencil, but we use the box 202! stencil temporarily. 203 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_ARRAY, PETSC_NULL_INTEGER_ARRAY, user%da, ierr)) 204 PetscCallA(DMSetFromOptions(user%da, ierr)) 205 PetscCallA(DMSetUp(user%da, ierr)) 206 PetscCallA(DMDAGetInfo(user%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_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr)) 207 208! 209! Visualize the distribution of the array across the processors 210! 211! PetscCallA(DMView(user%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(user%da, x, ierr)) 216 PetscCallA(VecDuplicate(x, r, ierr)) 217 218! Get local grid boundaries (for 2-dimensional DMDA) 219 PetscCallA(DMDAGetCorners(user%da, user%xs, user%ys, PETSC_NULL_INTEGER, user%xm, user%ym, PETSC_NULL_INTEGER, ierr)) 220 PetscCallA(DMDAGetGhostCorners(user%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(mysnes, user, ierr)) 235 236! Set function evaluation routine and vector 237 PetscCallA(SNESSetFunction(mysnes, 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 matrix used to construct the preconditioner 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(options, PETSC_NULL_CHARACTER, '-snes_mf', matrix_free, ierr)) 264 if (.not. matrix_free) then 265 PetscCallA(DMSetMatType(user%da, MATAIJ, ierr)) 266 PetscCallA(DMCreateMatrix(user%da, J, ierr)) 267 PetscCallA(SNESSetJacobian(mysnes, J, J, FormJacobian, user, ierr)) 268 end if 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(SNESSetFromOptions(mysnes, ierr)) 275 276! Test Fortran90 wrapper for SNESSet/Get ApplicationContext() 277 PetscCallA(PetscOptionsGetBool(options, PETSC_NULL_CHARACTER, '-test_appctx', flg, set, ierr)) 278 if (flg) then 279 PetscCallA(SNESGetApplicationContext(mysnes, puser, ierr)) 280 end if 281 282! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 283! Evaluate initial guess; then solve nonlinear system. 284! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 285! Note: The user should initialize the vector, x, with the initial guess 286! for the nonlinear solver prior to calling SNESSolve(). In particular, 287! to employ an initial guess of zero, the user should explicitly set 288! this vector to zero by calling VecSet(). 289 290 PetscCallA(FormInitialGuess(mysnes, x, ierr)) 291 PetscCallA(SNESSolve(mysnes, PETSC_NULL_VEC, x, ierr)) 292 PetscCallA(SNESGetIterationNumber(mysnes, its, ierr)) 293 if (user%rank == 0) then 294 write (6, 100) its 295 end if 296100 format('Number of SNES iterations = ', i5) 297 298! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 299! Free work space. All PETSc objects should be destroyed when they 300! are no longer needed. 301! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 302 if (.not. matrix_free) PetscCallA(MatDestroy(J, ierr)) 303 PetscCallA(VecDestroy(x, ierr)) 304 PetscCallA(VecDestroy(r, ierr)) 305 PetscCallA(SNESDestroy(mysnes, ierr)) 306 PetscCallA(DMDestroy(user%da, ierr)) 307 308 PetscCallA(PetscFinalize(ierr)) 309end 310 311! --------------------------------------------------------------------- 312! 313! FormInitialGuess - Forms initial approximation. 314! 315! Input Parameters: 316! X - vector 317! 318! Output Parameter: 319! X - vector 320! 321! Notes: 322! This routine serves as a wrapper for the lower-level routine 323! "InitialGuessLocal", where the actual computations are 324! done using the standard Fortran style of treating the local 325! vector data as a multidimensional array over the local mesh. 326! This routine merely handles ghost point scatters and accesses 327! the local vector data via VecGetArray() and VecRestoreArray(). 328! 329subroutine FormInitialGuess(mysnes, X, ierr) 330 use petscsnes 331 use ex5f90tmodule 332 use f90moduleinterfacest 333! Input/output variables: 334 type(tSNES) mysnes 335 type(userctx), pointer:: puser 336 type(tVec) X 337 PetscErrorCode ierr 338 339! Declarations for use with local arrays: 340 PetscScalar, pointer :: lx_v(:) 341 342 ierr = 0 343 PetscCallA(SNESGetApplicationContext(mysnes, puser, ierr)) 344! Get a pointer to vector data. 345! - VecGetArray90() returns a pointer to the data array. 346! - You MUST call VecRestoreArray() when you no longer need access to 347! the array. 348 349 PetscCallA(VecGetArray(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(VecRestoreArray(X, lx_v, ierr)) 356 357! Insert values into global vector 358 359end 360 361! --------------------------------------------------------------------- 362! 363! InitialGuessLocal - Computes initial approximation, called by 364! the higher level routine FormInitialGuess(). 365! 366! Input Parameter: 367! x - local vector data 368! 369! Output Parameters: 370! x - local vector data 371! ierr - error code 372! 373! Notes: 374! This routine uses standard Fortran-style computations over a 2-dim array. 375! 376subroutine InitialGuessLocal(user, x, ierr) 377 use petscsys 378 use ex5f90tmodule 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 PetscScalar temp1, temp, hx, hy 387 PetscScalar one 388 389! Set parameters 390 391 ierr = 0 392 one = 1.0 393 hx = one/(PetscIntToReal(user%mx - 1)) 394 hy = one/(PetscIntToReal(user%my - 1)) 395 temp1 = user%lambda/(user%lambda + one) 396 397 do 20 j = user%ys, user%ye 398 temp = PetscIntToReal(min(j - 1, user%my - j))*hy 399 do 10 i = user%xs, user%xe 400 if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 401 x(i, j) = 0.0 402 else 403 x(i, j) = temp1*sqrt(min(PetscIntToReal(min(i - 1, user%mx - i)*hx), PetscIntToReal(temp))) 404 end if 40510 continue 40620 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 petscsys 427 use ex5f90tmodule 428! Input/output variables: 429 type(userctx) user 430 PetscScalar x(user%gxs:user%gxe, user%gys:user%gye) 431 PetscScalar f(user%xs:user%xe, user%ys:user%ye) 432 PetscErrorCode ierr 433 434! Local variables: 435 PetscScalar two, one, hx, hy, hxdhy, hydhx, sc 436 PetscScalar u, uxx, uyy 437 PetscInt i, j 438 439 one = 1.0 440 two = 2.0 441 hx = one/PetscIntToReal(user%mx - 1) 442 hy = one/PetscIntToReal(user%my - 1) 443 sc = hx*hy*user%lambda 444 hxdhy = hx/hy 445 hydhx = hy/hx 446 447! Compute function over the locally owned part of the grid 448 449 do 20 j = user%ys, user%ye 450 do 10 i = user%xs, user%xe 451 if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 452 f(i, j) = x(i, j) 453 else 454 u = x(i, j) 455 uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j)) 456 uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1)) 457 f(i, j) = uxx + uyy - sc*exp(u) 458 end if 45910 continue 46020 continue 461 ierr = 0 462 end 463 464! --------------------------------------------------------------------- 465! 466! FormJacobian - Evaluates Jacobian matrix. 467! 468! Input Parameters: 469! snes - the SNES context 470! x - input vector 471! dummy - optional user-defined context, as set by SNESSetJacobian() 472! (not used here) 473! 474! Output Parameters: 475! jac - Jacobian matrix 476! jac_prec - optionally different matrix used to construct the preconditioner (not used here) 477! 478! Notes: 479! This routine serves as a wrapper for the lower-level routine 480! "FormJacobianLocal", where the actual computations are 481! done using the standard Fortran style of treating the local 482! vector data as a multidimensional array over the local mesh. 483! This routine merely accesses the local vector data via 484! VecGetArray() and VecRestoreArray(). 485! 486! Notes: 487! Due to grid point reordering with DMDAs, we must always work 488! with the local grid points, and then transform them to the new 489! global numbering with the "ltog" mapping 490! We cannot work directly with the global numbers for the original 491! uniprocessor grid! 492! 493! Two methods are available for imposing this transformation 494! when setting matrix entries: 495! (A) MatSetValuesLocal(), using the local ordering (including 496! ghost points!) 497! - Set matrix entries using the local ordering 498! by calling MatSetValuesLocal() 499! (B) MatSetValues(), using the global ordering 500! - Use DMGetLocalToGlobalMapping() then 501! ISLocalToGlobalMappingGetIndices() to extract the local-to-global map 502! - Then apply this map explicitly yourself 503! - Set matrix entries using the global ordering by calling 504! MatSetValues() 505! Option (A) seems cleaner/easier in many cases, and is the procedure 506! used in this example. 507! 508 subroutine FormJacobian(mysnes, X, jac, jac_prec, user, ierr) 509 use petscsnes 510 use ex5f90tmodule 511! Input/output variables: 512 type(tSNES) mysnes 513 type(tVec) X 514 type(tMat) jac, jac_prec 515 type(userctx) user 516 PetscErrorCode ierr 517 518! Declarations for use with local arrays: 519 PetscScalar, pointer :: lx_v(:) 520 type(tVec) localX 521 522! Scatter ghost points to local vector, using the 2-step process 523! DMGlobalToLocalBegin(), DMGlobalToLocalEnd() 524! Computations can be done while messages are in transition, 525! by placing code between these two statements. 526 527 PetscCallA(DMGetLocalVector(user%da, localX, ierr)) 528 PetscCallA(DMGlobalToLocalBegin(user%da, X, INSERT_VALUES, localX, ierr)) 529 PetscCallA(DMGlobalToLocalEnd(user%da, X, INSERT_VALUES, localX, ierr)) 530 531! Get a pointer to vector data 532 PetscCallA(VecGetArray(localX, lx_v, ierr)) 533 534! Compute entries for the locally owned part of the Jacobian preconditioner. 535 PetscCallA(FormJacobianLocal(lx_v, jac_prec, user, ierr)) 536 537! Assemble matrix, using the 2-step process: 538! MatAssemblyBegin(), MatAssemblyEnd() 539! Computations can be done while messages are in transition, 540! by placing code between these two statements. 541 542 PetscCallA(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)) 543! if (jac .ne. jac_prec) then 544 PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 545! endif 546 PetscCallA(VecRestoreArray(localX, lx_v, ierr)) 547 PetscCallA(DMRestoreLocalVector(user%da, localX, ierr)) 548 PetscCallA(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)) 549! if (jac .ne. jac_prec) then 550 PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 551! endif 552 553! Tell the matrix we will never add a new nonzero location to the 554! matrix. If we do it will generate an error. 555 556 PetscCallA(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr)) 557 558 end 559 560! --------------------------------------------------------------------- 561! 562! FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner, 563! called by the higher level routine FormJacobian(). 564! 565! Input Parameters: 566! x - local vector data 567! 568! Output Parameters: 569! jac_prec - Jacobian matrix used to compute the preconditioner 570! ierr - error code 571! 572! Notes: 573! This routine uses standard Fortran-style computations over a 2-dim array. 574! 575! Notes: 576! Due to grid point reordering with DMDAs, we must always work 577! with the local grid points, and then transform them to the new 578! global numbering with the "ltog" mapping 579! We cannot work directly with the global numbers for the original 580! uniprocessor grid! 581! 582! Two methods are available for imposing this transformation 583! when setting matrix entries: 584! (A) MatSetValuesLocal(), using the local ordering (including 585! ghost points!) 586! - Set matrix entries using the local ordering 587! by calling MatSetValuesLocal() 588! (B) MatSetValues(), using the global ordering 589! - Set matrix entries using the global ordering by calling 590! MatSetValues() 591! Option (A) seems cleaner/easier in many cases, and is the procedure 592! used in this example. 593! 594 subroutine FormJacobianLocal(x, jac_prec, user, ierr) 595 use petscmat 596 use ex5f90tmodule 597! Input/output variables: 598 type(userctx) user 599 PetscScalar x(user%gxs:user%gxe, user%gys:user%gye) 600 type(tMat) jac_prec 601 PetscErrorCode ierr 602 603! Local variables: 604 PetscInt row, col(5), i, j 605 PetscInt ione, ifive 606 PetscScalar two, one, hx, hy, hxdhy 607 PetscScalar hydhx, sc, v(5) 608 609! Set parameters 610 ione = 1 611 ifive = 5 612 one = 1.0 613 two = 2.0 614 hx = one/PetscIntToReal(user%mx - 1) 615 hy = one/PetscIntToReal(user%my - 1) 616 sc = hx*hy 617 hxdhy = hx/hy 618 hydhx = hy/hx 619 620! Compute entries for the locally owned part of the Jacobian. 621! - Currently, all PETSc parallel matrix formats are partitioned by 622! contiguous chunks of rows across the processors. 623! - Each processor needs to insert only elements that it owns 624! locally (but any non-local elements will be sent to the 625! appropriate processor during matrix assembly). 626! - Here, we set all entries for a particular row at once. 627! - We can set matrix entries either using either 628! MatSetValuesLocal() or MatSetValues(), as discussed above. 629! - Note that MatSetValues() uses 0-based row and column numbers 630! in Fortran as well as in C. 631 632 do 20 j = user%ys, user%ye 633 row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1 634 do 10 i = user%xs, user%xe 635 row = row + 1 636! boundary points 637 if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 638 col(1) = row 639 v(1) = one 640 PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ione, col, v, INSERT_VALUES, ierr)) 641! interior grid points 642 else 643 v(1) = -hxdhy 644 v(2) = -hydhx 645 v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i, j)) 646 v(4) = -hydhx 647 v(5) = -hxdhy 648 col(1) = row - user%gxm 649 col(2) = row - 1 650 col(3) = row 651 col(4) = row + 1 652 col(5) = row + user%gxm 653 PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ifive, col, v, INSERT_VALUES, ierr)) 654 end if 65510 continue 65620 continue 657 end 658 659!/*TEST 660! 661! test: 662! nsize: 4 663! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 664! 665!TEST*/ 666