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 return 409 end subroutine FormInitialGuess 410 411! --------------------------------------------------------------------- 412! 413! InitialGuessLocal - Computes initial approximation, called by 414! the higher level routine FormInitialGuess(). 415! 416! Input Parameter: 417! X1 - local vector data 418! 419! Output Parameters: 420! x - local vector data 421! ierr - error code 422! 423! Notes: 424! This routine uses standard Fortran-style computations over a 2-dim array. 425! 426 subroutine InitialGuessLocal(solver,X1,ierr) 427#include <petsc/finclude/petscsys.h> 428 use petscsys 429 use ex73f90tmodule 430 implicit none 431! Input/output variables: 432 type (ex73f90tmodule_type) solver 433 Vec:: X1 434 PetscErrorCode ierr 435 436! Local variables: 437 PetscInt row,i,j,ione,low,high 438 PetscReal temp1,temp,hx,hy,v 439 PetscReal one 440 441! Set parameters 442 ione = 1 443 ierr = 0 444 one = 1.0 445 hx = one/(solver%mx-1) 446 hy = one/(solver%my-1) 447 temp1 = solver%lambda/(solver%lambda + one) + one 448 449 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 450 451 do 20 row=low,high-1 452 j = row/solver%mx 453 i = mod(row,solver%mx) 454 temp = min(j,solver%my-j+1)*hy 455 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 456 v = 0.0 457 else 458 v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp)) 459 endif 460 PetscCall(VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr)) 461 20 continue 462 463 return 464 end subroutine InitialGuessLocal 465 466! --------------------------------------------------------------------- 467! 468! FormJacobian - Evaluates Jacobian matrix. 469! 470! Input Parameters: 471! dummy - the SNES context 472! x - input vector 473! solver - solver data 474! 475! Output Parameters: 476! jac - Jacobian matrix 477! jac_prec - optionally different preconditioning matrix (not used here) 478! flag - flag indicating matrix structure 479! 480 subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr) 481#include <petsc/finclude/petscsnes.h> 482 use petscsnes 483 use ex73f90tmodule 484 implicit none 485! Input/output variables: 486 SNES:: dummy 487 Vec:: X 488 Mat:: jac,jac_prec 489 type(ex73f90tmodule_type) solver 490 PetscErrorCode ierr 491 492! Declarations for use with local arrays: 493 Vec:: Xsub(1) 494 Mat:: Amat 495 PetscInt ione 496 497 ione = 1 498 499 PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 500 501! Compute entries for the locally owned part of the Jacobian preconditioner. 502 PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr)) 503 504 PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr)) 505 PetscCall(MatDestroy(Amat,ierr)) ! discard our reference 506 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 507 508 ! the rest of the matrix is not touched 509 PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 510 PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 511 if (jac .ne. jac_prec) then 512 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 513 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 514 end if 515 516! Tell the matrix we will never add a new nonzero location to the 517! matrix. If we do it will generate an error. 518 PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 519 520 return 521 end subroutine FormJacobian 522 523! --------------------------------------------------------------------- 524! 525! FormJacobianLocal - Computes Jacobian preconditioner matrix, 526! called by the higher level routine FormJacobian(). 527! 528! Input Parameters: 529! x - local vector data 530! 531! Output Parameters: 532! jac - Jacobian preconditioner matrix 533! ierr - error code 534! 535! Notes: 536! This routine uses standard Fortran-style computations over a 2-dim array. 537! 538 subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr) 539#include <petsc/finclude/petscmat.h> 540 use petscmat 541 use ex73f90tmodule 542 implicit none 543! Input/output variables: 544 type (ex73f90tmodule_type) solver 545 Vec:: X1 546 Mat:: jac 547 logical add_nl_term 548 PetscErrorCode ierr 549 550! Local variables: 551 PetscInt irow,row(1),col(5),i,j 552 PetscInt ione,ifive,low,high,ii 553 PetscScalar two,one,hx,hy,hy2inv 554 PetscScalar hx2inv,sc,v(5) 555 PetscScalar,pointer :: lx_v(:) 556 557! Set parameters 558 ione = 1 559 ifive = 5 560 one = 1.0 561 two = 2.0 562 hx = one/(solver%mx-1) 563 hy = one/(solver%my-1) 564 sc = solver%lambda 565 hx2inv = one/(hx*hx) 566 hy2inv = one/(hy*hy) 567 568 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 569 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 570 571 ii = 0 572 do 20 irow=low,high-1 573 j = irow/solver%mx 574 i = mod(irow,solver%mx) 575 ii = ii + 1 ! one based local index 576! boundary points 577 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 578 col(1) = irow 579 row(1) = irow 580 v(1) = one 581 PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr)) 582! interior grid points 583 else 584 v(1) = -hy2inv 585 if (j-1==0) v(1) = 0.0 586 v(2) = -hx2inv 587 if (i-1==0) v(2) = 0.0 588 v(3) = two*(hx2inv + hy2inv) 589 if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii)) 590 v(4) = -hx2inv 591 if (i+1==solver%mx-1) v(4) = 0.0 592 v(5) = -hy2inv 593 if (j+1==solver%my-1) v(5) = 0.0 594 col(1) = irow - solver%mx 595 col(2) = irow - 1 596 col(3) = irow 597 col(4) = irow + 1 598 col(5) = irow + solver%mx 599 row(1) = irow 600 PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr)) 601 endif 602 20 continue 603 604 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 605 606 return 607 end subroutine FormJacobianLocal 608 609! --------------------------------------------------------------------- 610! 611! FormFunction - Evaluates nonlinear function, F(x). 612! 613! Input Parameters: 614! snes - the SNES context 615! X - input vector 616! dummy - optional user-defined context, as set by SNESSetFunction() 617! (not used here) 618! 619! Output Parameter: 620! F - function vector 621! 622 subroutine FormFunction(snesIn,X,F,solver,ierr) 623#include <petsc/finclude/petscsnes.h> 624 use petscsnes 625 use ex73f90tmodule 626 implicit none 627! Input/output variables: 628 SNES:: snesIn 629 Vec:: X,F 630 PetscErrorCode ierr 631 type (ex73f90tmodule_type) solver 632 633! Declarations for use with local arrays: 634 Vec:: Xsub(2),Fsub(2) 635 PetscInt itwo 636 637! Scatter ghost points to local vector, using the 2-step process 638! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 639! By placing code between these two statements, computations can 640! be done while messages are in transition. 641 642 itwo = 2 643 PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 644 PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 645 646 PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr)) 647 PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr)) 648 649! do rest of operator (linear) 650 PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr)) 651 PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr)) 652 PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr)) 653 654 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 655 PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 656 return 657 end subroutine formfunction 658 659! --------------------------------------------------------------------- 660! 661! FormFunctionNLTerm - Computes nonlinear function, called by 662! the higher level routine FormFunction(). 663! 664! Input Parameter: 665! x - local vector data 666! 667! Output Parameters: 668! f - local vector data, f(x) 669! ierr - error code 670! 671! Notes: 672! This routine uses standard Fortran-style computations over a 2-dim array. 673! 674 subroutine FormFunctionNLTerm(X1,F1,solver,ierr) 675#include <petsc/finclude/petscvec.h> 676 use petscvec 677 use ex73f90tmodule 678 implicit none 679! Input/output variables: 680 type (ex73f90tmodule_type) solver 681 Vec:: X1,F1 682 PetscErrorCode ierr 683! Local variables: 684 PetscScalar sc 685 PetscScalar u,v(1) 686 PetscInt i,j,low,high,ii,ione,irow,row(1) 687 PetscScalar,pointer :: lx_v(:) 688 689 sc = solver%lambda 690 ione = 1 691 692 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 693 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 694 695! Compute function over the locally owned part of the grid 696 ii = 0 697 do 20 irow=low,high-1 698 j = irow/solver%mx 699 i = mod(irow,solver%mx) 700 ii = ii + 1 ! one based local index 701 row(1) = irow 702 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 703 v(1) = 0.0 704 else 705 u = lx_v(ii) 706 v(1) = -sc*exp(u) 707 endif 708 PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr)) 709 20 continue 710 711 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 712 713 PetscCall(VecAssemblyBegin(F1,ierr)) 714 PetscCall(VecAssemblyEnd(F1,ierr)) 715 716 ierr = 0 717 return 718 end subroutine FormFunctionNLTerm 719 720!/*TEST 721! 722! build: 723! requires: !single !complex 724! 725! test: 726! nsize: 4 727! 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. 728! 729!TEST*/ 730