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_ARRAY,PETSC_NULL_INTEGER_ARRAY,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_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,PETSC_NULL_ENUM,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! 484! Notes: 485! This routine serves as a wrapper for the lower-level routine 486! "FormJacobianLocal", where the actual computations are 487! done using the standard Fortran style of treating the local 488! vector data as a multidimensional array over the local mesh. 489! This routine merely accesses the local vector data via 490! VecGetArrayF90() and VecRestoreArrayF90(). 491! 492! Notes: 493! Due to grid point reordering with DMDAs, we must always work 494! with the local grid points, and then transform them to the new 495! global numbering with the "ltog" mapping 496! We cannot work directly with the global numbers for the original 497! uniprocessor grid! 498! 499! Two methods are available for imposing this transformation 500! when setting matrix entries: 501! (A) MatSetValuesLocal(), using the local ordering (including 502! ghost points!) 503! - Set matrix entries using the local ordering 504! by calling MatSetValuesLocal() 505! (B) MatSetValues(), using the global ordering 506! - Use DMGetLocalToGlobalMapping() then 507! ISLocalToGlobalMappingGetIndicesF90() to extract the local-to-global map 508! - Then apply this map explicitly yourself 509! - Set matrix entries using the global ordering by calling 510! MatSetValues() 511! Option (A) seems cleaner/easier in many cases, and is the procedure 512! used in this example. 513! 514 subroutine FormJacobian(mysnes,X,jac,jac_prec,user,ierr) 515#include <petsc/finclude/petscsnes.h> 516 use petscsnes 517 use ex5f90tmodule 518! Input/output variables: 519 type(tSNES) mysnes 520 type(tVec) X 521 type(tMat) jac,jac_prec 522 type(userctx) user 523 PetscErrorCode ierr 524 525! Declarations for use with local arrays: 526 PetscScalar,pointer :: lx_v(:) 527 type(tVec) localX 528 529! Scatter ghost points to local vector, using the 2-step process 530! DMGlobalToLocalBegin(), DMGlobalToLocalEnd() 531! Computations can be done while messages are in transition, 532! by placing code between these two statements. 533 534 PetscCallA(DMGetLocalVector(user%da,localX,ierr)) 535 PetscCallA(DMGlobalToLocalBegin(user%da,X,INSERT_VALUES,localX,ierr)) 536 PetscCallA(DMGlobalToLocalEnd(user%da,X,INSERT_VALUES,localX,ierr)) 537 538! Get a pointer to vector data 539 PetscCallA(VecGetArrayF90(localX,lx_v,ierr)) 540 541! Compute entries for the locally owned part of the Jacobian preconditioner. 542 PetscCallA(FormJacobianLocal(lx_v,jac_prec,user,ierr)) 543 544! Assemble matrix, using the 2-step process: 545! MatAssemblyBegin(), MatAssemblyEnd() 546! Computations can be done while messages are in transition, 547! by placing code between these two statements. 548 549 PetscCallA(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 550! if (jac .ne. jac_prec) then 551 PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 552! endif 553 PetscCallA(VecRestoreArrayF90(localX,lx_v,ierr)) 554 PetscCallA(DMRestoreLocalVector(user%da,localX,ierr)) 555 PetscCallA(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 556! if (jac .ne. jac_prec) then 557 PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 558! endif 559 560! Tell the matrix we will never add a new nonzero location to the 561! matrix. If we do it will generate an error. 562 563 PetscCallA(MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 564 565 end 566 567! --------------------------------------------------------------------- 568! 569! FormJacobianLocal - Computes Jacobian preconditioner matrix, 570! called by the higher level routine FormJacobian(). 571! 572! Input Parameters: 573! x - local vector data 574! 575! Output Parameters: 576! jac_prec - Jacobian preconditioner matrix 577! ierr - error code 578! 579! Notes: 580! This routine uses standard Fortran-style computations over a 2-dim array. 581! 582! Notes: 583! Due to grid point reordering with DMDAs, we must always work 584! with the local grid points, and then transform them to the new 585! global numbering with the "ltog" mapping 586! We cannot work directly with the global numbers for the original 587! uniprocessor grid! 588! 589! Two methods are available for imposing this transformation 590! when setting matrix entries: 591! (A) MatSetValuesLocal(), using the local ordering (including 592! ghost points!) 593! - Set matrix entries using the local ordering 594! by calling MatSetValuesLocal() 595! (B) MatSetValues(), using the global ordering 596! - Set matrix entries using the global ordering by calling 597! MatSetValues() 598! Option (A) seems cleaner/easier in many cases, and is the procedure 599! used in this example. 600! 601 subroutine FormJacobianLocal(x,jac_prec,user,ierr) 602#include <petsc/finclude/petscmat.h> 603 use petscmat 604 use ex5f90tmodule 605! Input/output variables: 606 type (userctx) user 607 PetscScalar x(user%gxs:user%gxe,user%gys:user%gye) 608 type(tMat) jac_prec 609 PetscErrorCode ierr 610 611! Local variables: 612 PetscInt row,col(5),i,j 613 PetscInt ione,ifive 614 PetscScalar two,one,hx,hy,hxdhy 615 PetscScalar hydhx,sc,v(5) 616 617! Set parameters 618 ione = 1 619 ifive = 5 620 one = 1.0 621 two = 2.0 622 hx = one/PetscIntToReal(user%mx-1) 623 hy = one/PetscIntToReal(user%my-1) 624 sc = hx*hy 625 hxdhy = hx/hy 626 hydhx = hy/hx 627 628! Compute entries for the locally owned part of the Jacobian. 629! - Currently, all PETSc parallel matrix formats are partitioned by 630! contiguous chunks of rows across the processors. 631! - Each processor needs to insert only elements that it owns 632! locally (but any non-local elements will be sent to the 633! appropriate processor during matrix assembly). 634! - Here, we set all entries for a particular row at once. 635! - We can set matrix entries either using either 636! MatSetValuesLocal() or MatSetValues(), as discussed above. 637! - Note that MatSetValues() uses 0-based row and column numbers 638! in Fortran as well as in C. 639 640 do 20 j=user%ys,user%ye 641 row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1 642 do 10 i=user%xs,user%xe 643 row = row + 1 644! boundary points 645 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. user%mx .or. j .eq. user%my) then 646 col(1) = row 647 v(1) = one 648 PetscCallA(MatSetValuesLocal(jac_prec,ione,[row],ione,col,v,INSERT_VALUES,ierr)) 649! interior grid points 650 else 651 v(1) = -hxdhy 652 v(2) = -hydhx 653 v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i,j)) 654 v(4) = -hydhx 655 v(5) = -hxdhy 656 col(1) = row - user%gxm 657 col(2) = row - 1 658 col(3) = row 659 col(4) = row + 1 660 col(5) = row + user%gxm 661 PetscCallA(MatSetValuesLocal(jac_prec,ione,[row],ifive,col,v,INSERT_VALUES,ierr)) 662 endif 663 10 continue 664 20 continue 665 end 666 667!/*TEST 668! 669! test: 670! nsize: 4 671! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 672! 673!TEST*/ 674