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