1! 2! Description: 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 <parameter>, where <parameter> indicates the nonlinearity of the problem 7! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 8! 9! 10! -------------------------------------------------------------------------- 11! 12! Solid Fuel Ignition (SFI) problem. This problem is modeled by 13! the partial differential equation 14! 15! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 16! 17! with boundary conditions 18! 19! u = 0 for x = 0, x = 1, y = 0, y = 1. 20! 21! A finite difference approximation with the usual 5-point stencil 22! is used to discretize the boundary value problem to obtain a nonlinear 23! system of equations. 24! 25! The uniprocessor version of this code is snes/tutorials/ex4f.F 26! 27! -------------------------------------------------------------------------- 28! The following define must be used before including any PETSc include files 29! into a module or interface. This is because they can't handle declarations 30! in them 31! 32 33 module f90modulet 34#include <petsc/finclude/petscdm.h> 35 use petscdmdef 36 type userctx 37 type(tDM) da 38 PetscInt xs,xe,xm,gxs,gxe,gxm 39 PetscInt ys,ye,ym,gys,gye,gym 40 PetscInt mx,my 41 PetscMPIInt rank 42 PetscReal lambda 43 end type userctx 44 45 contains 46! --------------------------------------------------------------------- 47! 48! FormFunction - Evaluates nonlinear function, F(x). 49! 50! Input Parameters: 51! snes - the SNES context 52! X - input vector 53! dummy - optional user-defined context, as set by SNESSetFunction() 54! (not used here) 55! 56! Output Parameter: 57! F - function vector 58! 59! Notes: 60! This routine serves as a wrapper for the lower-level routine 61! "FormFunctionLocal", where the actual computations are 62! done using the standard Fortran style of treating the local 63! vector data as a multidimensional array over the local mesh. 64! This routine merely handles ghost point scatters and accesses 65! the local vector data via VecGetArrayF90() and VecRestoreArrayF90(). 66! 67 subroutine FormFunction(snesIn,X,F,user,ierr) 68#include <petsc/finclude/petscsnes.h> 69 use petscsnes 70 71! Input/output variables: 72 type(tSNES) snesIn 73 type(tVec) X,F 74 PetscErrorCode ierr 75 type (userctx) user 76 77! Declarations for use with local arrays: 78 PetscScalar,pointer :: lx_v(:),lf_v(:) 79 type(tVec) localX 80 81! Scatter ghost points to local vector, using the 2-step process 82! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 83! By placing code between these two statements, computations can 84! be done while messages are in transition. 85 call DMGetLocalVector(user%da,localX,ierr);CHKERRQ(ierr) 86 call DMGlobalToLocalBegin(user%da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr) 87 call DMGlobalToLocalEnd(user%da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr) 88 89! Get a pointer to vector data. 90! - For default PETSc vectors, VecGetArray90() returns a pointer to 91! the data array. Otherwise, the routine is implementation dependent. 92! - You MUST call VecRestoreArrayF90() when you no longer need access to 93! the array. 94! - Note that the interface to VecGetArrayF90() differs from VecGetArray(), 95! and is useable from Fortran-90 Only. 96 97 call VecGetArrayF90(localX,lx_v,ierr);CHKERRQ(ierr) 98 call VecGetArrayF90(F,lf_v,ierr);CHKERRQ(ierr) 99 100! Compute function over the locally owned part of the grid 101 call FormFunctionLocal(lx_v,lf_v,user,ierr);CHKERRQ(ierr) 102 103! Restore vectors 104 call VecRestoreArrayF90(localX,lx_v,ierr);CHKERRQ(ierr) 105 call VecRestoreArrayF90(F,lf_v,ierr);CHKERRQ(ierr) 106 107! Insert values into global vector 108 109 call DMRestoreLocalVector(user%da,localX,ierr);CHKERRQ(ierr) 110 call PetscLogFlops(11.0d0*user%ym*user%xm,ierr) 111 112! call VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr) 113! call VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr) 114 return 115 end subroutine formfunction 116 end module f90modulet 117 118 module f90moduleinterfacest 119 use f90modulet 120 121 Interface SNESSetApplicationContext 122 Subroutine SNESSetApplicationContext(snesIn,ctx,ierr) 123#include <petsc/finclude/petscsnes.h> 124 use petscsnes 125 use f90modulet 126 type(tSNES) snesIn 127 type(userctx) ctx 128 PetscErrorCode ierr 129 End Subroutine 130 End Interface SNESSetApplicationContext 131 132 Interface SNESGetApplicationContext 133 Subroutine SNESGetApplicationContext(snesIn,ctx,ierr) 134#include <petsc/finclude/petscsnes.h> 135 use petscsnes 136 use f90modulet 137 type(tSNES) snesIn 138 type(userctx), pointer :: ctx 139 PetscErrorCode ierr 140 End Subroutine 141 End Interface SNESGetApplicationContext 142 end module f90moduleinterfacest 143 144 program main 145#include <petsc/finclude/petscdm.h> 146#include <petsc/finclude/petscsnes.h> 147 use petscdmda 148 use petscdm 149 use petscsnes 150 use f90modulet 151 use f90moduleinterfacest 152 implicit none 153! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 154! Variable declarations 155! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156! 157! Variables: 158! mysnes - nonlinear solver 159! x, r - solution, residual vectors 160! J - Jacobian matrix 161! its - iterations for convergence 162! Nx, Ny - number of preocessors in x- and y- directions 163! matrix_free - flag - 1 indicates matrix-free version 164! 165 type(tSNES) mysnes 166 type(tVec) x,r 167 type(tMat) J 168 PetscErrorCode ierr 169 PetscInt its 170 PetscBool flg,matrix_free,set 171 PetscInt ione,nfour 172 PetscReal lambda_max,lambda_min 173 type(userctx) user 174 type(userctx), pointer:: puser 175 type(tPetscOptions) :: options 176 177! Note: Any user-defined Fortran routines (such as FormJacobian) 178! MUST be declared as external. 179 external FormInitialGuess,FormJacobian 180 181! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 182! Initialize program 183! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 185 if (ierr .ne. 0) then 186 print*,'Unable to initialize PETSc' 187 stop 188 endif 189 call MPI_Comm_rank(PETSC_COMM_WORLD,user%rank,ierr) 190 191! Initialize problem parameters 192 options%v = 0 193 lambda_max = 6.81 194 lambda_min = 0.0 195 user%lambda = 6.0 196 ione = 1 197 nfour = 4 198 call PetscOptionsGetReal(options,PETSC_NULL_CHARACTER,'-par',user%lambda,flg,ierr);CHKERRA(ierr) 199 if (user%lambda .ge. lambda_max .or. user%lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range '); endif 200 201! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202! Create nonlinear solver context 203! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 204 call SNESCreate(PETSC_COMM_WORLD,mysnes,ierr);CHKERRA(ierr) 205 206! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207! Create vector data structures; set function evaluation routine 208! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 209 210! Create distributed array (DMDA) to manage parallel grid and vectors 211 212! This really needs only the star-type stencil, but we use the box 213! stencil temporarily. 214 call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione, & 215 & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,user%da,ierr);CHKERRA(ierr) 216 call DMSetFromOptions(user%da,ierr);CHKERRA(ierr) 217 call DMSetUp(user%da,ierr);CHKERRA(ierr) 218 call DMDAGetInfo(user%da,PETSC_NULL_INTEGER,user%mx,user%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 219 & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr) 220 221! 222! Visualize the distribution of the array across the processors 223! 224! call DMView(user%da,PETSC_VIEWER_DRAW_WORLD,ierr) 225 226! Extract global and local vectors from DMDA; then duplicate for remaining 227! vectors that are the same types 228 call DMCreateGlobalVector(user%da,x,ierr);CHKERRA(ierr) 229 call VecDuplicate(x,r,ierr);CHKERRA(ierr) 230 231! Get local grid boundaries (for 2-dimensional DMDA) 232 call DMDAGetCorners(user%da,user%xs,user%ys,PETSC_NULL_INTEGER,user%xm,user%ym,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr) 233 call DMDAGetGhostCorners(user%da,user%gxs,user%gys,PETSC_NULL_INTEGER,user%gxm,user%gym,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr) 234 235! Here we shift the starting indices up by one so that we can easily 236! use the Fortran convention of 1-based indices (rather 0-based indices). 237 user%xs = user%xs+1 238 user%ys = user%ys+1 239 user%gxs = user%gxs+1 240 user%gys = user%gys+1 241 242 user%ye = user%ys+user%ym-1 243 user%xe = user%xs+user%xm-1 244 user%gye = user%gys+user%gym-1 245 user%gxe = user%gxs+user%gxm-1 246 247 call SNESSetApplicationContext(mysnes,user,ierr);CHKERRA(ierr) 248 249! Set function evaluation routine and vector 250 call SNESSetFunction(mysnes,r,FormFunction,user,ierr);CHKERRA(ierr) 251 252! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 253! Create matrix data structure; set Jacobian evaluation routine 254! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 255 256! Set Jacobian matrix data structure and default Jacobian evaluation 257! routine. User can override with: 258! -snes_fd : default finite differencing approximation of Jacobian 259! -snes_mf : matrix-free Newton-Krylov method with no preconditioning 260! (unless user explicitly sets preconditioner) 261! -snes_mf_operator : form preconditioning matrix as set by the user, 262! but use matrix-free approx for Jacobian-vector 263! products within Newton-Krylov method 264! 265! Note: For the parallel case, vectors and matrices MUST be partitioned 266! accordingly. When using distributed arrays (DMDAs) to create vectors, 267! the DMDAs determine the problem partitioning. We must explicitly 268! specify the local matrix dimensions upon its creation for compatibility 269! with the vector distribution. Thus, the generic MatCreate() routine 270! is NOT sufficient when working with distributed arrays. 271! 272! Note: Here we only approximately preallocate storage space for the 273! Jacobian. See the users manual for a discussion of better techniques 274! for preallocating matrix memory. 275 276 call PetscOptionsHasName(options,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr);CHKERRA(ierr) 277 if (.not. matrix_free) then 278 call DMSetMatType(user%da,MATAIJ,ierr);CHKERRA(ierr) 279 call DMCreateMatrix(user%da,J,ierr);CHKERRA(ierr) 280 call SNESSetJacobian(mysnes,J,J,FormJacobian,user,ierr);CHKERRA(ierr) 281 endif 282 283! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 284! Customize nonlinear solver; set runtime options 285! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 286! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 287 call SNESSetFromOptions(mysnes,ierr);CHKERRA(ierr) 288 289! Test Fortran90 wrapper for SNESSet/Get ApplicationContext() 290 call PetscOptionsGetBool(options,PETSC_NULL_CHARACTER,'-test_appctx',flg,set,ierr);CHKERRA(ierr) 291 if (flg) then 292 call SNESGetApplicationContext(mysnes,puser,ierr);CHKERRA(ierr) 293 endif 294 295! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 296! Evaluate initial guess; then solve nonlinear system. 297! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 298! Note: The user should initialize the vector, x, with the initial guess 299! for the nonlinear solver prior to calling SNESSolve(). In particular, 300! to employ an initial guess of zero, the user should explicitly set 301! this vector to zero by calling VecSet(). 302 303 call FormInitialGuess(mysnes,x,ierr);CHKERRA(ierr) 304 call SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr);CHKERRA(ierr) 305 call SNESGetIterationNumber(mysnes,its,ierr);CHKERRA(ierr) 306 if (user%rank .eq. 0) then 307 write(6,100) its 308 endif 309 100 format('Number of SNES iterations = ',i5) 310 311! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 312! Free work space. All PETSc objects should be destroyed when they 313! are no longer needed. 314! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 315 if (.not. matrix_free) call MatDestroy(J,ierr);CHKERRA(ierr) 316 call VecDestroy(x,ierr);CHKERRA(ierr) 317 call VecDestroy(r,ierr);CHKERRA(ierr) 318 call SNESDestroy(mysnes,ierr);CHKERRA(ierr) 319 call DMDestroy(user%da,ierr);CHKERRA(ierr) 320 321 call PetscFinalize(ierr) 322 end 323 324! --------------------------------------------------------------------- 325! 326! FormInitialGuess - Forms initial approximation. 327! 328! Input Parameters: 329! X - vector 330! 331! Output Parameter: 332! X - vector 333! 334! Notes: 335! This routine serves as a wrapper for the lower-level routine 336! "InitialGuessLocal", where the actual computations are 337! done using the standard Fortran style of treating the local 338! vector data as a multidimensional array over the local mesh. 339! This routine merely handles ghost point scatters and accesses 340! the local vector data via VecGetArrayF90() and VecRestoreArrayF90(). 341! 342 subroutine FormInitialGuess(mysnes,X,ierr) 343#include <petsc/finclude/petscsnes.h> 344 use petscsnes 345 use f90modulet 346 use f90moduleinterfacest 347! Input/output variables: 348 type(tSNES) mysnes 349 type(userctx), pointer:: puser 350 type(tVec) X 351 PetscErrorCode ierr 352 353! Declarations for use with local arrays: 354 PetscScalar,pointer :: lx_v(:) 355 356 ierr = 0 357 call SNESGetApplicationContext(mysnes,puser,ierr) 358! Get a pointer to vector data. 359! - For default PETSc vectors, VecGetArray90() returns a pointer to 360! the data array. Otherwise, the routine is implementation dependent. 361! - You MUST call VecRestoreArrayF90() when you no longer need access to 362! the array. 363! - Note that the interface to VecGetArrayF90() differs from VecGetArray(), 364! and is useable from Fortran-90 Only. 365 366 call VecGetArrayF90(X,lx_v,ierr) 367 368! Compute initial guess over the locally owned part of the grid 369 call InitialGuessLocal(puser,lx_v,ierr) 370 371! Restore vector 372 call VecRestoreArrayF90(X,lx_v,ierr) 373 374! Insert values into global vector 375 376 return 377 end 378 379! --------------------------------------------------------------------- 380! 381! InitialGuessLocal - Computes initial approximation, called by 382! the higher level routine FormInitialGuess(). 383! 384! Input Parameter: 385! x - local vector data 386! 387! Output Parameters: 388! x - local vector data 389! ierr - error code 390! 391! Notes: 392! This routine uses standard Fortran-style computations over a 2-dim array. 393! 394 subroutine InitialGuessLocal(user,x,ierr) 395#include <petsc/finclude/petscsys.h> 396 use petscsys 397 use f90modulet 398! Input/output variables: 399 type (userctx) user 400 PetscScalar x(user%xs:user%xe,user%ys:user%ye) 401 PetscErrorCode ierr 402 403! Local variables: 404 PetscInt i,j 405 PetscScalar temp1,temp,hx,hy 406 PetscScalar one 407 408! Set parameters 409 410 ierr = 0 411 one = 1.0 412 hx = one/(PetscIntToReal(user%mx-1)) 413 hy = one/(PetscIntToReal(user%my-1)) 414 temp1 = user%lambda/(user%lambda + one) 415 416 do 20 j=user%ys,user%ye 417 temp = PetscIntToReal(min(j-1,user%my-j))*hy 418 do 10 i=user%xs,user%xe 419 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then 420 x(i,j) = 0.0 421 else 422 x(i,j) = temp1 * sqrt(min(PetscIntToReal(min(i-1,user%mx-i)*hx),PetscIntToReal(temp))) 423 endif 424 10 continue 425 20 continue 426 427 return 428 end 429 430! --------------------------------------------------------------------- 431! 432! FormFunctionLocal - Computes nonlinear function, called by 433! the higher level routine FormFunction(). 434! 435! Input Parameter: 436! x - local vector data 437! 438! Output Parameters: 439! f - local vector data, f(x) 440! ierr - error code 441! 442! Notes: 443! This routine uses standard Fortran-style computations over a 2-dim array. 444! 445 subroutine FormFunctionLocal(x,f,user,ierr) 446#include <petsc/finclude/petscsys.h> 447 use petscsys 448 use f90modulet 449! Input/output variables: 450 type (userctx) user 451 PetscScalar x(user%gxs:user%gxe,user%gys:user%gye) 452 PetscScalar f(user%xs:user%xe,user%ys:user%ye) 453 PetscErrorCode ierr 454 455! Local variables: 456 PetscScalar two,one,hx,hy,hxdhy,hydhx,sc 457 PetscScalar u,uxx,uyy 458 PetscInt i,j 459 460 one = 1.0 461 two = 2.0 462 hx = one/PetscIntToReal(user%mx-1) 463 hy = one/PetscIntToReal(user%my-1) 464 sc = hx*hy*user%lambda 465 hxdhy = hx/hy 466 hydhx = hy/hx 467 468! Compute function over the locally owned part of the grid 469 470 do 20 j=user%ys,user%ye 471 do 10 i=user%xs,user%xe 472 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then 473 f(i,j) = x(i,j) 474 else 475 u = x(i,j) 476 uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j)) 477 uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1)) 478 f(i,j) = uxx + uyy - sc*exp(u) 479 endif 480 10 continue 481 20 continue 482 ierr = 0 483 return 484 end 485 486! --------------------------------------------------------------------- 487! 488! FormJacobian - Evaluates Jacobian matrix. 489! 490! Input Parameters: 491! snes - the SNES context 492! x - input vector 493! dummy - optional user-defined context, as set by SNESSetJacobian() 494! (not used here) 495! 496! Output Parameters: 497! jac - Jacobian matrix 498! jac_prec - optionally different preconditioning matrix (not used here) 499! flag - flag indicating matrix structure 500! 501! Notes: 502! This routine serves as a wrapper for the lower-level routine 503! "FormJacobianLocal", where the actual computations are 504! done using the standard Fortran style of treating the local 505! vector data as a multidimensional array over the local mesh. 506! This routine merely accesses the local vector data via 507! VecGetArrayF90() and VecRestoreArrayF90(). 508! 509! Notes: 510! Due to grid point reordering with DMDAs, we must always work 511! with the local grid points, and then transform them to the new 512! global numbering with the "ltog" mapping 513! We cannot work directly with the global numbers for the original 514! uniprocessor grid! 515! 516! Two methods are available for imposing this transformation 517! when setting matrix entries: 518! (A) MatSetValuesLocal(), using the local ordering (including 519! ghost points!) 520! - Set matrix entries using the local ordering 521! by calling MatSetValuesLocal() 522! (B) MatSetValues(), using the global ordering 523! - Use DMGetLocalToGlobalMapping() then 524! ISLocalToGlobalMappingGetIndicesF90() to extract the local-to-global map 525! - Then apply this map explicitly yourself 526! - Set matrix entries using the global ordering by calling 527! MatSetValues() 528! Option (A) seems cleaner/easier in many cases, and is the procedure 529! used in this example. 530! 531 subroutine FormJacobian(mysnes,X,jac,jac_prec,user,ierr) 532#include <petsc/finclude/petscsnes.h> 533 use petscsnes 534 use f90modulet 535! Input/output variables: 536 type(tSNES) mysnes 537 type(tVec) X 538 type(tMat) jac,jac_prec 539 type(userctx) user 540 PetscErrorCode ierr 541 542! Declarations for use with local arrays: 543 PetscScalar,pointer :: lx_v(:) 544 type(tVec) localX 545 546! Scatter ghost points to local vector, using the 2-step process 547! DMGlobalToLocalBegin(), DMGlobalToLocalEnd() 548! Computations can be done while messages are in transition, 549! by placing code between these two statements. 550 551 call DMGetLocalVector(user%da,localX,ierr) 552 call DMGlobalToLocalBegin(user%da,X,INSERT_VALUES,localX,ierr) 553 call DMGlobalToLocalEnd(user%da,X,INSERT_VALUES,localX,ierr) 554 555! Get a pointer to vector data 556 call VecGetArrayF90(localX,lx_v,ierr) 557 558! Compute entries for the locally owned part of the Jacobian preconditioner. 559 call FormJacobianLocal(lx_v,jac_prec,user,ierr) 560 561! Assemble matrix, using the 2-step process: 562! MatAssemblyBegin(), MatAssemblyEnd() 563! Computations can be done while messages are in transition, 564! by placing code between these two statements. 565 566 call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr) 567! if (jac .ne. jac_prec) then 568 call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr) 569! endif 570 call VecRestoreArrayF90(localX,lx_v,ierr) 571 call DMRestoreLocalVector(user%da,localX,ierr) 572 call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr) 573! if (jac .ne. jac_prec) then 574 call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr) 575! endif 576 577! Tell the matrix we will never add a new nonzero location to the 578! matrix. If we do it will generate an error. 579 580 call MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr) 581 582 return 583 end 584 585! --------------------------------------------------------------------- 586! 587! FormJacobianLocal - Computes Jacobian preconditioner matrix, 588! called by the higher level routine FormJacobian(). 589! 590! Input Parameters: 591! x - local vector data 592! 593! Output Parameters: 594! jac_prec - Jacobian preconditioner matrix 595! ierr - error code 596! 597! Notes: 598! This routine uses standard Fortran-style computations over a 2-dim array. 599! 600! Notes: 601! Due to grid point reordering with DMDAs, we must always work 602! with the local grid points, and then transform them to the new 603! global numbering with the "ltog" mapping 604! We cannot work directly with the global numbers for the original 605! uniprocessor grid! 606! 607! Two methods are available for imposing this transformation 608! when setting matrix entries: 609! (A) MatSetValuesLocal(), using the local ordering (including 610! ghost points!) 611! - Set matrix entries using the local ordering 612! by calling MatSetValuesLocal() 613! (B) MatSetValues(), using the global ordering 614! - Set matrix entries using the global ordering by calling 615! MatSetValues() 616! Option (A) seems cleaner/easier in many cases, and is the procedure 617! used in this example. 618! 619 subroutine FormJacobianLocal(x,jac_prec,user,ierr) 620#include <petsc/finclude/petscmat.h> 621 use petscmat 622 use f90modulet 623! Input/output variables: 624 type (userctx) user 625 PetscScalar x(user%gxs:user%gxe,user%gys:user%gye) 626 type(tMat) jac_prec 627 PetscErrorCode ierr 628 629! Local variables: 630 PetscInt row,col(5),i,j 631 PetscInt ione,ifive 632 PetscScalar two,one,hx,hy,hxdhy 633 PetscScalar hydhx,sc,v(5) 634 635! Set parameters 636 ione = 1 637 ifive = 5 638 one = 1.0 639 two = 2.0 640 hx = one/PetscIntToReal(user%mx-1) 641 hy = one/PetscIntToReal(user%my-1) 642 sc = hx*hy 643 hxdhy = hx/hy 644 hydhx = hy/hx 645 646! Compute entries for the locally owned part of the Jacobian. 647! - Currently, all PETSc parallel matrix formats are partitioned by 648! contiguous chunks of rows across the processors. 649! - Each processor needs to insert only elements that it owns 650! locally (but any non-local elements will be sent to the 651! appropriate processor during matrix assembly). 652! - Here, we set all entries for a particular row at once. 653! - We can set matrix entries either using either 654! MatSetValuesLocal() or MatSetValues(), as discussed above. 655! - Note that MatSetValues() uses 0-based row and column numbers 656! in Fortran as well as in C. 657 658 do 20 j=user%ys,user%ye 659 row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1 660 do 10 i=user%xs,user%xe 661 row = row + 1 662! boundary points 663 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then 664 col(1) = row 665 v(1) = one 666 call MatSetValuesLocal(jac_prec,ione,row,ione,col,v,INSERT_VALUES,ierr) 667! interior grid points 668 else 669 v(1) = -hxdhy 670 v(2) = -hydhx 671 v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i,j)) 672 v(4) = -hydhx 673 v(5) = -hxdhy 674 col(1) = row - user%gxm 675 col(2) = row - 1 676 col(3) = row 677 col(4) = row + 1 678 col(5) = row + user%gxm 679 call MatSetValuesLocal(jac_prec,ione,row,ifive,col,v,INSERT_VALUES,ierr) 680 endif 681 10 continue 682 20 continue 683 return 684 end 685 686!/*TEST 687! 688! test: 689! nsize: 4 690! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 691! 692!TEST*/ 693