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