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