1! 2! This example shows how to avoid Fortran line lengths larger than 132 characters. 3! It avoids used of certain macros such as PetscCallA() and PetscCheckA() that 4! generate very long lines 5! 6! We recommend starting from src/snes/tutorials/ex5f90.F90 instead of this example 7! because that does not have the restricted formatting that makes this version 8! more difficult to read 9! 10! Description: This example solves a nonlinear system in parallel with SNES. 11! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular 12! domain, using distributed arrays (DMDAs) to partition the parallel grid. 13! The command line options include: 14! -par <param>, where <param> indicates the nonlinearity of the problem 15! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 16! 17! -------------------------------------------------------------------------- 18! 19! Solid Fuel Ignition (SFI) problem. This problem is modeled by 20! the partial differential equation 21! 22! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 23! 24! with boundary conditions 25! 26! u = 0 for x = 0, x = 1, y = 0, y = 1. 27! 28! A finite difference approximation with the usual 5-point stencil 29! is used to discretize the boundary value problem to obtain a nonlinear 30! system of equations. 31! 32! -------------------------------------------------------------------------- 33#include <petsc/finclude/petscsnes.h> 34#include <petsc/finclude/petscdmda.h> 35module ex5fmodule 36 use petscsnes 37 use petscdmda 38 implicit none 39 PetscInt xs, xe, xm, gxs, gxe, gxm 40 PetscInt ys, ye, ym, gys, gye, gym 41 PetscInt mx, my 42 PetscMPIInt rank, size 43 PetscReal lambda 44contains 45! --------------------------------------------------------------------- 46! 47! FormInitialGuess - Forms initial approximation. 48! 49! Input Parameters: 50! X - vector 51! 52! Output Parameter: 53! X - vector 54! 55! Notes: 56! This routine serves as a wrapper for the lower-level routine 57! "ApplicationInitialGuess", where the actual computations are 58! done using the standard Fortran style of treating the local 59! vector data as a multidimensional array over the local mesh. 60! This routine merely handles ghost point scatters and accesses 61! the local vector data via VecGetArray() and VecRestoreArray(). 62! 63 subroutine FormInitialGuess(X, ierr) 64 65! Input/output variables: 66 Vec X 67 PetscErrorCode ierr 68! Declarations for use with local arrays: 69 PetscScalar, pointer :: lx_v(:) 70 71 ierr = 0 72 73! Get a pointer to vector data. 74! - For default PETSc vectors, VecGetArray() returns a pointer to 75! the data array. Otherwise, the routine is implementation dependent. 76! - You MUST call VecRestoreArray() when you no longer need access to 77! the array. 78! - Note that the Fortran interface to VecGetArray() differs from the 79! C version. See the users manual for details. 80 81 call VecGetArray(X, lx_v, ierr) 82 CHKERRQ(ierr) 83 84! Compute initial guess over the locally owned part of the grid 85 86 call InitialGuessLocal(lx_v, ierr) 87 CHKERRQ(ierr) 88 89! Restore vector 90 91 call VecRestoreArray(X, lx_v, ierr) 92 CHKERRQ(ierr) 93 94 end 95 96! --------------------------------------------------------------------- 97! 98! InitialGuessLocal - Computes initial approximation, called by 99! the higher level routine FormInitialGuess(). 100! 101! Input Parameter: 102! x - local vector data 103! 104! Output Parameters: 105! x - local vector data 106! ierr - error code 107! 108! Notes: 109! This routine uses standard Fortran-style computations over a 2-dim array. 110! 111 subroutine InitialGuessLocal(x, ierr) 112 113! Input/output variables: 114 PetscScalar x(xs:xe, ys:ye) 115 PetscErrorCode ierr 116 117! Local variables: 118 PetscInt i, j 119 PetscReal temp1, temp, one, hx, hy 120 121! Set parameters 122 123 ierr = 0 124 one = 1.0 125 hx = one/((real(mx) - 1)) 126 hy = one/((real(my) - 1)) 127 temp1 = lambda/(lambda + one) 128 129 do j = ys, ye 130 temp = (real(min(j - 1, my - j)))*hy 131 do i = xs, xe 132 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 133 x(i, j) = 0.0 134 else 135 x(i, j) = temp1*sqrt(min(real(min(i - 1, mx - i))*hx, (temp))) 136 end if 137 end do 138 end do 139 140 end 141 142! --------------------------------------------------------------------- 143! 144! FormFunctionLocal - Computes nonlinear function, called by 145! the higher level routine FormFunction(). 146! 147! Input Parameter: 148! x - local vector data 149! 150! Output Parameters: 151! f - local vector data, f(x) 152! ierr - error code 153! 154! Notes: 155! This routine uses standard Fortran-style computations over a 2-dim array. 156! 157! 158 subroutine FormFunctionLocal(info, x, f, da, ierr) 159 160 DM da 161 162! Input/output variables: 163 DMDALocalInfo info 164 PetscScalar x(gxs:gxe, gys:gye) 165 PetscScalar f(xs:xe, ys:ye) 166 PetscErrorCode ierr 167 168! Local variables: 169 PetscScalar two, one, hx, hy 170 PetscScalar hxdhy, hydhx, sc 171 PetscScalar u, uxx, uyy 172 PetscInt i, j 173 174 xs = info%XS + 1 175 xe = xs + info%XM - 1 176 ys = info%YS + 1 177 ye = ys + info%YM - 1 178 mx = info%MX 179 my = info%MY 180 181 one = 1.0 182 two = 2.0 183 hx = one/(real(mx) - 1) 184 hy = one/(real(my) - 1) 185 sc = hx*hy*lambda 186 hxdhy = hx/hy 187 hydhx = hy/hx 188 189! Compute function over the locally owned part of the grid 190 191 do j = ys, ye 192 do i = xs, xe 193 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 194 f(i, j) = x(i, j) 195 else 196 u = x(i, j) 197 uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j)) 198 uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1)) 199 f(i, j) = uxx + uyy - sc*exp(u) 200 end if 201 end do 202 end do 203 204 call PetscLogFlops(11.0d0*ym*xm, ierr) 205 CHKERRQ(ierr) 206 207 end 208 209! --------------------------------------------------------------------- 210! 211! FormJacobianLocal - Computes Jacobian matrix, called by 212! the higher level routine FormJacobian(). 213! 214! Input Parameters: 215! x - local vector data 216! 217! Output Parameters: 218! jac - Jacobian matrix 219! jac_prec - optionally different matrix used to construct the preconditioner (not used here) 220! ierr - error code 221! 222! Notes: 223! This routine uses standard Fortran-style computations over a 2-dim array. 224! 225! Notes: 226! Due to grid point reordering with DMDAs, we must always work 227! with the local grid points, and then transform them to the new 228! global numbering with the "ltog" mapping 229! We cannot work directly with the global numbers for the original 230! uniprocessor grid! 231! 232! Two methods are available for imposing this transformation 233! when setting matrix entries: 234! (A) MatSetValuesLocal(), using the local ordering (including 235! ghost points!) 236! by calling MatSetValuesLocal() 237! (B) MatSetValues(), using the global ordering 238! - Use DMDAGetGlobalIndices() to extract the local-to-global map 239! - Then apply this map explicitly yourself 240! - Set matrix entries using the global ordering by calling 241! MatSetValues() 242! Option (A) seems cleaner/easier in many cases, and is the procedure 243! used in this example. 244! 245 subroutine FormJacobianLocal(info, x, A, jac, da, ierr) 246 247 DM da 248 249! Input/output variables: 250 PetscScalar x(gxs:gxe, gys:gye) 251 Mat A, jac 252 PetscErrorCode ierr 253 DMDALocalInfo info 254 255! Local variables: 256 PetscInt row, col(5), i, j, i1, i5 257 PetscScalar two, one, hx, hy, v(5) 258 PetscScalar hxdhy, hydhx, sc 259 260! Set parameters 261 262 i1 = 1 263 i5 = 5 264 one = 1.0 265 two = 2.0 266 hx = one/(real(mx) - 1) 267 hy = one/(real(my) - 1) 268 sc = hx*hy 269 hxdhy = hx/hy 270 hydhx = hy/hx 271! -Wmaybe-uninitialized 272 v = 0.0 273 col = 0 274 275! Compute entries for the locally owned part of the Jacobian. 276! - Currently, all PETSc parallel matrix formats are partitioned by 277! contiguous chunks of rows across the processors. 278! - Each processor needs to insert only elements that it owns 279! locally (but any non-local elements will be sent to the 280! appropriate processor during matrix assembly). 281! - Here, we set all entries for a particular row at once. 282! - We can set matrix entries either using either 283! MatSetValuesLocal() or MatSetValues(), as discussed above. 284! - Note that MatSetValues() uses 0-based row and column numbers 285! in Fortran as well as in C. 286 287 do j = ys, ye 288 row = (j - gys)*gxm + xs - gxs - 1 289 do i = xs, xe 290 row = row + 1 291! boundary points 292 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 293! Some f90 compilers need 4th arg to be of same type in both calls 294 col(1) = row 295 v(1) = one 296 call MatSetValuesLocal(jac, i1, [row], i1, [col], [v], INSERT_VALUES, ierr) 297 CHKERRQ(ierr) 298! interior grid points 299 else 300 v(1) = -hxdhy 301 v(2) = -hydhx 302 v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i, j)) 303 v(4) = -hydhx 304 v(5) = -hxdhy 305 col(1) = row - gxm 306 col(2) = row - 1 307 col(3) = row 308 col(4) = row + 1 309 col(5) = row + gxm 310 call MatSetValuesLocal(jac, i1, [row], i5, [col], [v], INSERT_VALUES, ierr) 311 CHKERRQ(ierr) 312 end if 313 end do 314 end do 315 call MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr) 316 CHKERRQ(ierr) 317 call MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr) 318 CHKERRQ(ierr) 319 if (A /= jac) then 320 call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr) 321 CHKERRQ(ierr) 322 call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr) 323 CHKERRQ(ierr) 324 end if 325 end 326 327! 328! Simple convergence test based on the infinity norm of the residual being small 329! 330 subroutine MySNESConverged(snes, it, xnorm, snorm, fnorm, reason, dummy, ierr) 331 332 SNES snes 333 PetscInt it, dummy 334 PetscReal xnorm, snorm, fnorm, nrm 335 SNESConvergedReason reason 336 Vec f 337 PetscErrorCode ierr 338 339 call SNESGetFunction(snes, f, PETSC_NULL_FUNCTION, dummy, ierr) 340 CHKERRQ(ierr) 341 call VecNorm(f, NORM_INFINITY, nrm, ierr) 342 CHKERRQ(ierr) 343 if (nrm <= 1.e-5) reason = SNES_CONVERGED_FNORM_ABS 344 345 end 346 347end module ex5fmodule 348 349program main 350 use ex5fmodule 351 implicit none 352 353! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 354! Variable declarations 355! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 356! 357! Variables: 358! snes - nonlinear solver 359! x, r - solution, residual vectors 360! its - iterations for convergence 361! 362! See additional variable declarations in the file ex5f.h 363! 364 SNES snes 365 Vec x, r 366 PetscInt its, i1, i4 367 PetscErrorCode ierr 368 PetscReal lambda_max, lambda_min 369 PetscBool flg 370 DM da 371 372! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 373! Initialize program 374! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 375 376 call PetscInitialize(ierr) 377 CHKERRA(ierr) 378 call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr) 379 CHKERRMPIA(ierr) 380 call MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr) 381 CHKERRMPIA(ierr) 382! Initialize problem parameters 383 384 i1 = 1 385 i4 = 4 386 lambda_max = 6.81 387 lambda_min = 0.0 388 lambda = 6.0 389 call PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', lambda, PETSC_NULL_BOOL, ierr) 390 CHKERRA(ierr) 391 392! this statement is split into multiple-lines to keep lines under 132 char limit - required by 'make check' 393 if (lambda >= lambda_max .or. lambda <= lambda_min) then 394 ierr = PETSC_ERR_ARG_OUTOFRANGE 395 SETERRA(PETSC_COMM_WORLD, ierr, 'Lambda') 396 end if 397 398! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 399! Create nonlinear solver context 400! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 401 402 call SNESCreate(PETSC_COMM_WORLD, snes, ierr) 403 CHKERRA(ierr) 404 405! Set convergence test routine if desired 406 407 call PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my_snes_convergence', flg, ierr) 408 CHKERRA(ierr) 409 if (flg) then 410 call SNESSetConvergenceTest(snes, MySNESConverged, 0, PETSC_NULL_FUNCTION, ierr) 411 CHKERRA(ierr) 412 end if 413 414! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 415! Create vector data structures; set function evaluation routine 416! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 417 418! Create distributed array (DMDA) to manage parallel grid and vectors 419 420! This really needs only the star-type stencil, but we use the box stencil 421 422 call DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, i4, i4, PETSC_DECIDE, PETSC_DECIDE, & 423 i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr) 424 CHKERRA(ierr) 425 call DMSetFromOptions(da, ierr) 426 CHKERRA(ierr) 427 call DMSetUp(da, ierr) 428 CHKERRA(ierr) 429 430! Extract global and local vectors from DMDA; then duplicate for remaining 431! vectors that are the same types 432 433 call DMCreateGlobalVector(da, x, ierr) 434 CHKERRA(ierr) 435 call VecDuplicate(x, r, ierr) 436 CHKERRA(ierr) 437 438! Get local grid boundaries (for 2-dimensional DMDA) 439 440 call DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, & 441 PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, & 442 PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr) 443 CHKERRA(ierr) 444 call DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr) 445 CHKERRA(ierr) 446 call DMDAGetGhostCorners(da, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr) 447 CHKERRA(ierr) 448 449! Here we shift the starting indices up by one so that we can easily 450! use the Fortran convention of 1-based indices (rather 0-based indices). 451 452 xs = xs + 1 453 ys = ys + 1 454 gxs = gxs + 1 455 gys = gys + 1 456 457 ye = ys + ym - 1 458 xe = xs + xm - 1 459 gye = gys + gym - 1 460 gxe = gxs + gxm - 1 461 462! Set function evaluation routine and vector 463 464 call DMDASNESSetFunctionLocal(da, INSERT_VALUES, FormFunctionLocal, da, ierr) 465 CHKERRA(ierr) 466 call DMDASNESSetJacobianLocal(da, FormJacobianLocal, da, ierr) 467 CHKERRA(ierr) 468 call SNESSetDM(snes, da, ierr) 469 CHKERRA(ierr) 470 471! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 472! Customize nonlinear solver; set runtime options 473! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 474 475! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 476 477 call SNESSetFromOptions(snes, ierr) 478 CHKERRA(ierr) 479! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 480! Evaluate initial guess; then solve nonlinear system. 481! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 482 483! Note: The user should initialize the vector, x, with the initial guess 484! for the nonlinear solver prior to calling SNESSolve(). In particular, 485! to employ an initial guess of zero, the user should explicitly set 486! this vector to zero by calling VecSet(). 487 488 call FormInitialGuess(x, ierr) 489 CHKERRA(ierr) 490 call SNESSolve(snes, PETSC_NULL_VEC, x, ierr) 491 CHKERRA(ierr) 492 call SNESGetIterationNumber(snes, its, ierr) 493 CHKERRA(ierr) 494 if (rank == 0) then 495 write (6, 100) its 496 end if 497100 format('Number of SNES iterations = ', i5) 498 499! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 500! Free work space. All PETSc objects should be destroyed when they 501! are no longer needed. 502! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 503 504 call VecDestroy(x, ierr) 505 CHKERRA(ierr) 506 call VecDestroy(r, ierr) 507 CHKERRA(ierr) 508 call SNESDestroy(snes, ierr) 509 CHKERRA(ierr) 510 call DMDestroy(da, ierr) 511 CHKERRA(ierr) 512 call PetscFinalize(ierr) 513 CHKERRA(ierr) 514end 515!/*TEST 516! 517! build: 518! requires: !complex !single 519! 520! test: 521! nsize: 4 522! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short \ 523! -ksp_gmres_cgs_refinement_type refine_always 524! 525! test: 526! suffix: 2 527! nsize: 4 528! args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 529! 530! test: 531! suffix: 3 532! nsize: 3 533! args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 534! 535! test: 536! suffix: 6 537! nsize: 1 538! args: -snes_monitor_short -my_snes_convergence 539! 540!TEST*/ 541