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/petscdmda.h> 43#include <petsc/finclude/petscdmcomposite.h> 44#include <petsc/finclude/petscmat.h> 45 use petscdmda 46 use petscdmcomposite 47 use petscmat 48 type ex73f90tmodule_type 49 DM::da 50! temp A block stuff 51 PetscInt mx,my 52 PetscMPIInt rank 53 PetscReal lambda 54! Mats 55 Mat::Amat,AmatLin,Bmat,CMat,Dmat 56 IS::isPhi,isLambda 57 end type ex73f90tmodule_type 58 59 end module ex73f90tmodule 60 61 module ex73f90tmodule_interfaces 62 use ex73f90tmodule 63 64 Interface SNESSetApplicationContext 65 Subroutine SNESSetApplicationContext(snesIn,ctx,ierr) 66#include <petsc/finclude/petscsnes.h> 67 use petscsnes 68 use ex73f90tmodule 69 SNES:: snesIn 70 type(ex73f90tmodule_type) ctx 71 PetscErrorCode ierr 72 End Subroutine 73 End Interface SNESSetApplicationContext 74 75 Interface SNESGetApplicationContext 76 Subroutine SNESGetApplicationContext(snesIn,ctx,ierr) 77#include <petsc/finclude/petscsnes.h> 78 use petscsnes 79 use ex73f90tmodule 80 SNES:: snesIn 81 type(ex73f90tmodule_type), pointer :: ctx 82 PetscErrorCode ierr 83 End Subroutine 84 End Interface SNESGetApplicationContext 85 end module ex73f90tmodule_interfaces 86 87 subroutine MyObjective(snes, x, result, ctx, ierr ) 88#include <petsc/finclude/petsc.h> 89 use petsc 90 implicit none 91 PetscInt ctx 92 Vec x, f 93 SNES snes 94 PetscErrorCode ierr 95 PetscScalar result 96 PetscReal fnorm 97 98 PetscCall(VecDuplicate(x,f,ierr)) 99 PetscCall(SNESComputeFunction(snes,x,f,ierr)) 100 PetscCall(VecNorm(f,NORM_2,fnorm,ierr)) 101 result = .5*fnorm*fnorm 102 PetscCall(VecDestroy(f,ierr)) 103 end subroutine MyObjective 104 105 program main 106#include <petsc/finclude/petscdmda.h> 107#include <petsc/finclude/petscsnes.h> 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, pointer :: isglobal(:) 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_DMBOUNDARYTYPE,PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMDASTENCILTYPE,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 PetscCallA(PetscObjectReference(solver%isPhi,ierr)) 299 solver%isLambda = isglobal(2) 300 PetscCallA(PetscObjectReference(solver%isLambda,ierr)) 301 302! cache matrices 303 solver%Amat = Amat 304 solver%Bmat = Bmat 305 solver%Cmat = Cmat 306 solver%Dmat = Dmat 307 308 matArray(1) = Amat 309 matArray(2) = Bmat 310 matArray(3) = Cmat 311 matArray(4) = Dmat 312 313 PetscCallA(MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr)) 314 PetscCallA(DMCompositeRestoreGlobalISs(solver%da,isglobal,ierr)) 315 PetscCallA(MatSetFromOptions(KKTmat,ierr)) 316 317! Extract global and local vectors from DMDA; then duplicate for remaining 318! vectors that are the same types 319 PetscCallA(MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr)) 320 PetscCallA(VecDuplicate(x,r,ierr)) 321 322 PetscCallA(SNESCreate(PETSC_COMM_WORLD,mysnes,ierr)) 323 324 PetscCallA(SNESSetDM(mysnes,solver%da,ierr)) 325 326 PetscCallA(SNESSetApplicationContext(mysnes,solver,ierr)) 327 328 PetscCallA(SNESSetDM(mysnes,solver%da,ierr)) 329 330! Set function evaluation routine and vector 331 PetscCallA(SNESSetFunction(mysnes, r, FormFunction, solver,ierr)) 332 if (useobjective .eqv. PETSC_TRUE) then 333 PetscCallA(SNESSetObjective(mysnes, MyObjective, solver, ierr)) 334 endif 335 PetscCallA(SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr)) 336 337! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 338! Customize nonlinear solver; set runtime options 339! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 340! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 341 PetscCallA(SNESSetFromOptions(mysnes,ierr)) 342 343! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 344! Evaluate initial guess; then solve nonlinear system. 345! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 346! Note: The user should initialize the vector, x, with the initial guess 347! for the nonlinear solver prior to calling SNESSolve(). In particular, 348! to employ an initial guess of zero, the user should explicitly set 349! this vector to zero by calling VecSet(). 350 351 PetscCallA(FormInitialGuess(mysnes,x,ierr)) 352 PetscCallA(SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr)) 353 PetscCallA(SNESGetIterationNumber(mysnes,its,ierr)) 354 if (solver%rank .eq. 0) then 355 write(6,100) its 356 endif 357 100 format('Number of SNES iterations = ',i5) 358 359! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 360! Free work space. All PETSc objects should be destroyed when they 361! are no longer needed. 362! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 363 PetscCallA(MatDestroy(KKTmat,ierr)) 364 PetscCallA(MatDestroy(Amat,ierr)) 365 PetscCallA(MatDestroy(Dmat,ierr)) 366 PetscCallA(MatDestroy(Bmat,ierr)) 367 PetscCallA(MatDestroy(Cmat,ierr)) 368 PetscCallA(MatDestroy(solver%AmatLin,ierr)) 369 PetscCallA(ISDestroy(solver%isPhi,ierr)) 370 PetscCallA(ISDestroy(solver%isLambda,ierr)) 371 PetscCallA(VecDestroy(x,ierr)) 372 PetscCallA(VecDestroy(x2,ierr)) 373 PetscCallA(VecDestroy(x1,ierr)) 374 PetscCallA(VecDestroy(x1loc,ierr)) 375 PetscCallA(VecDestroy(x2loc,ierr)) 376 PetscCallA(VecDestroy(r,ierr)) 377 PetscCallA(SNESDestroy(mysnes,ierr)) 378 PetscCallA(DMDestroy(solver%da,ierr)) 379 PetscCallA(DMDestroy(daphi,ierr)) 380 PetscCallA(DMDestroy(dalam,ierr)) 381 382 PetscCallA(PetscFinalize(ierr)) 383 end 384 385! --------------------------------------------------------------------- 386! 387! FormInitialGuess - Forms initial approximation. 388! 389! Input Parameters: 390! X - vector 391! 392! Output Parameter: 393! X - vector 394! 395! Notes: 396! This routine serves as a wrapper for the lower-level routine 397! "InitialGuessLocal", where the actual computations are 398! done using the standard Fortran style of treating the local 399! vector data as a multidimensional array over the local mesh. 400! This routine merely handles ghost point scatters and accesses 401! the local vector data via VecGetArray() and VecRestoreArray(). 402! 403 subroutine FormInitialGuess(mysnes,Xnest,ierr) 404#include <petsc/finclude/petscsnes.h> 405 use petscsnes 406 use ex73f90tmodule 407 use ex73f90tmodule_interfaces 408 implicit none 409! Input/output variables: 410 SNES:: mysnes 411 Vec:: Xnest 412 PetscErrorCode ierr 413 414! Declarations for use with local arrays: 415 type(ex73f90tmodule_type), pointer:: solver 416 Vec:: Xsub(2) 417 PetscInt:: izero,ione,itwo 418 419 izero = 0 420 ione = 1 421 itwo = 2 422 ierr = 0 423 PetscCall(SNESGetApplicationContext(mysnes,solver,ierr)) 424 PetscCall(DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr)) 425 426 PetscCall(InitialGuessLocal(solver,Xsub(1),ierr)) 427 PetscCall(VecAssemblyBegin(Xsub(1),ierr)) 428 PetscCall(VecAssemblyEnd(Xsub(1),ierr)) 429 430! zero out lambda 431 PetscCall(VecZeroEntries(Xsub(2),ierr)) 432 PetscCall(DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr)) 433 434 end subroutine FormInitialGuess 435 436! --------------------------------------------------------------------- 437! 438! InitialGuessLocal - Computes initial approximation, called by 439! the higher level routine FormInitialGuess(). 440! 441! Input Parameter: 442! X1 - local vector data 443! 444! Output Parameters: 445! x - local vector data 446! ierr - error code 447! 448! Notes: 449! This routine uses standard Fortran-style computations over a 2-dim array. 450! 451 subroutine InitialGuessLocal(solver,X1,ierr) 452#include <petsc/finclude/petscsys.h> 453 use petscsys 454 use ex73f90tmodule 455 implicit none 456! Input/output variables: 457 type (ex73f90tmodule_type) solver 458 Vec:: X1 459 PetscErrorCode ierr 460 461! Local variables: 462 PetscInt row,i,j,ione,low,high 463 PetscReal temp1,temp,hx,hy,v 464 PetscReal one 465 466! Set parameters 467 ione = 1 468 ierr = 0 469 one = 1.0 470 hx = one/(solver%mx-1) 471 hy = one/(solver%my-1) 472 temp1 = solver%lambda/(solver%lambda + one) + one 473 474 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 475 476 do 20 row=low,high-1 477 j = row/solver%mx 478 i = mod(row,solver%mx) 479 temp = min(j,solver%my-j+1)*hy 480 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 481 v = 0.0 482 else 483 v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp)) 484 endif 485 PetscCall(VecSetValues(X1,ione,[row],[v],INSERT_VALUES,ierr)) 486 20 continue 487 488 end subroutine InitialGuessLocal 489 490! --------------------------------------------------------------------- 491! 492! FormJacobian - Evaluates Jacobian matrix. 493! 494! Input Parameters: 495! dummy - the SNES context 496! x - input vector 497! solver - solver data 498! 499! Output Parameters: 500! jac - Jacobian matrix 501! jac_prec - optionally different matrix used to construct the preconditioner (not used here) 502! 503 subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr) 504#include <petsc/finclude/petscsnes.h> 505 use petscsnes 506 use ex73f90tmodule 507 implicit none 508! Input/output variables: 509 SNES:: dummy 510 Vec:: X 511 Mat:: jac,jac_prec 512 type(ex73f90tmodule_type) solver 513 PetscErrorCode ierr 514 515! Declarations for use with local arrays: 516 Vec:: Xsub(1) 517 Mat:: Amat 518 PetscInt ione 519 520 ione = 1 521 522 PetscCall(DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr)) 523 524! Compute entries for the locally owned part of the Jacobian preconditioner. 525 PetscCall(MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr)) 526 527 PetscCall(FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr)) 528 PetscCall(MatDestroy(Amat,ierr)) ! discard our reference 529 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr)) 530 531 ! the rest of the matrix is not touched 532 PetscCall(MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 533 PetscCall(MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr)) 534 if (jac .ne. jac_prec) then 535 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 536 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 537 end if 538 539! Tell the matrix we will never add a new nonzero location to the 540! matrix. If we do it will generate an error. 541 PetscCall(MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr)) 542 543 end subroutine FormJacobian 544 545! --------------------------------------------------------------------- 546! 547! FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner, 548! called by the higher level routine FormJacobian(). 549! 550! Input Parameters: 551! x - local vector data 552! 553! Output Parameters: 554! jac - Jacobian matrix used to compute the preconditioner 555! ierr - error code 556! 557! Notes: 558! This routine uses standard Fortran-style computations over a 2-dim array. 559! 560 subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr) 561#include <petsc/finclude/petscmat.h> 562 use ex73f90tmodule 563 implicit none 564! Input/output variables: 565 type (ex73f90tmodule_type) solver 566 Vec:: X1 567 Mat:: jac 568 logical add_nl_term 569 PetscErrorCode ierr 570 571! Local variables: 572 PetscInt irow,row(1),col(5),i,j 573 PetscInt ione,ifive,low,high,ii 574 PetscScalar two,one,hx,hy,hy2inv 575 PetscScalar hx2inv,sc,v(5) 576 PetscScalar,pointer :: lx_v(:) 577 578! Set parameters 579 ione = 1 580 ifive = 5 581 one = 1.0 582 two = 2.0 583 hx = one/(solver%mx-1) 584 hy = one/(solver%my-1) 585 sc = solver%lambda 586 hx2inv = one/(hx*hx) 587 hy2inv = one/(hy*hy) 588 589 PetscCall(VecGetOwnershipRange(X1,low,high,ierr)) 590 PetscCall(VecGetArrayRead(X1,lx_v,ierr)) 591 592 ii = 0 593 do 20 irow=low,high-1 594 j = irow/solver%mx 595 i = mod(irow,solver%mx) 596 ii = ii + 1 ! one based local index 597! boundary points 598 if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 599 col(1) = irow 600 row(1) = irow 601 v(1) = one 602 PetscCall(MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr)) 603! interior grid points 604 else 605 v(1) = -hy2inv 606 if (j-1==0) v(1) = 0.0 607 v(2) = -hx2inv 608 if (i-1==0) v(2) = 0.0 609 v(3) = two*(hx2inv + hy2inv) 610 if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii)) 611 v(4) = -hx2inv 612 if (i+1==solver%mx-1) v(4) = 0.0 613 v(5) = -hy2inv 614 if (j+1==solver%my-1) v(5) = 0.0 615 col(1) = irow - solver%mx 616 col(2) = irow - 1 617 col(3) = irow 618 col(4) = irow + 1 619 col(5) = irow + solver%mx 620 row(1) = irow 621 PetscCall(MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr)) 622 endif 623 20 continue 624 625 PetscCall(VecRestoreArrayRead(X1,lx_v,ierr)) 626 627 end subroutine FormJacobianLocal 628 629! --------------------------------------------------------------------- 630! 631! FormFunction - Evaluates nonlinear function, F(x). 632! 633! Input Parameters: 634! snes - the SNES context 635! X - input vector 636! dummy - optional user-defined context, as set by SNESSetFunction() 637! (not used here) 638! 639! Output Parameter: 640! F - function vector 641! 642 subroutine FormFunction(snesIn,X,F,solver,ierr) 643#include <petsc/finclude/petscsnes.h> 644 use petscsnes 645 use ex73f90tmodule 646 implicit none 647! Input/output variables: 648 SNES:: snesIn 649 Vec:: X,F 650 PetscErrorCode ierr 651 type (ex73f90tmodule_type) solver 652 653! Declarations for use with local arrays: 654 Vec:: Xsub(2),Fsub(2) 655 PetscInt itwo 656 657! Scatter ghost points to local vector, using the 2-step process 658! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 659! By placing code between these two statements, computations can 660! be done while messages are in transition. 661 662 itwo = 2 663 PetscCall(DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr)) 664 PetscCall(DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER_ARRAY,Fsub,ierr)) 665 666 PetscCall(FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr)) 667 PetscCall(MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr)) 668 669! do rest of operator (linear) 670 PetscCall(MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr)) 671 PetscCall(MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr)) 672 PetscCall(MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr)) 673 674 PetscCall(DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER_ARRAY,Xsub,ierr)) 675 PetscCall(DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER_ARRAY,Fsub,ierr)) 676 end subroutine formfunction 677 678! --------------------------------------------------------------------- 679! 680! FormFunctionNLTerm - Computes nonlinear function, called by 681! the higher level routine FormFunction(). 682! 683! Input Parameter: 684! x - local vector data 685! 686! Output Parameters: 687! f - local vector data, f(x) 688! ierr - error code 689! 690! Notes: 691! This routine uses standard Fortran-style computations over a 2-dim array. 692! 693 subroutine FormFunctionNLTerm(X1,F1,solver,ierr) 694#include <petsc/finclude/petscvec.h> 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(VecGetArrayRead(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(VecRestoreArrayRead(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 {{secant 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