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! This system (A) is augmented with constraints: 10! 11! A -B * phi = rho 12! -C I lam = 0 13! 14! where I is the identity, A is the "normal" Poisson equation, B is the "distributor" of the 15! total flux (the first block equation is the flux surface averaging equation). The second 16! equation lambda = C * x enforces the surface flux auxiliary equation. B and C have all 17! positive entries, areas in C and fraction of area in B. 18! 19 20! 21! -------------------------------------------------------------------------- 22! 23! Solid Fuel Ignition (SFI) problem. This problem is modeled by 24! the partial differential equation 25! 26! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 27! 28! with boundary conditions 29! 30! u = 0 for x = 0, x = 1, y = 0, y = 1. 31! 32! A finite difference approximation with the usual 5-point stencil 33! is used to discretize the boundary value problem to obtain a nonlinear 34! system of equations. 35! 36! -------------------------------------------------------------------------- 37! The following define must be used before including any PETSc include files 38! into a module or interface. This is because they can't handle declarations 39! in them 40! 41 module ex73f90tmodule 42#include <petsc/finclude/petscdm.h> 43#include <petsc/finclude/petscmat.h> 44 use petscdm 45 use petscmat 46 type ex73f90tmodule_type 47 DM::da 48! temp A block stuff 49 PetscInt mx,my 50 PetscMPIInt rank 51 PetscReal lambda 52! Mats 53 Mat::Amat,AmatLin,Bmat,CMat,Dmat 54 IS::isPhi,isLambda 55 end type ex73f90tmodule_type 56 57 end module ex73f90tmodule 58 59 module ex73f90tmodule_interfaces 60 use ex73f90tmodule 61 62 Interface SNESSetApplicationContext 63 Subroutine SNESSetApplicationContext(snesIn,ctx,ierr) 64#include <petsc/finclude/petscsnes.h> 65 use petscsnes 66 use ex73f90tmodule 67 SNES:: snesIn 68 type(ex73f90tmodule_type) ctx 69 PetscErrorCode ierr 70 End Subroutine 71 End Interface SNESSetApplicationContext 72 73 Interface SNESGetApplicationContext 74 Subroutine SNESGetApplicationContext(snesIn,ctx,ierr) 75#include <petsc/finclude/petscsnes.h> 76 use petscsnes 77 use ex73f90tmodule 78 SNES:: snesIn 79 type(ex73f90tmodule_type), pointer :: ctx 80 PetscErrorCode ierr 81 End Subroutine 82 End Interface SNESGetApplicationContext 83 end module ex73f90tmodule_interfaces 84 85 program main 86#include <petsc/finclude/petscdm.h> 87#include <petsc/finclude/petscsnes.h> 88 use petscdm 89 use petscdmda 90 use petscsnes 91 use ex73f90tmodule 92 use ex73f90tmodule_interfaces 93 implicit none 94! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95! Variable declarations 96! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97! 98! Variables: 99! mysnes - nonlinear solver 100! x, r - solution, residual vectors 101! J - Jacobian matrix 102! its - iterations for convergence 103! Nx, Ny - number of preocessors in x- and y- directions 104! 105 SNES:: mysnes 106 Vec:: x,r,x2,x1,x1loc,x2loc 107 Mat:: Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4) 108! Mat:: tmat 109 DM:: daphi,dalam 110 IS:: isglobal(2) 111 PetscErrorCode ierr 112 PetscInt its,N1,N2,i,j,irow,row(1) 113 PetscInt col(1),low,high,lamlow,lamhigh 114 PetscBool flg 115 PetscInt ione,nfour,itwo,nloc,nloclam 116 PetscReal lambda_max,lambda_min 117 type(ex73f90tmodule_type) solver 118 PetscScalar bval(1),cval(1),one 119 120! Note: Any user-defined Fortran routines (such as FormJacobian) 121! MUST be declared as external. 122 external FormInitialGuess,FormJacobian,FormFunction 123 124! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 125! Initialize program 126! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 PetscCallA(PetscInitialize(ierr)) 128 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr)) 129 130! Initialize problem parameters 131 lambda_max = 6.81_PETSC_REAL_KIND 132 lambda_min = 0.0 133 solver%lambda = 6.0 134 ione = 1 135 nfour = 4 136 itwo = 2 137 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr)) 138 if (solver%lambda .ge. lambda_max .or. solver%lambda .lt. lambda_min) then 139 SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range') 140 endif 141 142! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143! Create vector data structures; set function evaluation routine 144! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 146! just get size 147 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,daphi,ierr)) 148 PetscCallA(DMSetFromOptions(daphi,ierr)) 149 PetscCallA(DMSetUp(daphi,ierr)) 150 PetscCallA(DMDAGetInfo(daphi,PETSC_NULL_INTEGER,solver%mx,solver%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)) 151 N1 = solver%my*solver%mx 152 N2 = solver%my 153 flg = .false. 154 PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr)) 155 if (flg) then 156 N2 = 0 157 endif 158 159 PetscCallA(DMDestroy(daphi,ierr)) 160 161! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 162! Create matrix data structure; set Jacobian evaluation routine 163! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 164 PetscCallA(DMShellCreate(PETSC_COMM_WORLD,daphi,ierr)) 165 PetscCallA(DMSetOptionsPrefix(daphi,'phi_',ierr)) 166 PetscCallA(DMSetFromOptions(daphi,ierr)) 167 168 PetscCallA(VecCreate(PETSC_COMM_WORLD,x1,ierr)) 169 PetscCallA(VecSetSizes(x1,PETSC_DECIDE,N1,ierr)) 170 PetscCallA(VecSetFromOptions(x1,ierr)) 171 172 PetscCallA(VecGetOwnershipRange(x1,low,high,ierr)) 173 nloc = high - low 174 175 PetscCallA(MatCreate(PETSC_COMM_WORLD,Amat,ierr)) 176 PetscCallA(MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr)) 177 PetscCallA(MatSetUp(Amat,ierr)) 178 179 PetscCallA(MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr)) 180 PetscCallA(MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr)) 181 PetscCallA(MatSetUp(solver%AmatLin,ierr)) 182 183 PetscCallA(FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr)) 184 PetscCallA(MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr)) 185 PetscCallA(MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr)) 186 187 PetscCallA(DMShellSetGlobalVector(daphi,x1,ierr)) 188 PetscCallA(DMShellSetMatrix(daphi,Amat,ierr)) 189 190 PetscCallA(VecCreate(PETSC_COMM_SELF,x1loc,ierr)) 191 PetscCallA(VecSetSizes(x1loc,nloc,nloc,ierr)) 192 PetscCallA(VecSetFromOptions(x1loc,ierr)) 193 PetscCallA(DMShellSetLocalVector(daphi,x1loc,ierr)) 194 195! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196! Create B, C, & D matrices 197! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 198 PetscCallA(MatCreate(PETSC_COMM_WORLD,Cmat,ierr)) 199 PetscCallA(MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr)) 200 PetscCallA(MatSetUp(Cmat,ierr)) 201! create data for C and B 202 PetscCallA(MatCreate(PETSC_COMM_WORLD,Bmat,ierr)) 203 PetscCallA(MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr)) 204 PetscCallA(MatSetUp(Bmat,ierr)) 205! create data for D 206 PetscCallA(MatCreate(PETSC_COMM_WORLD,Dmat,ierr)) 207 PetscCallA(MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr)) 208 PetscCallA(MatSetUp(Dmat,ierr)) 209 210 PetscCallA(VecCreate(PETSC_COMM_WORLD,x2,ierr)) 211 PetscCallA(VecSetSizes(x2,PETSC_DECIDE,N2,ierr)) 212 PetscCallA(VecSetFromOptions(x2,ierr)) 213 214 PetscCallA(VecGetOwnershipRange(x2,lamlow,lamhigh,ierr)) 215 nloclam = lamhigh-lamlow 216 217! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 218! Set fake B and C 219! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 220 one = 1.0 221 if (N2 .gt. 0) then 222 bval(1) = -one/(solver%mx-2) 223! cval = -one/(solver%my*solver%mx) 224 cval(1) = -one 225 do 20 irow=low,high-1 226 j = irow/solver%mx ! row in domain 227 i = mod(irow,solver%mx) 228 row(1) = irow 229 col(1) = j 230 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 231 ! no op 232 else 233 PetscCallA(MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr)) 234 endif 235 row(1) = j 236 PetscCallA(MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr)) 237 20 continue 238 endif 239 PetscCallA(MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr)) 240 PetscCallA(MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr)) 241 PetscCallA(MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr)) 242 PetscCallA(MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr)) 243 244! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 245! Set D (indentity) 246! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 247 do 30 j=lamlow,lamhigh-1 248 row(1) = j 249 cval(1) = one 250 PetscCallA(MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr)) 251 30 continue 252 PetscCallA(MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr)) 253 PetscCallA(MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr)) 254 255! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 256! DM for lambda (dalam) : temp driver for A block, setup A block solver data 257! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 258 PetscCallA(DMShellCreate(PETSC_COMM_WORLD,dalam,ierr)) 259 PetscCallA(DMShellSetGlobalVector(dalam,x2,ierr)) 260 PetscCallA(DMShellSetMatrix(dalam,Dmat,ierr)) 261 262 PetscCallA(VecCreate(PETSC_COMM_SELF,x2loc,ierr)) 263 PetscCallA(VecSetSizes(x2loc,nloclam,nloclam,ierr)) 264 PetscCallA(VecSetFromOptions(x2loc,ierr)) 265 PetscCallA(DMShellSetLocalVector(dalam,x2loc,ierr)) 266 267 PetscCallA(DMSetOptionsPrefix(dalam,'lambda_',ierr)) 268 PetscCallA(DMSetFromOptions(dalam,ierr)) 269! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 270! Create field split DA 271! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 272 PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr)) 273 PetscCallA(DMCompositeAddDM(solver%da,daphi,ierr)) 274 PetscCallA(DMCompositeAddDM(solver%da,dalam,ierr)) 275 PetscCallA(DMSetFromOptions(solver%da,ierr)) 276 PetscCallA(DMSetUp(solver%da,ierr)) 277 PetscCallA(DMCompositeGetGlobalISs(solver%da,isglobal,ierr)) 278 solver%isPhi = isglobal(1) 279 solver%isLambda = isglobal(2) 280 281! cache matrices 282 solver%Amat = Amat 283 solver%Bmat = Bmat 284 solver%Cmat = Cmat 285 solver%Dmat = Dmat 286 287 matArray(1) = Amat 288 matArray(2) = Bmat 289 matArray(3) = Cmat 290 matArray(4) = Dmat 291 292 PetscCallA(MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr)) 293 PetscCallA(MatSetFromOptions(KKTmat,ierr)) 294 295! Extract global and local vectors from DMDA; then duplicate for remaining 296! vectors that are the same types 297 PetscCallA(MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr)) 298 PetscCallA(VecDuplicate(x,r,ierr)) 299 300 PetscCallA(SNESCreate(PETSC_COMM_WORLD,mysnes,ierr)) 301 302 PetscCallA(SNESSetDM(mysnes,solver%da,ierr)) 303 304 PetscCallA(SNESSetApplicationContext(mysnes,solver,ierr)) 305 306 PetscCallA(SNESSetDM(mysnes,solver%da,ierr)) 307 308! Set function evaluation routine and vector 309 PetscCallA(SNESSetFunction(mysnes,r,FormFunction,solver,ierr)) 310 311 PetscCallA(SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr)) 312 313! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 314! Customize nonlinear solver; set runtime options 315! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 316! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 317 PetscCallA(SNESSetFromOptions(mysnes,ierr)) 318 319! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 320! Evaluate initial guess; then solve nonlinear system. 321! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 322! Note: The user should initialize the vector, x, with the initial guess 323! for the nonlinear solver prior to calling SNESSolve(). In particular, 324! to employ an initial guess of zero, the user should explicitly set 325! this vector to zero by calling VecSet(). 326 327 PetscCallA(FormInitialGuess(mysnes,x,ierr)) 328 PetscCallA(SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr)) 329 PetscCallA(SNESGetIterationNumber(mysnes,its,ierr)) 330 if (solver%rank .eq. 0) then 331 write(6,100) its 332 endif 333 100 format('Number of SNES iterations = ',i5) 334 335! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 336! Free work space. All PETSc objects should be destroyed when they 337! are no longer needed. 338! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 339 PetscCallA(MatDestroy(KKTmat,ierr)) 340 PetscCallA(MatDestroy(Amat,ierr)) 341 PetscCallA(MatDestroy(Dmat,ierr)) 342 PetscCallA(MatDestroy(Bmat,ierr)) 343 PetscCallA(MatDestroy(Cmat,ierr)) 344 PetscCallA(MatDestroy(solver%AmatLin,ierr)) 345 PetscCallA(ISDestroy(solver%isPhi,ierr)) 346 PetscCallA(ISDestroy(solver%isLambda,ierr)) 347 PetscCallA(VecDestroy(x,ierr)) 348 PetscCallA(VecDestroy(x2,ierr)) 349 PetscCallA(VecDestroy(x1,ierr)) 350 PetscCallA(VecDestroy(x1loc,ierr)) 351 PetscCallA(VecDestroy(x2loc,ierr)) 352 PetscCallA(VecDestroy(r,ierr)) 353 PetscCallA(SNESDestroy(mysnes,ierr)) 354 PetscCallA(DMDestroy(solver%da,ierr)) 355 PetscCallA(DMDestroy(daphi,ierr)) 356 PetscCallA(DMDestroy(dalam,ierr)) 357 358 PetscCallA(PetscFinalize(ierr)) 359 end 360 361! --------------------------------------------------------------------- 362! 363! FormInitialGuess - Forms initial approximation. 364! 365! Input Parameters: 366! X - vector 367! 368! Output Parameter: 369! X - vector 370! 371! Notes: 372! This routine serves as a wrapper for the lower-level routine 373! "InitialGuessLocal", where the actual computations are 374! done using the standard Fortran style of treating the local 375! vector data as a multidimensional array over the local mesh. 376! This routine merely handles ghost point scatters and accesses 377! the local vector data via VecGetArrayF90() and VecRestoreArrayF90(). 378! 379 subroutine FormInitialGuess(mysnes,Xnest,ierr) 380#include <petsc/finclude/petscsnes.h> 381 use petscsnes 382 use ex73f90tmodule 383 use ex73f90tmodule_interfaces 384 implicit none 385! Input/output variables: 386 SNES:: mysnes 387 Vec:: Xnest 388 PetscErrorCode ierr 389 390! Declarations for use with local arrays: 391 type(ex73f90tmodule_type), pointer:: solver 392 Vec:: Xsub(2) 393 PetscInt:: izero,ione,itwo 394 395 izero = 0 396 ione = 1 397 itwo = 2 398 ierr = 0 399 PetscCall(SNESGetApplicationContext(mysnes,solver,ierr)) 400 PetscCall(DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 401 402 PetscCall(InitialGuessLocal(solver,Xsub(1),ierr)) 403 PetscCall(VecAssemblyBegin(Xsub(1),ierr)) 404 PetscCall(VecAssemblyEnd(Xsub(1),ierr)) 405 406! zero out lambda 407 PetscCall(VecZeroEntries(Xsub(2),ierr)) 408 PetscCall(DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 409 410 return 411 end subroutine FormInitialGuess 412 413! --------------------------------------------------------------------- 414! 415! InitialGuessLocal - Computes initial approximation, called by 416! the higher level routine FormInitialGuess(). 417! 418! Input Parameter: 419! X1 - local vector data 420! 421! Output Parameters: 422! x - local vector data 423! ierr - error code 424! 425! Notes: 426! This routine uses standard Fortran-style computations over a 2-dim array. 427! 428 subroutine InitialGuessLocal(solver,X1,ierr) 429#include <petsc/finclude/petscsys.h> 430 use petscsys 431 use ex73f90tmodule 432 implicit none 433! Input/output variables: 434 type (ex73f90tmodule_type) solver 435 Vec:: X1 436 PetscErrorCode ierr 437 438! Local variables: 439 PetscInt row,i,j,ione,low,high 440 PetscReal temp1,temp,hx,hy,v 441 PetscReal one 442 443! Set parameters 444 ione = 1 445 ierr = 0 446 one = 1.0 447 hx = one/(solver%mx-1) 448 hy = one/(solver%my-1) 449 temp1 = solver%lambda/(solver%lambda + one) + one 450 451 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 452 453 do 20 row=low,high-1 454 j = row/solver%mx 455 i = mod(row,solver%mx) 456 temp = min(j,solver%my-j+1)*hy 457 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 458 v = 0.0 459 else 460 v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp)) 461 endif 462 PetscCall(VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr)) 463 20 continue 464 465 return 466 end subroutine InitialGuessLocal 467 468! --------------------------------------------------------------------- 469! 470! FormJacobian - Evaluates Jacobian matrix. 471! 472! Input Parameters: 473! dummy - the SNES context 474! x - input vector 475! solver - solver data 476! 477! Output Parameters: 478! jac - Jacobian matrix 479! jac_prec - optionally different preconditioning matrix (not used here) 480! flag - flag indicating matrix structure 481! 482 subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr) 483#include <petsc/finclude/petscsnes.h> 484 use petscsnes 485 use ex73f90tmodule 486 implicit none 487! Input/output variables: 488 SNES:: dummy 489 Vec:: X 490 Mat:: jac,jac_prec 491 type(ex73f90tmodule_type) solver 492 PetscErrorCode ierr 493 494! Declarations for use with local arrays: 495 Vec:: Xsub(1) 496 Mat:: Amat 497 PetscInt ione 498 499 ione = 1 500 501 PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 502 503! Compute entries for the locally owned part of the Jacobian preconditioner. 504 PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr)) 505 506 PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr)) 507 PetscCall(MatDestroy(Amat,ierr)) ! discard our reference 508 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 509 510 ! the rest of the matrix is not touched 511 PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 512 PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 513 if (jac .ne. jac_prec) then 514 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 515 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 516 end if 517 518! Tell the matrix we will never add a new nonzero location to the 519! matrix. If we do it will generate an error. 520 PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 521 522 return 523 end subroutine FormJacobian 524 525! --------------------------------------------------------------------- 526! 527! FormJacobianLocal - Computes Jacobian preconditioner matrix, 528! called by the higher level routine FormJacobian(). 529! 530! Input Parameters: 531! x - local vector data 532! 533! Output Parameters: 534! jac - Jacobian preconditioner matrix 535! ierr - error code 536! 537! Notes: 538! This routine uses standard Fortran-style computations over a 2-dim array. 539! 540 subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr) 541#include <petsc/finclude/petscmat.h> 542 use petscmat 543 use ex73f90tmodule 544 implicit none 545! Input/output variables: 546 type (ex73f90tmodule_type) solver 547 Vec:: X1 548 Mat:: jac 549 logical add_nl_term 550 PetscErrorCode ierr 551 552! Local variables: 553 PetscInt irow,row(1),col(5),i,j 554 PetscInt ione,ifive,low,high,ii 555 PetscScalar two,one,hx,hy,hy2inv 556 PetscScalar hx2inv,sc,v(5) 557 PetscScalar,pointer :: lx_v(:) 558 559! Set parameters 560 ione = 1 561 ifive = 5 562 one = 1.0 563 two = 2.0 564 hx = one/(solver%mx-1) 565 hy = one/(solver%my-1) 566 sc = solver%lambda 567 hx2inv = one/(hx*hx) 568 hy2inv = one/(hy*hy) 569 570 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 571 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 572 573 ii = 0 574 do 20 irow=low,high-1 575 j = irow/solver%mx 576 i = mod(irow,solver%mx) 577 ii = ii + 1 ! one based local index 578! boundary points 579 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 580 col(1) = irow 581 row(1) = irow 582 v(1) = one 583 PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr)) 584! interior grid points 585 else 586 v(1) = -hy2inv 587 if (j-1==0) v(1) = 0.0 588 v(2) = -hx2inv 589 if (i-1==0) v(2) = 0.0 590 v(3) = two*(hx2inv + hy2inv) 591 if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii)) 592 v(4) = -hx2inv 593 if (i+1==solver%mx-1) v(4) = 0.0 594 v(5) = -hy2inv 595 if (j+1==solver%my-1) v(5) = 0.0 596 col(1) = irow - solver%mx 597 col(2) = irow - 1 598 col(3) = irow 599 col(4) = irow + 1 600 col(5) = irow + solver%mx 601 row(1) = irow 602 PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr)) 603 endif 604 20 continue 605 606 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 607 608 return 609 end subroutine FormJacobianLocal 610 611! --------------------------------------------------------------------- 612! 613! FormFunction - Evaluates nonlinear function, F(x). 614! 615! Input Parameters: 616! snes - the SNES context 617! X - input vector 618! dummy - optional user-defined context, as set by SNESSetFunction() 619! (not used here) 620! 621! Output Parameter: 622! F - function vector 623! 624 subroutine FormFunction(snesIn,X,F,solver,ierr) 625#include <petsc/finclude/petscsnes.h> 626 use petscsnes 627 use ex73f90tmodule 628 implicit none 629! Input/output variables: 630 SNES:: snesIn 631 Vec:: X,F 632 PetscErrorCode ierr 633 type (ex73f90tmodule_type) solver 634 635! Declarations for use with local arrays: 636 Vec:: Xsub(2),Fsub(2) 637 PetscInt itwo 638 639! Scatter ghost points to local vector, using the 2-step process 640! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 641! By placing code between these two statements, computations can 642! be done while messages are in transition. 643 644 itwo = 2 645 PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 646 PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 647 648 PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr)) 649 PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr)) 650 651! do rest of operator (linear) 652 PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr)) 653 PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr)) 654 PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr)) 655 656 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 657 PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 658 return 659 end subroutine formfunction 660 661! --------------------------------------------------------------------- 662! 663! FormFunctionNLTerm - Computes nonlinear function, called by 664! the higher level routine FormFunction(). 665! 666! Input Parameter: 667! x - local vector data 668! 669! Output Parameters: 670! f - local vector data, f(x) 671! ierr - error code 672! 673! Notes: 674! This routine uses standard Fortran-style computations over a 2-dim array. 675! 676 subroutine FormFunctionNLTerm(X1,F1,solver,ierr) 677#include <petsc/finclude/petscvec.h> 678 use petscvec 679 use ex73f90tmodule 680 implicit none 681! Input/output variables: 682 type (ex73f90tmodule_type) solver 683 Vec:: X1,F1 684 PetscErrorCode ierr 685! Local variables: 686 PetscScalar sc 687 PetscScalar u,v(1) 688 PetscInt i,j,low,high,ii,ione,irow,row(1) 689 PetscScalar,pointer :: lx_v(:) 690 691 sc = solver%lambda 692 ione = 1 693 694 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 695 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 696 697! Compute function over the locally owned part of the grid 698 ii = 0 699 do 20 irow=low,high-1 700 j = irow/solver%mx 701 i = mod(irow,solver%mx) 702 ii = ii + 1 ! one based local index 703 row(1) = irow 704 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 705 v(1) = 0.0 706 else 707 u = lx_v(ii) 708 v(1) = -sc*exp(u) 709 endif 710 PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr)) 711 20 continue 712 713 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 714 715 PetscCall(VecAssemblyBegin(F1,ierr)) 716 PetscCall(VecAssemblyEnd(F1,ierr)) 717 718 ierr = 0 719 return 720 end subroutine FormFunctionNLTerm 721 722!/*TEST 723! 724! build: 725! requires: !single !complex 726! 727! test: 728! nsize: 4 729! args: -par 5.0 -da_grid_x 10 -da_grid_y 10 -snes_monitor_short -snes_linesearch_type basic -snes_converged_reason -ksp_type fgmres -ksp_norm_type unpreconditioned -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type upper -ksp_monitor_short -fieldsplit_lambda_ksp_type preonly -fieldsplit_lambda_pc_type jacobi -fieldsplit_phi_pc_type gamg -fieldsplit_phi_pc_gamg_esteig_ksp_type cg -fieldsplit_phi_pc_gamg_esteig_ksp_max_it 10 -fieldsplit_phi_pc_gamg_agg_nsmooths 1 -fieldsplit_phi_pc_gamg_threshold 0. 730! 731!TEST*/ 732