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