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