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! ApplicationInitialGuess - Computes initial approximation, called by 346! the higher level routine FormInitialGuess(). 347! 348! Input Parameter: 349! x - local vector data 350! 351! Output Parameters: 352! f - local vector data, f(x) 353! ierr - error code 354! 355! Notes: 356! This routine uses standard Fortran-style computations over a 2-dim array. 357! 358 subroutine ApplicationInitialGuess(x,ierr) 359 use petscksp 360 implicit none 361 362! Common blocks: 363 PetscReal lambda 364 PetscInt mx,my 365 PetscBool fd_coloring 366 common /params/ lambda,mx,my,fd_coloring 367 368! Input/output variables: 369 PetscScalar x(mx,my) 370 PetscErrorCode ierr 371 372! Local variables: 373 PetscInt i,j 374 PetscReal temp1,temp,hx,hy,one 375 376! Set parameters 377 378 ierr = 0 379 one = 1.0 380 hx = one/(mx-1) 381 hy = one/(my-1) 382 temp1 = lambda/(lambda + one) 383 384 do 20 j=1,my 385 temp = min(j-1,my-j)*hy 386 do 10 i=1,mx 387 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 388 x(i,j) = 0.0 389 else 390 x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp)) 391 endif 392 10 continue 393 20 continue 394 395 return 396 end 397 398! --------------------------------------------------------------------- 399! 400! FormFunction - Evaluates nonlinear function, F(x). 401! 402! Input Parameters: 403! snes - the SNES context 404! X - input vector 405! dummy - optional user-defined context, as set by SNESSetFunction() 406! (not used here) 407! 408! Output Parameter: 409! F - vector with newly computed function 410! 411! Notes: 412! This routine serves as a wrapper for the lower-level routine 413! "ApplicationFunction", where the actual computations are 414! done using the standard Fortran style of treating the local 415! vector data as a multidimensional array over the local mesh. 416! This routine merely accesses the local vector data via 417! VecGetArray() and VecRestoreArray(). 418! 419 subroutine FormFunction(snes,X,F,fdcoloring,ierr) 420 use petscsnes 421 implicit none 422 423! Input/output variables: 424 SNES snes 425 Vec X,F 426 PetscErrorCode ierr 427 MatFDColoring fdcoloring 428 429! Common blocks: 430 PetscReal lambda 431 PetscInt mx,my 432 PetscBool fd_coloring 433 common /params/ lambda,mx,my,fd_coloring 434 435! Declarations for use with local arrays: 436 PetscScalar lx_v(0:1),lf_v(0:1) 437 PetscOffset lx_i,lf_i 438 439 PetscInt, pointer :: indices(:) 440 441! Get pointers to vector data. 442! - For default PETSc vectors, VecGetArray() returns a pointer to 443! the data array. Otherwise, the routine is implementation dependent. 444! - You MUST call VecRestoreArray() when you no longer need access to 445! the array. 446! - Note that the Fortran interface to VecGetArray() differs from the 447! C version. See the Fortran chapter of the users manual for details. 448 449 PetscCallA(VecGetArrayRead(X,lx_v,lx_i,ierr)) 450 PetscCallA(VecGetArray(F,lf_v,lf_i,ierr)) 451 452! Compute function 453 454 PetscCallA(ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr)) 455 456! Restore vectors 457 458 PetscCallA(VecRestoreArrayRead(X,lx_v,lx_i,ierr)) 459 PetscCallA(VecRestoreArray(F,lf_v,lf_i,ierr)) 460 461 PetscCallA(PetscLogFlops(11.0d0*mx*my,ierr)) 462! 463! fdcoloring is in the common block and used here ONLY to test the 464! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90() 465! 466 if (fd_coloring) then 467 PetscCallA(MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr)) 468 print*,'Indices from GetPerturbedColumnsF90' 469 write(*,1000) indices 470 1000 format(50i4) 471 PetscCallA(MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr)) 472 endif 473 return 474 end 475 476! --------------------------------------------------------------------- 477! 478! ApplicationFunction - Computes nonlinear function, called by 479! the higher level routine FormFunction(). 480! 481! Input Parameter: 482! x - local vector data 483! 484! Output Parameters: 485! f - local vector data, f(x) 486! ierr - error code 487! 488! Notes: 489! This routine uses standard Fortran-style computations over a 2-dim array. 490! 491 subroutine ApplicationFunction(x,f,ierr) 492 use petscsnes 493 implicit none 494 495! Common blocks: 496 PetscReal lambda 497 PetscInt mx,my 498 PetscBool fd_coloring 499 common /params/ lambda,mx,my,fd_coloring 500 501! Input/output variables: 502 PetscScalar x(mx,my),f(mx,my) 503 PetscErrorCode ierr 504 505! Local variables: 506 PetscScalar two,one,hx,hy 507 PetscScalar hxdhy,hydhx,sc 508 PetscScalar u,uxx,uyy 509 PetscInt i,j 510 511 ierr = 0 512 one = 1.0 513 two = 2.0 514 hx = one/(mx-1) 515 hy = one/(my-1) 516 sc = hx*hy*lambda 517 hxdhy = hx/hy 518 hydhx = hy/hx 519 520! Compute function 521 522 do 20 j=1,my 523 do 10 i=1,mx 524 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 525 f(i,j) = x(i,j) 526 else 527 u = x(i,j) 528 uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j)) 529 uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1)) 530 f(i,j) = uxx + uyy - sc*exp(u) 531 endif 532 10 continue 533 20 continue 534 535 return 536 end 537 538! --------------------------------------------------------------------- 539! 540! FormJacobian - Evaluates Jacobian matrix. 541! 542! Input Parameters: 543! snes - the SNES context 544! x - input vector 545! dummy - optional user-defined context, as set by SNESSetJacobian() 546! (not used here) 547! 548! Output Parameters: 549! jac - Jacobian matrix 550! jac_prec - optionally different preconditioning matrix (not used here) 551! flag - flag indicating matrix structure 552! 553! Notes: 554! This routine serves as a wrapper for the lower-level routine 555! "ApplicationJacobian", where the actual computations are 556! done using the standard Fortran style of treating the local 557! vector data as a multidimensional array over the local mesh. 558! This routine merely accesses the local vector data via 559! VecGetArray() and VecRestoreArray(). 560! 561 subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr) 562 use petscsnes 563 implicit none 564 565! Input/output variables: 566 SNES snes 567 Vec X 568 Mat jac,jac_prec 569 PetscErrorCode ierr 570 integer dummy 571 572! Common blocks: 573 PetscReal lambda 574 PetscInt mx,my 575 PetscBool fd_coloring 576 common /params/ lambda,mx,my,fd_coloring 577 578! Declarations for use with local array: 579 PetscScalar lx_v(0:1) 580 PetscOffset lx_i 581 582! Get a pointer to vector data 583 584 PetscCallA(VecGetArrayRead(X,lx_v,lx_i,ierr)) 585 586! Compute Jacobian entries 587 588 PetscCallA(ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr)) 589 590! Restore vector 591 592 PetscCallA(VecRestoreArrayRead(X,lx_v,lx_i,ierr)) 593 594! Assemble matrix 595 596 PetscCallA(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 597 PetscCallA(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 598 599 return 600 end 601 602! --------------------------------------------------------------------- 603! 604! ApplicationJacobian - Computes Jacobian matrix, called by 605! the higher level routine FormJacobian(). 606! 607! Input Parameters: 608! x - local vector data 609! 610! Output Parameters: 611! jac - Jacobian matrix 612! jac_prec - optionally different preconditioning matrix (not used here) 613! ierr - error code 614! 615! Notes: 616! This routine uses standard Fortran-style computations over a 2-dim array. 617! 618 subroutine ApplicationJacobian(x,jac,jac_prec,ierr) 619 use petscsnes 620 implicit none 621 622! Common blocks: 623 PetscReal lambda 624 PetscInt mx,my 625 PetscBool fd_coloring 626 common /params/ lambda,mx,my,fd_coloring 627 628! Input/output variables: 629 PetscScalar x(mx,my) 630 Mat jac,jac_prec 631 PetscErrorCode ierr 632 633! Local variables: 634 PetscInt i,j,row(1),col(5),i1,i5 635 PetscScalar two,one, hx,hy 636 PetscScalar hxdhy,hydhx,sc,v(5) 637 638! Set parameters 639 640 i1 = 1 641 i5 = 5 642 one = 1.0 643 two = 2.0 644 hx = one/(mx-1) 645 hy = one/(my-1) 646 sc = hx*hy 647 hxdhy = hx/hy 648 hydhx = hy/hx 649 650! Compute entries of the Jacobian matrix 651! - Here, we set all entries for a particular row at once. 652! - Note that MatSetValues() uses 0-based row and column numbers 653! in Fortran as well as in C. 654 655 do 20 j=1,my 656 row(1) = (j-1)*mx - 1 657 do 10 i=1,mx 658 row(1) = row(1) + 1 659! boundary points 660 if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 661 PetscCallA(MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr)) 662! interior grid points 663 else 664 v(1) = -hxdhy 665 v(2) = -hydhx 666 v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j)) 667 v(4) = -hydhx 668 v(5) = -hxdhy 669 col(1) = row(1) - mx 670 col(2) = row(1) - 1 671 col(3) = row(1) 672 col(4) = row(1) + 1 673 col(5) = row(1) + mx 674 PetscCallA(MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr)) 675 endif 676 10 continue 677 20 continue 678 679 return 680 end 681 682! 683!/*TEST 684! 685! build: 686! requires: !single 687! 688! test: 689! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always 690! 691! test: 692! suffix: 2 693! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always 694! 695! test: 696! suffix: 3 697! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always 698! filter: sort -b 699! filter_output: sort -b 700! 701! test: 702! suffix: 4 703! args: -pc -par 6.807 -nox 704! 705!TEST*/ 706