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