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 PetscCheckA(solver%lambda .le. lambda_max .and. solver%lambda .ge. lambda_min,PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range') 139 140! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141! Create vector data structures; set function evaluation routine 142! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 144! just get size 145 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)) 146 PetscCallA(DMSetFromOptions(daphi,ierr)) 147 PetscCallA(DMSetUp(daphi,ierr)) 148 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)) 149 N1 = solver%my*solver%mx 150 N2 = solver%my 151 flg = .false. 152 PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr)) 153 if (flg) then 154 N2 = 0 155 endif 156 157 PetscCallA(DMDestroy(daphi,ierr)) 158 159! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160! Create matrix data structure; set Jacobian evaluation routine 161! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 162 PetscCallA(DMShellCreate(PETSC_COMM_WORLD,daphi,ierr)) 163 PetscCallA(DMSetOptionsPrefix(daphi,'phi_',ierr)) 164 PetscCallA(DMSetFromOptions(daphi,ierr)) 165 166 PetscCallA(VecCreate(PETSC_COMM_WORLD,x1,ierr)) 167 PetscCallA(VecSetSizes(x1,PETSC_DECIDE,N1,ierr)) 168 PetscCallA(VecSetFromOptions(x1,ierr)) 169 170 PetscCallA(VecGetOwnershipRange(x1,low,high,ierr)) 171 nloc = high - low 172 173 PetscCallA(MatCreate(PETSC_COMM_WORLD,Amat,ierr)) 174 PetscCallA(MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr)) 175 PetscCallA(MatSetUp(Amat,ierr)) 176 177 PetscCallA(MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr)) 178 PetscCallA(MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr)) 179 PetscCallA(MatSetUp(solver%AmatLin,ierr)) 180 181 PetscCallA(FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr)) 182 PetscCallA(MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr)) 183 PetscCallA(MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr)) 184 185 PetscCallA(DMShellSetGlobalVector(daphi,x1,ierr)) 186 PetscCallA(DMShellSetMatrix(daphi,Amat,ierr)) 187 188 PetscCallA(VecCreate(PETSC_COMM_SELF,x1loc,ierr)) 189 PetscCallA(VecSetSizes(x1loc,nloc,nloc,ierr)) 190 PetscCallA(VecSetFromOptions(x1loc,ierr)) 191 PetscCallA(DMShellSetLocalVector(daphi,x1loc,ierr)) 192 193! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194! Create B, C, & D matrices 195! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196 PetscCallA(MatCreate(PETSC_COMM_WORLD,Cmat,ierr)) 197 PetscCallA(MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr)) 198 PetscCallA(MatSetUp(Cmat,ierr)) 199! create data for C and B 200 PetscCallA(MatCreate(PETSC_COMM_WORLD,Bmat,ierr)) 201 PetscCallA(MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr)) 202 PetscCallA(MatSetUp(Bmat,ierr)) 203! create data for D 204 PetscCallA(MatCreate(PETSC_COMM_WORLD,Dmat,ierr)) 205 PetscCallA(MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr)) 206 PetscCallA(MatSetUp(Dmat,ierr)) 207 208 PetscCallA(VecCreate(PETSC_COMM_WORLD,x2,ierr)) 209 PetscCallA(VecSetSizes(x2,PETSC_DECIDE,N2,ierr)) 210 PetscCallA(VecSetFromOptions(x2,ierr)) 211 212 PetscCallA(VecGetOwnershipRange(x2,lamlow,lamhigh,ierr)) 213 nloclam = lamhigh-lamlow 214 215! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 216! Set fake B and C 217! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 218 one = 1.0 219 if (N2 .gt. 0) then 220 bval(1) = -one/(solver%mx-2) 221! cval = -one/(solver%my*solver%mx) 222 cval(1) = -one 223 do 20 irow=low,high-1 224 j = irow/solver%mx ! row in domain 225 i = mod(irow,solver%mx) 226 row(1) = irow 227 col(1) = j 228 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 229 ! no op 230 else 231 PetscCallA(MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr)) 232 endif 233 row(1) = j 234 PetscCallA(MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr)) 235 20 continue 236 endif 237 PetscCallA(MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr)) 238 PetscCallA(MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr)) 239 PetscCallA(MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr)) 240 PetscCallA(MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr)) 241 242! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243! Set D (identity) 244! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 245 do 30 j=lamlow,lamhigh-1 246 row(1) = j 247 cval(1) = one 248 PetscCallA(MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr)) 249 30 continue 250 PetscCallA(MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr)) 251 PetscCallA(MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr)) 252 253! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 254! DM for lambda (dalam) : temp driver for A block, setup A block solver data 255! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 256 PetscCallA(DMShellCreate(PETSC_COMM_WORLD,dalam,ierr)) 257 PetscCallA(DMShellSetGlobalVector(dalam,x2,ierr)) 258 PetscCallA(DMShellSetMatrix(dalam,Dmat,ierr)) 259 260 PetscCallA(VecCreate(PETSC_COMM_SELF,x2loc,ierr)) 261 PetscCallA(VecSetSizes(x2loc,nloclam,nloclam,ierr)) 262 PetscCallA(VecSetFromOptions(x2loc,ierr)) 263 PetscCallA(DMShellSetLocalVector(dalam,x2loc,ierr)) 264 265 PetscCallA(DMSetOptionsPrefix(dalam,'lambda_',ierr)) 266 PetscCallA(DMSetFromOptions(dalam,ierr)) 267! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 268! Create field split DA 269! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 270 PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr)) 271 PetscCallA(DMCompositeAddDM(solver%da,daphi,ierr)) 272 PetscCallA(DMCompositeAddDM(solver%da,dalam,ierr)) 273 PetscCallA(DMSetFromOptions(solver%da,ierr)) 274 PetscCallA(DMSetUp(solver%da,ierr)) 275 PetscCallA(DMCompositeGetGlobalISs(solver%da,isglobal,ierr)) 276 solver%isPhi = isglobal(1) 277 solver%isLambda = isglobal(2) 278 279! cache matrices 280 solver%Amat = Amat 281 solver%Bmat = Bmat 282 solver%Cmat = Cmat 283 solver%Dmat = Dmat 284 285 matArray(1) = Amat 286 matArray(2) = Bmat 287 matArray(3) = Cmat 288 matArray(4) = Dmat 289 290 PetscCallA(MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr)) 291 PetscCallA(MatSetFromOptions(KKTmat,ierr)) 292 293! Extract global and local vectors from DMDA; then duplicate for remaining 294! vectors that are the same types 295 PetscCallA(MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr)) 296 PetscCallA(VecDuplicate(x,r,ierr)) 297 298 PetscCallA(SNESCreate(PETSC_COMM_WORLD,mysnes,ierr)) 299 300 PetscCallA(SNESSetDM(mysnes,solver%da,ierr)) 301 302 PetscCallA(SNESSetApplicationContext(mysnes,solver,ierr)) 303 304 PetscCallA(SNESSetDM(mysnes,solver%da,ierr)) 305 306! Set function evaluation routine and vector 307 PetscCallA(SNESSetFunction(mysnes,r,FormFunction,solver,ierr)) 308 309 PetscCallA(SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr)) 310 311! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 312! Customize nonlinear solver; set runtime options 313! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 314! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 315 PetscCallA(SNESSetFromOptions(mysnes,ierr)) 316 317! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 318! Evaluate initial guess; then solve nonlinear system. 319! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 320! Note: The user should initialize the vector, x, with the initial guess 321! for the nonlinear solver prior to calling SNESSolve(). In particular, 322! to employ an initial guess of zero, the user should explicitly set 323! this vector to zero by calling VecSet(). 324 325 PetscCallA(FormInitialGuess(mysnes,x,ierr)) 326 PetscCallA(SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr)) 327 PetscCallA(SNESGetIterationNumber(mysnes,its,ierr)) 328 if (solver%rank .eq. 0) then 329 write(6,100) its 330 endif 331 100 format('Number of SNES iterations = ',i5) 332 333! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 334! Free work space. All PETSc objects should be destroyed when they 335! are no longer needed. 336! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 337 PetscCallA(MatDestroy(KKTmat,ierr)) 338 PetscCallA(MatDestroy(Amat,ierr)) 339 PetscCallA(MatDestroy(Dmat,ierr)) 340 PetscCallA(MatDestroy(Bmat,ierr)) 341 PetscCallA(MatDestroy(Cmat,ierr)) 342 PetscCallA(MatDestroy(solver%AmatLin,ierr)) 343 PetscCallA(ISDestroy(solver%isPhi,ierr)) 344 PetscCallA(ISDestroy(solver%isLambda,ierr)) 345 PetscCallA(VecDestroy(x,ierr)) 346 PetscCallA(VecDestroy(x2,ierr)) 347 PetscCallA(VecDestroy(x1,ierr)) 348 PetscCallA(VecDestroy(x1loc,ierr)) 349 PetscCallA(VecDestroy(x2loc,ierr)) 350 PetscCallA(VecDestroy(r,ierr)) 351 PetscCallA(SNESDestroy(mysnes,ierr)) 352 PetscCallA(DMDestroy(solver%da,ierr)) 353 PetscCallA(DMDestroy(daphi,ierr)) 354 PetscCallA(DMDestroy(dalam,ierr)) 355 356 PetscCallA(PetscFinalize(ierr)) 357 end 358 359! --------------------------------------------------------------------- 360! 361! FormInitialGuess - Forms initial approximation. 362! 363! Input Parameters: 364! X - vector 365! 366! Output Parameter: 367! X - vector 368! 369! Notes: 370! This routine serves as a wrapper for the lower-level routine 371! "InitialGuessLocal", where the actual computations are 372! done using the standard Fortran style of treating the local 373! vector data as a multidimensional array over the local mesh. 374! This routine merely handles ghost point scatters and accesses 375! the local vector data via VecGetArrayF90() and VecRestoreArrayF90(). 376! 377 subroutine FormInitialGuess(mysnes,Xnest,ierr) 378#include <petsc/finclude/petscsnes.h> 379 use petscsnes 380 use ex73f90tmodule 381 use ex73f90tmodule_interfaces 382 implicit none 383! Input/output variables: 384 SNES:: mysnes 385 Vec:: Xnest 386 PetscErrorCode ierr 387 388! Declarations for use with local arrays: 389 type(ex73f90tmodule_type), pointer:: solver 390 Vec:: Xsub(2) 391 PetscInt:: izero,ione,itwo 392 393 izero = 0 394 ione = 1 395 itwo = 2 396 ierr = 0 397 PetscCall(SNESGetApplicationContext(mysnes,solver,ierr)) 398 PetscCall(DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 399 400 PetscCall(InitialGuessLocal(solver,Xsub(1),ierr)) 401 PetscCall(VecAssemblyBegin(Xsub(1),ierr)) 402 PetscCall(VecAssemblyEnd(Xsub(1),ierr)) 403 404! zero out lambda 405 PetscCall(VecZeroEntries(Xsub(2),ierr)) 406 PetscCall(DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 407 408 end subroutine FormInitialGuess 409 410! --------------------------------------------------------------------- 411! 412! InitialGuessLocal - Computes initial approximation, called by 413! the higher level routine FormInitialGuess(). 414! 415! Input Parameter: 416! X1 - local vector data 417! 418! Output Parameters: 419! x - local vector data 420! ierr - error code 421! 422! Notes: 423! This routine uses standard Fortran-style computations over a 2-dim array. 424! 425 subroutine InitialGuessLocal(solver,X1,ierr) 426#include <petsc/finclude/petscsys.h> 427 use petscsys 428 use ex73f90tmodule 429 implicit none 430! Input/output variables: 431 type (ex73f90tmodule_type) solver 432 Vec:: X1 433 PetscErrorCode ierr 434 435! Local variables: 436 PetscInt row,i,j,ione,low,high 437 PetscReal temp1,temp,hx,hy,v 438 PetscReal one 439 440! Set parameters 441 ione = 1 442 ierr = 0 443 one = 1.0 444 hx = one/(solver%mx-1) 445 hy = one/(solver%my-1) 446 temp1 = solver%lambda/(solver%lambda + one) + one 447 448 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 449 450 do 20 row=low,high-1 451 j = row/solver%mx 452 i = mod(row,solver%mx) 453 temp = min(j,solver%my-j+1)*hy 454 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 455 v = 0.0 456 else 457 v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp)) 458 endif 459 PetscCall(VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr)) 460 20 continue 461 462 end subroutine InitialGuessLocal 463 464! --------------------------------------------------------------------- 465! 466! FormJacobian - Evaluates Jacobian matrix. 467! 468! Input Parameters: 469! dummy - the SNES context 470! x - input vector 471! solver - solver data 472! 473! Output Parameters: 474! jac - Jacobian matrix 475! jac_prec - optionally different preconditioning matrix (not used here) 476! flag - flag indicating matrix structure 477! 478 subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr) 479#include <petsc/finclude/petscsnes.h> 480 use petscsnes 481 use ex73f90tmodule 482 implicit none 483! Input/output variables: 484 SNES:: dummy 485 Vec:: X 486 Mat:: jac,jac_prec 487 type(ex73f90tmodule_type) solver 488 PetscErrorCode ierr 489 490! Declarations for use with local arrays: 491 Vec:: Xsub(1) 492 Mat:: Amat 493 PetscInt ione 494 495 ione = 1 496 497 PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 498 499! Compute entries for the locally owned part of the Jacobian preconditioner. 500 PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr)) 501 502 PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr)) 503 PetscCall(MatDestroy(Amat,ierr)) ! discard our reference 504 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 505 506 ! the rest of the matrix is not touched 507 PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 508 PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 509 if (jac .ne. jac_prec) then 510 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 511 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 512 end if 513 514! Tell the matrix we will never add a new nonzero location to the 515! matrix. If we do it will generate an error. 516 PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 517 518 end subroutine FormJacobian 519 520! --------------------------------------------------------------------- 521! 522! FormJacobianLocal - Computes Jacobian preconditioner matrix, 523! called by the higher level routine FormJacobian(). 524! 525! Input Parameters: 526! x - local vector data 527! 528! Output Parameters: 529! jac - Jacobian preconditioner matrix 530! ierr - error code 531! 532! Notes: 533! This routine uses standard Fortran-style computations over a 2-dim array. 534! 535 subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr) 536#include <petsc/finclude/petscmat.h> 537 use petscmat 538 use ex73f90tmodule 539 implicit none 540! Input/output variables: 541 type (ex73f90tmodule_type) solver 542 Vec:: X1 543 Mat:: jac 544 logical add_nl_term 545 PetscErrorCode ierr 546 547! Local variables: 548 PetscInt irow,row(1),col(5),i,j 549 PetscInt ione,ifive,low,high,ii 550 PetscScalar two,one,hx,hy,hy2inv 551 PetscScalar hx2inv,sc,v(5) 552 PetscScalar,pointer :: lx_v(:) 553 554! Set parameters 555 ione = 1 556 ifive = 5 557 one = 1.0 558 two = 2.0 559 hx = one/(solver%mx-1) 560 hy = one/(solver%my-1) 561 sc = solver%lambda 562 hx2inv = one/(hx*hx) 563 hy2inv = one/(hy*hy) 564 565 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 566 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 567 568 ii = 0 569 do 20 irow=low,high-1 570 j = irow/solver%mx 571 i = mod(irow,solver%mx) 572 ii = ii + 1 ! one based local index 573! boundary points 574 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 575 col(1) = irow 576 row(1) = irow 577 v(1) = one 578 PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr)) 579! interior grid points 580 else 581 v(1) = -hy2inv 582 if (j-1==0) v(1) = 0.0 583 v(2) = -hx2inv 584 if (i-1==0) v(2) = 0.0 585 v(3) = two*(hx2inv + hy2inv) 586 if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii)) 587 v(4) = -hx2inv 588 if (i+1==solver%mx-1) v(4) = 0.0 589 v(5) = -hy2inv 590 if (j+1==solver%my-1) v(5) = 0.0 591 col(1) = irow - solver%mx 592 col(2) = irow - 1 593 col(3) = irow 594 col(4) = irow + 1 595 col(5) = irow + solver%mx 596 row(1) = irow 597 PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr)) 598 endif 599 20 continue 600 601 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 602 603 end subroutine FormJacobianLocal 604 605! --------------------------------------------------------------------- 606! 607! FormFunction - Evaluates nonlinear function, F(x). 608! 609! Input Parameters: 610! snes - the SNES context 611! X - input vector 612! dummy - optional user-defined context, as set by SNESSetFunction() 613! (not used here) 614! 615! Output Parameter: 616! F - function vector 617! 618 subroutine FormFunction(snesIn,X,F,solver,ierr) 619#include <petsc/finclude/petscsnes.h> 620 use petscsnes 621 use ex73f90tmodule 622 implicit none 623! Input/output variables: 624 SNES:: snesIn 625 Vec:: X,F 626 PetscErrorCode ierr 627 type (ex73f90tmodule_type) solver 628 629! Declarations for use with local arrays: 630 Vec:: Xsub(2),Fsub(2) 631 PetscInt itwo 632 633! Scatter ghost points to local vector, using the 2-step process 634! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 635! By placing code between these two statements, computations can 636! be done while messages are in transition. 637 638 itwo = 2 639 PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 640 PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 641 642 PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr)) 643 PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr)) 644 645! do rest of operator (linear) 646 PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr)) 647 PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr)) 648 PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr)) 649 650 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 651 PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 652 end subroutine formfunction 653 654! --------------------------------------------------------------------- 655! 656! FormFunctionNLTerm - Computes nonlinear function, called by 657! the higher level routine FormFunction(). 658! 659! Input Parameter: 660! x - local vector data 661! 662! Output Parameters: 663! f - local vector data, f(x) 664! ierr - error code 665! 666! Notes: 667! This routine uses standard Fortran-style computations over a 2-dim array. 668! 669 subroutine FormFunctionNLTerm(X1,F1,solver,ierr) 670#include <petsc/finclude/petscvec.h> 671 use petscvec 672 use ex73f90tmodule 673 implicit none 674! Input/output variables: 675 type (ex73f90tmodule_type) solver 676 Vec:: X1,F1 677 PetscErrorCode ierr 678! Local variables: 679 PetscScalar sc 680 PetscScalar u,v(1) 681 PetscInt i,j,low,high,ii,ione,irow,row(1) 682 PetscScalar,pointer :: lx_v(:) 683 684 sc = solver%lambda 685 ione = 1 686 687 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 688 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 689 690! Compute function over the locally owned part of the grid 691 ii = 0 692 do 20 irow=low,high-1 693 j = irow/solver%mx 694 i = mod(irow,solver%mx) 695 ii = ii + 1 ! one based local index 696 row(1) = irow 697 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 698 v(1) = 0.0 699 else 700 u = lx_v(ii) 701 v(1) = -sc*exp(u) 702 endif 703 PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr)) 704 20 continue 705 706 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 707 708 PetscCall(VecAssemblyBegin(F1,ierr)) 709 PetscCall(VecAssemblyEnd(F1,ierr)) 710 711 ierr = 0 712 end subroutine FormFunctionNLTerm 713 714!/*TEST 715! 716! build: 717! requires: !single !complex 718! 719! test: 720! nsize: 4 721! 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. 722! 723!TEST*/ 724