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