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