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