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