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