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