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