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