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