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