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