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_ARRAY,PETSC_NULL_INTEGER_ARRAY,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_ENUM,PETSC_NULL_ENUM, & 146 PETSC_NULL_ENUM,PETSC_NULL_ENUM,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 end 272 273! --------------------------------------------------------------------- 274! 275! InitialGuessLocal - Computes initial approximation, called by 276! the higher level routine FormInitialGuess(). 277! 278! Input Parameter: 279! x - local vector data 280! 281! Output Parameters: 282! x - local vector data 283! ierr - error code 284! 285! Notes: 286! This routine uses standard Fortran-style computations over a 2-dim array. 287! 288 subroutine InitialGuessLocal(x,ierr) 289 use ex5fmodule 290 implicit none 291 292! Input/output variables: 293 PetscScalar x(xs:xe,ys:ye) 294 PetscErrorCode ierr 295 296! Local variables: 297 PetscInt i,j 298 PetscReal temp1,temp,one,hx,hy 299 300! Set parameters 301 302 ierr = 0 303 one = 1.0 304 hx = one/((real(mx)-1)) 305 hy = one/((real(my)-1)) 306 temp1 = lambda/(lambda + one) 307 308 do 20 j=ys,ye 309 temp = (real(min(j-1,my-j)))*hy 310 do 10 i=xs,xe 311 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 312 x(i,j) = 0.0 313 else 314 x(i,j) = temp1 * sqrt(min(real(min(i-1,mx-i))*hx,(temp))) 315 endif 316 10 continue 317 20 continue 318 319 end 320 321! --------------------------------------------------------------------- 322! 323! FormFunctionLocal - Computes nonlinear function, called by 324! the higher level routine FormFunction(). 325! 326! Input Parameter: 327! x - local vector data 328! 329! Output Parameters: 330! f - local vector data, f(x) 331! ierr - error code 332! 333! Notes: 334! This routine uses standard Fortran-style computations over a 2-dim array. 335! 336! 337 subroutine FormFunctionLocal(info,x,f,da,ierr) 338 use ex5fmodule 339 implicit none 340 341 DM da 342 343! Input/output variables: 344 DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE) 345 PetscScalar x(gxs:gxe,gys:gye) 346 PetscScalar f(xs:xe,ys:ye) 347 PetscErrorCode ierr 348 349! Local variables: 350 PetscScalar two,one,hx,hy 351 PetscScalar hxdhy,hydhx,sc 352 PetscScalar u,uxx,uyy 353 PetscInt i,j 354 355 xs = info(DMDA_LOCAL_INFO_XS)+1 356 xe = xs+info(DMDA_LOCAL_INFO_XM)-1 357 ys = info(DMDA_LOCAL_INFO_YS)+1 358 ye = ys+info(DMDA_LOCAL_INFO_YM)-1 359 mx = info(DMDA_LOCAL_INFO_MX) 360 my = info(DMDA_LOCAL_INFO_MY) 361 362 one = 1.0 363 two = 2.0 364 hx = one/(real(mx)-1) 365 hy = one/(real(my)-1) 366 sc = hx*hy*lambda 367 hxdhy = hx/hy 368 hydhx = hy/hx 369 370! Compute function over the locally owned part of the grid 371 372 do 20 j=ys,ye 373 do 10 i=xs,xe 374 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 375 f(i,j) = x(i,j) 376 else 377 u = x(i,j) 378 uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j)) 379 uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1)) 380 f(i,j) = uxx + uyy - sc*exp(u) 381 endif 382 10 continue 383 20 continue 384 385 call PetscLogFlops(11.0d0*ym*xm,ierr) 386 CHKERRQ(ierr) 387 388 end 389 390! --------------------------------------------------------------------- 391! 392! FormJacobianLocal - Computes Jacobian matrix, called by 393! the higher level routine FormJacobian(). 394! 395! Input Parameters: 396! x - local vector data 397! 398! Output Parameters: 399! jac - Jacobian matrix 400! jac_prec - optionally different preconditioning matrix (not used here) 401! ierr - error code 402! 403! Notes: 404! This routine uses standard Fortran-style computations over a 2-dim array. 405! 406! Notes: 407! Due to grid point reordering with DMDAs, we must always work 408! with the local grid points, and then transform them to the new 409! global numbering with the "ltog" mapping 410! We cannot work directly with the global numbers for the original 411! uniprocessor grid! 412! 413! Two methods are available for imposing this transformation 414! when setting matrix entries: 415! (A) MatSetValuesLocal(), using the local ordering (including 416! ghost points!) 417! by calling MatSetValuesLocal() 418! (B) MatSetValues(), using the global ordering 419! - Use DMDAGetGlobalIndices() to extract the local-to-global map 420! - Then apply this map explicitly yourself 421! - Set matrix entries using the global ordering by calling 422! MatSetValues() 423! Option (A) seems cleaner/easier in many cases, and is the procedure 424! used in this example. 425! 426 subroutine FormJacobianLocal(info,x,A,jac,da,ierr) 427 use ex5fmodule 428 implicit none 429 430 DM da 431 432! Input/output variables: 433 PetscScalar x(gxs:gxe,gys:gye) 434 Mat A,jac 435 PetscErrorCode ierr 436 DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE) 437 438! Local variables: 439 PetscInt row,col(5),i,j,i1,i5 440 PetscScalar two,one,hx,hy,v(5) 441 PetscScalar hxdhy,hydhx,sc 442 443! Set parameters 444 445 i1 = 1 446 i5 = 5 447 one = 1.0 448 two = 2.0 449 hx = one/(real(mx)-1) 450 hy = one/(real(my)-1) 451 sc = hx*hy 452 hxdhy = hx/hy 453 hydhx = hy/hx 454! -Wmaybe-uninitialized 455 v = 0.0 456 col = 0 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 end 509 510! 511! Simple convergence test based on the infinity norm of the residual being small 512! 513 subroutine MySNESConverged(snes,it,xnorm,snorm,fnorm,reason,dummy,ierr) 514 use ex5fmodule 515 implicit none 516 517 SNES snes 518 PetscInt it,dummy 519 PetscReal xnorm,snorm,fnorm,nrm 520 SNESConvergedReason reason 521 Vec f 522 PetscErrorCode ierr 523 524 call SNESGetFunction(snes,f,PETSC_NULL_FUNCTION,dummy,ierr) 525 CHKERRQ(ierr) 526 call VecNorm(f,NORM_INFINITY,nrm,ierr) 527 CHKERRQ(ierr) 528 if (nrm .le. 1.e-5) reason = SNES_CONVERGED_FNORM_ABS 529 530 end 531 532!/*TEST 533! 534! build: 535! requires: !complex !single 536! 537! test: 538! nsize: 4 539! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short \ 540! -ksp_gmres_cgs_refinement_type refine_always 541! 542! test: 543! suffix: 2 544! nsize: 4 545! args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 546! 547! test: 548! suffix: 3 549! nsize: 3 550! args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 551! 552! test: 553! suffix: 6 554! nsize: 1 555! args: -snes_monitor_short -my_snes_convergence 556! 557!TEST*/ 558