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