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