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