1! 2! Description: This example solves a nonlinear system on 1 processor with SNES. 3! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular 4! domain. The command line options include: 5! -par <parameter>, where <parameter> indicates the nonlinearity of the problem 6! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 7! -mx <xg>, where <xg> = number of grid points in the x-direction 8! -my <yg>, where <yg> = number of grid points in the y-direction 9! 10 11! 12! -------------------------------------------------------------------------- 13! 14! Solid Fuel Ignition (SFI) problem. This problem is modeled by 15! the partial differential equation 16! 17! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 18! 19! with boundary conditions 20! 21! u = 0 for x = 0, x = 1, y = 0, y = 1. 22! 23! A finite difference approximation with the usual 5-point stencil 24! is used to discretize the boundary value problem to obtain a nonlinear 25! system of equations. 26! 27! The parallel version of this code is snes/tutorials/ex5f.F 28! 29! -------------------------------------------------------------------------- 30 subroutine postcheck(snes,x,y,w,changed_y,changed_w,ctx,ierr) 31#include <petsc/finclude/petscsnes.h> 32 use petscsnes 33 implicit none 34 SNES snes 35 PetscReal norm 36 Vec tmp,x,y,w 37 PetscBool changed_w,changed_y 38 PetscErrorCode ierr 39 PetscInt ctx 40 PetscScalar mone 41 42 PetscCallA(VecDuplicate(x,tmp,ierr)) 43 mone = -1.0 44 PetscCallA(VecWAXPY(tmp,mone,x,w,ierr)) 45 PetscCallA(VecNorm(tmp,NORM_2,norm,ierr)) 46 PetscCallA(VecDestroy(tmp,ierr)) 47 print*, 'Norm of search step ',norm 48 changed_y = PETSC_FALSE 49 changed_w = PETSC_FALSE 50 return 51 end 52 53 program main 54#include <petsc/finclude/petscdraw.h> 55 use petscsnes 56 implicit none 57 interface SNESSetJacobian 58 subroutine SNESSetJacobian1(a,b,c,d,e,z) 59 use petscsnes 60 SNES a 61 Mat b 62 Mat c 63 external d 64 MatFDColoring e 65 PetscErrorCode z 66 end subroutine 67 subroutine SNESSetJacobian2(a,b,c,d,e,z) 68 use petscsnes 69 SNES a 70 Mat b 71 Mat c 72 external d 73 integer e 74 PetscErrorCode z 75 end subroutine 76 end interface 77! 78! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79! Variable declarations 80! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 81! 82! Variables: 83! snes - nonlinear solver 84! x, r - solution, residual vectors 85! J - Jacobian matrix 86! its - iterations for convergence 87! matrix_free - flag - 1 indicates matrix-free version 88! lambda - nonlinearity parameter 89! draw - drawing context 90! 91 SNES snes 92 MatColoring mc 93 Vec x,r 94 PetscDraw draw 95 Mat J 96 PetscBool matrix_free,flg,fd_coloring 97 PetscErrorCode ierr 98 PetscInt its,N, mx,my,i5 99 PetscMPIInt size,rank 100 PetscReal lambda_max,lambda_min,lambda 101 MatFDColoring fdcoloring 102 ISColoring iscoloring 103 PetscBool pc 104 external postcheck 105 106 PetscScalar lx_v(0:1) 107 PetscOffset lx_i 108 109! Store parameters in common block 110 111 common /params/ lambda,mx,my,fd_coloring 112 113! Note: Any user-defined Fortran routines (such as FormJacobian) 114! MUST be declared as external. 115 116 external FormFunction,FormInitialGuess,FormJacobian 117 118! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119! Initialize program 120! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 121 122 PetscCallA(PetscInitialize(ierr)) 123 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)) 124 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)) 125 126 if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif 127 128! Initialize problem parameters 129 i5 = 5 130 lambda_max = 6.81 131 lambda_min = 0.0 132 lambda = 6.0 133 mx = 4 134 my = 4 135 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)) 136 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)) 137 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)) 138 if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif 139 N = mx*my 140 pc = PETSC_FALSE 141 PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr)) 142 143! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 144! Create nonlinear solver context 145! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 147 PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr)) 148 149 if (pc .eqv. PETSC_TRUE) then 150 PetscCallA(SNESSetType(snes,SNESNEWTONTR,ierr)) 151 PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr)) 152 endif 153 154! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155! Create vector data structures; set function evaluation routine 156! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 157 158 PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr)) 159 PetscCallA(VecSetSizes(x,PETSC_DECIDE,N,ierr)) 160 PetscCallA(VecSetFromOptions(x,ierr)) 161 PetscCallA(VecDuplicate(x,r,ierr)) 162 163! Set function evaluation routine and vector. Whenever the nonlinear 164! solver needs to evaluate the nonlinear function, it will call this 165! routine. 166! - Note that the final routine argument is the user-defined 167! context that provides application-specific data for the 168! function evaluation routine. 169 170 PetscCallA(SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr)) 171 172! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173! Create matrix data structure; set Jacobian evaluation routine 174! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 175 176! Create matrix. Here we only approximately preallocate storage space 177! for the Jacobian. See the users manual for a discussion of better 178! techniques for preallocating matrix memory. 179 180 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr)) 181 if (.not. matrix_free) then 182 PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr)) 183 endif 184 185! 186! This option will cause the Jacobian to be computed via finite differences 187! efficiently using a coloring of the columns of the matrix. 188! 189 fd_coloring = .false. 190 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr)) 191 if (fd_coloring) then 192 193! 194! This initializes the nonzero structure of the Jacobian. This is artificial 195! because clearly if we had a routine to compute the Jacobian we won't need 196! to use finite differences. 197! 198 PetscCallA(FormJacobian(snes,x,J,J,0,ierr)) 199! 200! Color the matrix, i.e. determine groups of columns that share no common 201! rows. These columns in the Jacobian can all be computed simultaneously. 202! 203 PetscCallA(MatColoringCreate(J,mc,ierr)) 204 PetscCallA(MatColoringSetType(mc,MATCOLORINGNATURAL,ierr)) 205 PetscCallA(MatColoringSetFromOptions(mc,ierr)) 206 PetscCallA(MatColoringApply(mc,iscoloring,ierr)) 207 PetscCallA(MatColoringDestroy(mc,ierr)) 208! 209! Create the data structure that SNESComputeJacobianDefaultColor() uses 210! to compute the actual Jacobians via finite differences. 211! 212 PetscCallA(MatFDColoringCreate(J,iscoloring,fdcoloring,ierr)) 213 PetscCallA(MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr)) 214 PetscCallA(MatFDColoringSetFromOptions(fdcoloring,ierr)) 215 PetscCallA(MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr)) 216! 217! Tell SNES to use the routine SNESComputeJacobianDefaultColor() 218! to compute Jacobians. 219! 220 PetscCallA(SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr)) 221 PetscCallA(ISColoringDestroy(iscoloring,ierr)) 222 223 else if (.not. matrix_free) then 224 225! Set Jacobian matrix data structure and default Jacobian evaluation 226! routine. Whenever the nonlinear solver needs to compute the 227! Jacobian matrix, it will call this routine. 228! - Note that the final routine argument is the user-defined 229! context that provides application-specific data for the 230! Jacobian evaluation routine. 231! - The user can override with: 232! -snes_fd : default finite differencing approximation of Jacobian 233! -snes_mf : matrix-free Newton-Krylov method with no preconditioning 234! (unless user explicitly sets preconditioner) 235! -snes_mf_operator : form preconditioning matrix as set by the user, 236! but use matrix-free approx for Jacobian-vector 237! products within Newton-Krylov method 238! 239 PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)) 240 endif 241 242! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243! Customize nonlinear solver; set runtime options 244! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 245 246! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 247 248 PetscCallA(SNESSetFromOptions(snes,ierr)) 249 250! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 251! Evaluate initial guess; then solve nonlinear system. 252! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 253 254! Note: The user should initialize the vector, x, with the initial guess 255! for the nonlinear solver prior to calling SNESSolve(). In particular, 256! to employ an initial guess of zero, the user should explicitly set 257! this vector to zero by calling VecSet(). 258 259 PetscCallA(FormInitialGuess(x,ierr)) 260 PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr)) 261 PetscCallA(SNESGetIterationNumber(snes,its,ierr)) 262 if (rank .eq. 0) then 263 write(6,100) its 264 endif 265 100 format('Number of SNES iterations = ',i1) 266 267! PetscDraw contour plot of solution 268 269 PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr)) 270 PetscCallA(PetscDrawSetFromOptions(draw,ierr)) 271 272 PetscCallA(VecGetArrayRead(x,lx_v,lx_i,ierr)) 273 PetscCallA(PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr)) 274 PetscCallA(VecRestoreArrayRead(x,lx_v,lx_i,ierr)) 275 276! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 277! Free work space. All PETSc objects should be destroyed when they 278! are no longer needed. 279! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 280 281 if (.not. matrix_free) PetscCallA(MatDestroy(J,ierr)) 282 if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring,ierr)) 283 284 PetscCallA(VecDestroy(x,ierr)) 285 PetscCallA(VecDestroy(r,ierr)) 286 PetscCallA(SNESDestroy(snes,ierr)) 287 PetscCallA(PetscDrawDestroy(draw,ierr)) 288 PetscCallA(PetscFinalize(ierr)) 289 end 290 291! --------------------------------------------------------------------- 292! 293! FormInitialGuess - Forms initial approximation. 294! 295! Input Parameter: 296! X - vector 297! 298! Output Parameters: 299! X - vector 300! ierr - error code 301! 302! Notes: 303! This routine serves as a wrapper for the lower-level routine 304! "ApplicationInitialGuess", where the actual computations are 305! done using the standard Fortran style of treating the local 306! vector data as a multidimensional array over the local mesh. 307! This routine merely accesses the local vector data via 308! VecGetArray() and VecRestoreArray(). 309! 310 subroutine FormInitialGuess(X,ierr) 311 use petscsnes 312 implicit none 313 314! Input/output variables: 315 Vec X 316 PetscErrorCode ierr 317 318! Declarations for use with local arrays: 319 PetscScalar lx_v(0:1) 320 PetscOffset lx_i 321 322 ierr = 0 323 324! Get a pointer to vector data. 325! - For default PETSc vectors, VecGetArray() returns a pointer to 326! the data array. Otherwise, the routine is implementation dependent. 327! - You MUST call VecRestoreArray() when you no longer need access to 328! the array. 329! - Note that the Fortran interface to VecGetArray() differs from the 330! C version. See the users manual for details. 331 332 PetscCallA(VecGetArray(X,lx_v,lx_i,ierr)) 333 334! Compute initial guess 335 336 PetscCallA(ApplicationInitialGuess(lx_v(lx_i),ierr)) 337 338! Restore vector 339 340 PetscCallA(VecRestoreArray(X,lx_v,lx_i,ierr)) 341 342 return 343 end 344 345! --------------------------------------------------------------------- 346! 347! ApplicationInitialGuess - Computes initial approximation, called by 348! the higher level routine FormInitialGuess(). 349! 350! Input Parameter: 351! x - local vector data 352! 353! Output Parameters: 354! f - local vector data, f(x) 355! ierr - error code 356! 357! Notes: 358! This routine uses standard Fortran-style computations over a 2-dim array. 359! 360 subroutine ApplicationInitialGuess(x,ierr) 361 use petscksp 362 implicit none 363 364! Common blocks: 365 PetscReal lambda 366 PetscInt mx,my 367 PetscBool fd_coloring 368 common /params/ lambda,mx,my,fd_coloring 369 370! Input/output variables: 371 PetscScalar x(mx,my) 372 PetscErrorCode ierr 373 374! Local variables: 375 PetscInt i,j 376 PetscReal temp1,temp,hx,hy,one 377 378! Set parameters 379 380 ierr = 0 381 one = 1.0 382 hx = one/(mx-1) 383 hy = one/(my-1) 384 temp1 = lambda/(lambda + one) 385 386 do 20 j=1,my 387 temp = min(j-1,my-j)*hy 388 do 10 i=1,mx 389 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 390 x(i,j) = 0.0 391 else 392 x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp)) 393 endif 394 10 continue 395 20 continue 396 397 return 398 end 399 400! --------------------------------------------------------------------- 401! 402! FormFunction - Evaluates nonlinear function, F(x). 403! 404! Input Parameters: 405! snes - the SNES context 406! X - input vector 407! dummy - optional user-defined context, as set by SNESSetFunction() 408! (not used here) 409! 410! Output Parameter: 411! F - vector with newly computed function 412! 413! Notes: 414! This routine serves as a wrapper for the lower-level routine 415! "ApplicationFunction", where the actual computations are 416! done using the standard Fortran style of treating the local 417! vector data as a multidimensional array over the local mesh. 418! This routine merely accesses the local vector data via 419! VecGetArray() and VecRestoreArray(). 420! 421 subroutine FormFunction(snes,X,F,fdcoloring,ierr) 422 use petscsnes 423 implicit none 424 425! Input/output variables: 426 SNES snes 427 Vec X,F 428 PetscErrorCode ierr 429 MatFDColoring fdcoloring 430 431! Common blocks: 432 PetscReal lambda 433 PetscInt mx,my 434 PetscBool fd_coloring 435 common /params/ lambda,mx,my,fd_coloring 436 437! Declarations for use with local arrays: 438 PetscScalar lx_v(0:1),lf_v(0:1) 439 PetscOffset lx_i,lf_i 440 441 PetscInt, pointer :: indices(:) 442 443! Get pointers to vector data. 444! - For default PETSc vectors, VecGetArray() returns a pointer to 445! the data array. Otherwise, the routine is implementation dependent. 446! - You MUST call VecRestoreArray() when you no longer need access to 447! the array. 448! - Note that the Fortran interface to VecGetArray() differs from the 449! C version. See the Fortran chapter of the users manual for details. 450 451 PetscCallA(VecGetArrayRead(X,lx_v,lx_i,ierr)) 452 PetscCallA(VecGetArray(F,lf_v,lf_i,ierr)) 453 454! Compute function 455 456 PetscCallA(ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)) 457 458! Restore vectors 459 460 PetscCallA(VecRestoreArrayRead(X,lx_v,lx_i,ierr)) 461 PetscCallA(VecRestoreArray(F,lf_v,lf_i,ierr)) 462 463 PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr)) 464! 465! fdcoloring is in the common block and used here ONLY to test the 466! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90() 467! 468 if (fd_coloring) then 469 PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)) 470 print*,'Indices from GetPerturbedColumnsF90' 471 write(*,1000) indices 472 1000 format(50i4) 473 PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)) 474 endif 475 return 476 end 477 478! --------------------------------------------------------------------- 479! 480! ApplicationFunction - Computes nonlinear function, called by 481! the higher level routine FormFunction(). 482! 483! Input Parameter: 484! x - local vector data 485! 486! Output Parameters: 487! f - local vector data, f(x) 488! ierr - error code 489! 490! Notes: 491! This routine uses standard Fortran-style computations over a 2-dim array. 492! 493 subroutine ApplicationFunction(x,f,ierr) 494 use petscsnes 495 implicit none 496 497! Common blocks: 498 PetscReal lambda 499 PetscInt mx,my 500 PetscBool fd_coloring 501 common /params/ lambda,mx,my,fd_coloring 502 503! Input/output variables: 504 PetscScalar x(mx,my),f(mx,my) 505 PetscErrorCode ierr 506 507! Local variables: 508 PetscScalar two,one,hx,hy 509 PetscScalar hxdhy,hydhx,sc 510 PetscScalar u,uxx,uyy 511 PetscInt i,j 512 513 ierr = 0 514 one = 1.0 515 two = 2.0 516 hx = one/(mx-1) 517 hy = one/(my-1) 518 sc = hx*hy*lambda 519 hxdhy = hx/hy 520 hydhx = hy/hx 521 522! Compute function 523 524 do 20 j=1,my 525 do 10 i=1,mx 526 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 527 f(i,j) = x(i,j) 528 else 529 u = x(i,j) 530 uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j)) 531 uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1)) 532 f(i,j) = uxx + uyy - sc*exp(u) 533 endif 534 10 continue 535 20 continue 536 537 return 538 end 539 540! --------------------------------------------------------------------- 541! 542! FormJacobian - Evaluates Jacobian matrix. 543! 544! Input Parameters: 545! snes - the SNES context 546! x - input vector 547! dummy - optional user-defined context, as set by SNESSetJacobian() 548! (not used here) 549! 550! Output Parameters: 551! jac - Jacobian matrix 552! jac_prec - optionally different preconditioning matrix (not used here) 553! flag - flag indicating matrix structure 554! 555! Notes: 556! This routine serves as a wrapper for the lower-level routine 557! "ApplicationJacobian", where the actual computations are 558! done using the standard Fortran style of treating the local 559! vector data as a multidimensional array over the local mesh. 560! This routine merely accesses the local vector data via 561! VecGetArray() and VecRestoreArray(). 562! 563 subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr) 564 use petscsnes 565 implicit none 566 567! Input/output variables: 568 SNES snes 569 Vec X 570 Mat jac,jac_prec 571 PetscErrorCode ierr 572 integer dummy 573 574! Common blocks: 575 PetscReal lambda 576 PetscInt mx,my 577 PetscBool fd_coloring 578 common /params/ lambda,mx,my,fd_coloring 579 580! Declarations for use with local array: 581 PetscScalar lx_v(0:1) 582 PetscOffset lx_i 583 584! Get a pointer to vector data 585 586 PetscCallA(VecGetArrayRead(X,lx_v,lx_i,ierr)) 587 588! Compute Jacobian entries 589 590 PetscCallA(ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)) 591 592! Restore vector 593 594 PetscCallA(VecRestoreArrayRead(X,lx_v,lx_i,ierr)) 595 596! Assemble matrix 597 598 PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 599 PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 600 601 return 602 end 603 604! --------------------------------------------------------------------- 605! 606! ApplicationJacobian - Computes Jacobian matrix, called by 607! the higher level routine FormJacobian(). 608! 609! Input Parameters: 610! x - local vector data 611! 612! Output Parameters: 613! jac - Jacobian matrix 614! jac_prec - optionally different preconditioning matrix (not used here) 615! ierr - error code 616! 617! Notes: 618! This routine uses standard Fortran-style computations over a 2-dim array. 619! 620 subroutine ApplicationJacobian(x,jac,jac_prec,ierr) 621 use petscsnes 622 implicit none 623 624! Common blocks: 625 PetscReal lambda 626 PetscInt mx,my 627 PetscBool fd_coloring 628 common /params/ lambda,mx,my,fd_coloring 629 630! Input/output variables: 631 PetscScalar x(mx,my) 632 Mat jac,jac_prec 633 PetscErrorCode ierr 634 635! Local variables: 636 PetscInt i,j,row(1),col(5),i1,i5 637 PetscScalar two,one, hx,hy 638 PetscScalar hxdhy,hydhx,sc,v(5) 639 640! Set parameters 641 642 i1 = 1 643 i5 = 5 644 one = 1.0 645 two = 2.0 646 hx = one/(mx-1) 647 hy = one/(my-1) 648 sc = hx*hy 649 hxdhy = hx/hy 650 hydhx = hy/hx 651 652! Compute entries of the Jacobian matrix 653! - Here, we set all entries for a particular row at once. 654! - Note that MatSetValues() uses 0-based row and column numbers 655! in Fortran as well as in C. 656 657 do 20 j=1,my 658 row(1) = (j-1)*mx - 1 659 do 10 i=1,mx 660 row(1) = row(1) + 1 661! boundary points 662 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 663 PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)) 664! interior grid points 665 else 666 v(1) = -hxdhy 667 v(2) = -hydhx 668 v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j)) 669 v(4) = -hydhx 670 v(5) = -hxdhy 671 col(1) = row(1) - mx 672 col(2) = row(1) - 1 673 col(3) = row(1) 674 col(4) = row(1) + 1 675 col(5) = row(1) + mx 676 PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)) 677 endif 678 10 continue 679 20 continue 680 681 return 682 end 683 684! 685!/*TEST 686! 687! build: 688! requires: !single 689! 690! test: 691! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always 692! 693! test: 694! suffix: 2 695! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always 696! 697! test: 698! suffix: 3 699! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always 700! filter: sort -b 701! filter_output: sort -b 702! 703! test: 704! suffix: 4 705! args: -pc -par 6.807 -nox 706! 707!TEST*/ 708