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