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,PETSC_NULL_INTEGER,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_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,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,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,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! flag - flag indicating matrix structure 500! 501 subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr) 502#include <petsc/finclude/petscsnes.h> 503 use petscsnes 504 use ex73f90tmodule 505 implicit none 506! Input/output variables: 507 SNES:: dummy 508 Vec:: X 509 Mat:: jac,jac_prec 510 type(ex73f90tmodule_type) solver 511 PetscErrorCode ierr 512 513! Declarations for use with local arrays: 514 Vec:: Xsub(1) 515 Mat:: Amat 516 PetscInt ione 517 518 ione = 1 519 520 PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 521 522! Compute entries for the locally owned part of the Jacobian preconditioner. 523 PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr)) 524 525 PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr)) 526 PetscCall(MatDestroy(Amat,ierr)) ! discard our reference 527 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr)) 528 529 ! the rest of the matrix is not touched 530 PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 531 PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 532 if (jac .ne. jac_prec) then 533 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 534 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 535 end if 536 537! Tell the matrix we will never add a new nonzero location to the 538! matrix. If we do it will generate an error. 539 PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 540 541 end subroutine FormJacobian 542 543! --------------------------------------------------------------------- 544! 545! FormJacobianLocal - Computes Jacobian preconditioner matrix, 546! called by the higher level routine FormJacobian(). 547! 548! Input Parameters: 549! x - local vector data 550! 551! Output Parameters: 552! jac - Jacobian preconditioner matrix 553! ierr - error code 554! 555! Notes: 556! This routine uses standard Fortran-style computations over a 2-dim array. 557! 558 subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr) 559#include <petsc/finclude/petscmat.h> 560 use petscmat 561 use ex73f90tmodule 562 implicit none 563! Input/output variables: 564 type (ex73f90tmodule_type) solver 565 Vec:: X1 566 Mat:: jac 567 logical add_nl_term 568 PetscErrorCode ierr 569 570! Local variables: 571 PetscInt irow,row(1),col(5),i,j 572 PetscInt ione,ifive,low,high,ii 573 PetscScalar two,one,hx,hy,hy2inv 574 PetscScalar hx2inv,sc,v(5) 575 PetscScalar,pointer :: lx_v(:) 576 577! Set parameters 578 ione = 1 579 ifive = 5 580 one = 1.0 581 two = 2.0 582 hx = one/(solver%mx-1) 583 hy = one/(solver%my-1) 584 sc = solver%lambda 585 hx2inv = one/(hx*hx) 586 hy2inv = one/(hy*hy) 587 588 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 589 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 590 591 ii = 0 592 do 20 irow=low,high-1 593 j = irow/solver%mx 594 i = mod(irow,solver%mx) 595 ii = ii + 1 ! one based local index 596! boundary points 597 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 598 col(1) = irow 599 row(1) = irow 600 v(1) = one 601 PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr)) 602! interior grid points 603 else 604 v(1) = -hy2inv 605 if (j-1==0) v(1) = 0.0 606 v(2) = -hx2inv 607 if (i-1==0) v(2) = 0.0 608 v(3) = two*(hx2inv + hy2inv) 609 if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii)) 610 v(4) = -hx2inv 611 if (i+1==solver%mx-1) v(4) = 0.0 612 v(5) = -hy2inv 613 if (j+1==solver%my-1) v(5) = 0.0 614 col(1) = irow - solver%mx 615 col(2) = irow - 1 616 col(3) = irow 617 col(4) = irow + 1 618 col(5) = irow + solver%mx 619 row(1) = irow 620 PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr)) 621 endif 622 20 continue 623 624 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 625 626 end subroutine FormJacobianLocal 627 628! --------------------------------------------------------------------- 629! 630! FormFunction - Evaluates nonlinear function, F(x). 631! 632! Input Parameters: 633! snes - the SNES context 634! X - input vector 635! dummy - optional user-defined context, as set by SNESSetFunction() 636! (not used here) 637! 638! Output Parameter: 639! F - function vector 640! 641 subroutine FormFunction(snesIn,X,F,solver,ierr) 642#include <petsc/finclude/petscsnes.h> 643 use petscsnes 644 use ex73f90tmodule 645 implicit none 646! Input/output variables: 647 SNES:: snesIn 648 Vec:: X,F 649 PetscErrorCode ierr 650 type (ex73f90tmodule_type) solver 651 652! Declarations for use with local arrays: 653 Vec:: Xsub(2),Fsub(2) 654 PetscInt itwo 655 656! Scatter ghost points to local vector, using the 2-step process 657! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 658! By placing code between these two statements, computations can 659! be done while messages are in transition. 660 661 itwo = 2 662 PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 663 PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 664 665 PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr)) 666 PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr)) 667 668! do rest of operator (linear) 669 PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr)) 670 PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr)) 671 PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr)) 672 673 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr)) 674 PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr)) 675 end subroutine formfunction 676 677! --------------------------------------------------------------------- 678! 679! FormFunctionNLTerm - Computes nonlinear function, called by 680! the higher level routine FormFunction(). 681! 682! Input Parameter: 683! x - local vector data 684! 685! Output Parameters: 686! f - local vector data, f(x) 687! ierr - error code 688! 689! Notes: 690! This routine uses standard Fortran-style computations over a 2-dim array. 691! 692 subroutine FormFunctionNLTerm(X1,F1,solver,ierr) 693#include <petsc/finclude/petscvec.h> 694 use petscvec 695 use ex73f90tmodule 696 implicit none 697! Input/output variables: 698 type (ex73f90tmodule_type) solver 699 Vec:: X1,F1 700 PetscErrorCode ierr 701! Local variables: 702 PetscScalar sc 703 PetscScalar u,v(1) 704 PetscInt i,j,low,high,ii,ione,irow,row(1) 705 PetscScalar,pointer :: lx_v(:) 706 707 sc = solver%lambda 708 ione = 1 709 710 PetscCall(VecGetArrayReadF90(X1,lx_v,ierr)) 711 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 712 713! Compute function over the locally owned part of the grid 714 ii = 0 715 do 20 irow=low,high-1 716 j = irow/solver%mx 717 i = mod(irow,solver%mx) 718 ii = ii + 1 ! one based local index 719 row(1) = irow 720 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 721 v(1) = 0.0 722 else 723 u = lx_v(ii) 724 v(1) = -sc*exp(u) 725 endif 726 PetscCall(VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr)) 727 20 continue 728 729 PetscCall(VecRestoreArrayReadF90(X1,lx_v,ierr)) 730 731 PetscCall(VecAssemblyBegin(F1,ierr)) 732 PetscCall(VecAssemblyEnd(F1,ierr)) 733 734 ierr = 0 735 end subroutine FormFunctionNLTerm 736 737!/*TEST 738! 739! build: 740! requires: !single !complex 741! 742! test: 743! nsize: 4 744! 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. 745! 746! test: 747! args: -snes_linesearch_type {{l2 cp}separate output} -objective {{false true}shared output} 748! 749! test: 750! args: -snes_linesearch_type bt -objective {{false true}separate output} 751! 752!TEST*/ 753