1c4762a1bSJed Brown! 2c4762a1bSJed Brown! Description: Solves a nonlinear system in parallel with SNES. 3c4762a1bSJed Brown! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular 4c4762a1bSJed Brown! domain, using distributed arrays (DMDAs) to partition the parallel grid. 5c4762a1bSJed Brown! The command line options include: 6c4762a1bSJed Brown! -par <parameter>, where <parameter> indicates the nonlinearity of the problem 7c4762a1bSJed Brown! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 8c4762a1bSJed Brown! 9c4762a1bSJed Brown 10c4762a1bSJed Brown! 11c4762a1bSJed Brown! -------------------------------------------------------------------------- 12c4762a1bSJed Brown! 13c4762a1bSJed Brown! Solid Fuel Ignition (SFI) problem. This problem is modeled by 14c4762a1bSJed Brown! the partial differential equation 15c4762a1bSJed Brown! 16c4762a1bSJed Brown! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 17c4762a1bSJed Brown! 18c4762a1bSJed Brown! with boundary conditions 19c4762a1bSJed Brown! 20c4762a1bSJed Brown! u = 0 for x = 0, x = 1, y = 0, y = 1. 21c4762a1bSJed Brown! 22c4762a1bSJed Brown! A finite difference approximation with the usual 5-point stencil 23c4762a1bSJed Brown! is used to discretize the boundary value problem to obtain a nonlinear 24c4762a1bSJed Brown! system of equations. 25c4762a1bSJed Brown! 26c4762a1bSJed Brown! The uniprocessor version of this code is snes/tutorials/ex4f.F 27c4762a1bSJed Brown! 28c4762a1bSJed Brown! -------------------------------------------------------------------------- 29c4762a1bSJed Brown! The following define must be used before including any PETSc include files 30c4762a1bSJed Brown! into a module or interface. This is because they can't handle declarations 31c4762a1bSJed Brown! in them 32c4762a1bSJed Brown! 33ce78bad3SBarry Smith#include <petsc/finclude/petscsnes.h> 34ce78bad3SBarry Smith#include <petsc/finclude/petscdmda.h> 35c5e229c2SMartin Diehlmodule ex5f90module 36c4762a1bSJed Brown use petscsnes 37dfbbaf82SBarry Smith use petscdmda 38*e7a95102SMartin Diehl implicit none 39c4762a1bSJed Brown type userctx 40c4762a1bSJed Brown PetscInt xs, xe, xm, gxs, gxe, gxm 41c4762a1bSJed Brown PetscInt ys, ye, ym, gys, gye, gym 42c4762a1bSJed Brown PetscInt mx, my 43c4762a1bSJed Brown PetscMPIInt rank 44c4762a1bSJed Brown PetscReal lambda 45c4762a1bSJed Brown end type userctx 46c4762a1bSJed Brown 47*e7a95102SMartin Diehl interface 48*e7a95102SMartin Diehl subroutine SNESSetApplicationContext(snes, ctx, ierr) 49*e7a95102SMartin Diehl use petscsnes 50*e7a95102SMartin Diehl import userctx 51*e7a95102SMartin Diehl implicit none 52*e7a95102SMartin Diehl SNES snes 53*e7a95102SMartin Diehl type(userctx) ctx 54*e7a95102SMartin Diehl PetscErrorCode ierr 55*e7a95102SMartin Diehl end subroutine 56*e7a95102SMartin Diehl subroutine SNESGetApplicationContext(snes, ctx, ierr) 57*e7a95102SMartin Diehl use petscsnes 58*e7a95102SMartin Diehl import userctx 59*e7a95102SMartin Diehl implicit none 60*e7a95102SMartin Diehl SNES snes 61*e7a95102SMartin Diehl type(userctx), pointer :: ctx 62*e7a95102SMartin Diehl PetscErrorCode ierr 63*e7a95102SMartin Diehl end subroutine 64*e7a95102SMartin Diehl end interface 65*e7a95102SMartin Diehl 66c4762a1bSJed Browncontains 67c4762a1bSJed Brown! --------------------------------------------------------------------- 68c4762a1bSJed Brown! 69c4762a1bSJed Brown! FormFunction - Evaluates nonlinear function, F(x). 70c4762a1bSJed Brown! 71c4762a1bSJed Brown! Input Parameters: 72c4762a1bSJed Brown! snes - the SNES context 73c4762a1bSJed Brown! X - input vector 74c4762a1bSJed Brown! dummy - optional user-defined context, as set by SNESSetFunction() 75c4762a1bSJed Brown! (not used here) 76c4762a1bSJed Brown! 77c4762a1bSJed Brown! Output Parameter: 78c4762a1bSJed Brown! F - function vector 79c4762a1bSJed Brown! 80c4762a1bSJed Brown! Notes: 81c4762a1bSJed Brown! This routine serves as a wrapper for the lower-level routine 82c4762a1bSJed Brown! "FormFunctionLocal", where the actual computations are 83c4762a1bSJed Brown! done using the standard Fortran style of treating the local 84c4762a1bSJed Brown! vector data as a multidimensional array over the local mesh. 85c4762a1bSJed Brown! This routine merely handles ghost point scatters and accesses 86ce78bad3SBarry Smith! the local vector data via VecGetArray() and VecRestoreArray(). 87c4762a1bSJed Brown! 88c4762a1bSJed Brown subroutine FormFunction(snes, X, F, user, ierr) 89c4762a1bSJed Brown implicit none 90c4762a1bSJed Brown 91c4762a1bSJed Brown! Input/output variables: 92c4762a1bSJed Brown SNES snes 93c4762a1bSJed Brown Vec X, F 94c4762a1bSJed Brown PetscErrorCode ierr 95c4762a1bSJed Brown type(userctx) user 96c4762a1bSJed Brown DM da 97c4762a1bSJed Brown 98c4762a1bSJed Brown! Declarations for use with local arrays: 99c4762a1bSJed Brown PetscScalar, pointer :: lx_v(:), lf_v(:) 100c4762a1bSJed Brown Vec localX 101c4762a1bSJed Brown 102c4762a1bSJed Brown! Scatter ghost points to local vector, using the 2-step process 103c4762a1bSJed Brown! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 104c4762a1bSJed Brown! By placing code between these two statements, computations can 105c4762a1bSJed Brown! be done while messages are in transition. 106d8606c27SBarry Smith PetscCall(SNESGetDM(snes, da, ierr)) 107d8606c27SBarry Smith PetscCall(DMGetLocalVector(da, localX, ierr)) 108d8606c27SBarry Smith PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr)) 109d8606c27SBarry Smith PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr)) 110c4762a1bSJed Brown 111c4762a1bSJed Brown! Get a pointer to vector data. 112ce78bad3SBarry Smith! - For default PETSc vectors, VecGetArray() returns a pointer to 113c4762a1bSJed Brown! the data array. Otherwise, the routine is implementation dependent. 114ce78bad3SBarry Smith! - You MUST call VecRestoreArray() when you no longer need access to 115c4762a1bSJed Brown! the array. 116ce78bad3SBarry Smith! - Note that the interface to VecGetArray() differs from VecGetArray(). 117c4762a1bSJed Brown 118ce78bad3SBarry Smith PetscCall(VecGetArray(localX, lx_v, ierr)) 119ce78bad3SBarry Smith PetscCall(VecGetArray(F, lf_v, ierr)) 120c4762a1bSJed Brown 121c4762a1bSJed Brown! Compute function over the locally owned part of the grid 122d8606c27SBarry Smith PetscCall(FormFunctionLocal(lx_v, lf_v, user, ierr)) 123c4762a1bSJed Brown 124c4762a1bSJed Brown! Restore vectors 125ce78bad3SBarry Smith PetscCall(VecRestoreArray(localX, lx_v, ierr)) 126ce78bad3SBarry Smith PetscCall(VecRestoreArray(F, lf_v, ierr)) 127c4762a1bSJed Brown 128c4762a1bSJed Brown! Insert values into global vector 129c4762a1bSJed Brown 130d8606c27SBarry Smith PetscCall(DMRestoreLocalVector(da, localX, ierr)) 131d8606c27SBarry Smith PetscCall(PetscLogFlops(11.0d0*user%ym*user%xm, ierr)) 132c4762a1bSJed Brown 133d8606c27SBarry Smith! PetscCallA(VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)) 134d8606c27SBarry Smith! PetscCallA(VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr)) 135c4762a1bSJed Brown end subroutine formfunction 136*e7a95102SMartin Diehl 137*e7a95102SMartin Diehl! --------------------------------------------------------------------- 138*e7a95102SMartin Diehl! 139*e7a95102SMartin Diehl! FormInitialGuess - Forms initial approximation. 140*e7a95102SMartin Diehl! 141*e7a95102SMartin Diehl! Input Parameters: 142*e7a95102SMartin Diehl! X - vector 143*e7a95102SMartin Diehl! 144*e7a95102SMartin Diehl! Output Parameter: 145*e7a95102SMartin Diehl! X - vector 146*e7a95102SMartin Diehl! 147*e7a95102SMartin Diehl! Notes: 148*e7a95102SMartin Diehl! This routine serves as a wrapper for the lower-level routine 149*e7a95102SMartin Diehl! "InitialGuessLocal", where the actual computations are 150*e7a95102SMartin Diehl! done using the standard Fortran style of treating the local 151*e7a95102SMartin Diehl! vector data as a multidimensional array over the local mesh. 152*e7a95102SMartin Diehl! This routine merely handles ghost point scatters and accesses 153*e7a95102SMartin Diehl! the local vector data via VecGetArray() and VecRestoreArray(). 154*e7a95102SMartin Diehl! 155*e7a95102SMartin Diehl subroutine FormInitialGuess(snes, X, ierr) 156*e7a95102SMartin Diehl! Input/output variables: 157*e7a95102SMartin Diehl SNES snes 158*e7a95102SMartin Diehl type(userctx), pointer:: puser 159*e7a95102SMartin Diehl Vec X 160*e7a95102SMartin Diehl PetscErrorCode ierr 161*e7a95102SMartin Diehl DM da 162*e7a95102SMartin Diehl 163*e7a95102SMartin Diehl! Declarations for use with local arrays: 164*e7a95102SMartin Diehl PetscScalar, pointer :: lx_v(:) 165*e7a95102SMartin Diehl 166*e7a95102SMartin Diehl ierr = 0 167*e7a95102SMartin Diehl PetscCallA(SNESGetDM(snes, da, ierr)) 168*e7a95102SMartin Diehl PetscCallA(SNESGetApplicationContext(snes, puser, ierr)) 169*e7a95102SMartin Diehl! Get a pointer to vector data. 170*e7a95102SMartin Diehl! - For default PETSc vectors, VecGetArray() returns a pointer to 171*e7a95102SMartin Diehl! the data array. Otherwise, the routine is implementation dependent. 172*e7a95102SMartin Diehl! - You MUST call VecRestoreArray() when you no longer need access to 173*e7a95102SMartin Diehl! the array. 174*e7a95102SMartin Diehl! - Note that the interface to VecGetArray() differs from VecGetArray(). 175*e7a95102SMartin Diehl 176*e7a95102SMartin Diehl PetscCallA(VecGetArray(X, lx_v, ierr)) 177*e7a95102SMartin Diehl 178*e7a95102SMartin Diehl! Compute initial guess over the locally owned part of the grid 179*e7a95102SMartin Diehl PetscCallA(InitialGuessLocal(puser, lx_v, ierr)) 180*e7a95102SMartin Diehl 181*e7a95102SMartin Diehl! Restore vector 182*e7a95102SMartin Diehl PetscCallA(VecRestoreArray(X, lx_v, ierr)) 183*e7a95102SMartin Diehl 184*e7a95102SMartin Diehl! Insert values into global vector 185*e7a95102SMartin Diehl 186*e7a95102SMartin Diehl end 187*e7a95102SMartin Diehl 188*e7a95102SMartin Diehl! --------------------------------------------------------------------- 189*e7a95102SMartin Diehl! 190*e7a95102SMartin Diehl! InitialGuessLocal - Computes initial approximation, called by 191*e7a95102SMartin Diehl! the higher level routine FormInitialGuess(). 192*e7a95102SMartin Diehl! 193*e7a95102SMartin Diehl! Input Parameter: 194*e7a95102SMartin Diehl! x - local vector data 195*e7a95102SMartin Diehl! 196*e7a95102SMartin Diehl! Output Parameters: 197*e7a95102SMartin Diehl! x - local vector data 198*e7a95102SMartin Diehl! ierr - error code 199*e7a95102SMartin Diehl! 200*e7a95102SMartin Diehl! Notes: 201*e7a95102SMartin Diehl! This routine uses standard Fortran-style computations over a 2-dim array. 202*e7a95102SMartin Diehl! 203*e7a95102SMartin Diehl subroutine InitialGuessLocal(user, x, ierr) 204*e7a95102SMartin Diehl! Input/output variables: 205*e7a95102SMartin Diehl type(userctx) user 206*e7a95102SMartin Diehl PetscScalar x(user%xs:user%xe, user%ys:user%ye) 207*e7a95102SMartin Diehl PetscErrorCode ierr 208*e7a95102SMartin Diehl 209*e7a95102SMartin Diehl! Local variables: 210*e7a95102SMartin Diehl PetscInt i, j 211*e7a95102SMartin Diehl PetscReal temp1, temp, hx, hy 212*e7a95102SMartin Diehl PetscReal one 213*e7a95102SMartin Diehl 214*e7a95102SMartin Diehl! Set parameters 215*e7a95102SMartin Diehl 216*e7a95102SMartin Diehl ierr = 0 217*e7a95102SMartin Diehl one = 1.0 218*e7a95102SMartin Diehl hx = one/(user%mx - 1) 219*e7a95102SMartin Diehl hy = one/(user%my - 1) 220*e7a95102SMartin Diehl temp1 = user%lambda/(user%lambda + one) 221*e7a95102SMartin Diehl 222*e7a95102SMartin Diehl do j = user%ys, user%ye 223*e7a95102SMartin Diehl temp = min(j - 1, user%my - j)*hy 224*e7a95102SMartin Diehl do i = user%xs, user%xe 225*e7a95102SMartin Diehl if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 226*e7a95102SMartin Diehl x(i, j) = 0.0 227*e7a95102SMartin Diehl else 228*e7a95102SMartin Diehl x(i, j) = temp1*sqrt(min(hx*min(i - 1, user%mx - i), temp)) 229*e7a95102SMartin Diehl end if 230*e7a95102SMartin Diehl end do 231*e7a95102SMartin Diehl end do 232*e7a95102SMartin Diehl 233*e7a95102SMartin Diehl end 234*e7a95102SMartin Diehl 235*e7a95102SMartin Diehl! --------------------------------------------------------------------- 236*e7a95102SMartin Diehl! 237*e7a95102SMartin Diehl! FormFunctionLocal - Computes nonlinear function, called by 238*e7a95102SMartin Diehl! the higher level routine FormFunction(). 239*e7a95102SMartin Diehl! 240*e7a95102SMartin Diehl! Input Parameter: 241*e7a95102SMartin Diehl! x - local vector data 242*e7a95102SMartin Diehl! 243*e7a95102SMartin Diehl! Output Parameters: 244*e7a95102SMartin Diehl! f - local vector data, f(x) 245*e7a95102SMartin Diehl! ierr - error code 246*e7a95102SMartin Diehl! 247*e7a95102SMartin Diehl! Notes: 248*e7a95102SMartin Diehl! This routine uses standard Fortran-style computations over a 2-dim array. 249*e7a95102SMartin Diehl! 250*e7a95102SMartin Diehl subroutine FormFunctionLocal(x, f, user, ierr) 251*e7a95102SMartin Diehl! Input/output variables: 252*e7a95102SMartin Diehl type(userctx) user 253*e7a95102SMartin Diehl PetscScalar x(user%gxs:user%gxe, user%gys:user%gye) 254*e7a95102SMartin Diehl PetscScalar f(user%xs:user%xe, user%ys:user%ye) 255*e7a95102SMartin Diehl PetscErrorCode ierr 256*e7a95102SMartin Diehl 257*e7a95102SMartin Diehl! Local variables: 258*e7a95102SMartin Diehl PetscScalar two, one, hx, hy, hxdhy, hydhx, sc 259*e7a95102SMartin Diehl PetscScalar u, uxx, uyy 260*e7a95102SMartin Diehl PetscInt i, j 261*e7a95102SMartin Diehl 262*e7a95102SMartin Diehl one = 1.0 263*e7a95102SMartin Diehl two = 2.0 264*e7a95102SMartin Diehl hx = one/(user%mx - 1) 265*e7a95102SMartin Diehl hy = one/(user%my - 1) 266*e7a95102SMartin Diehl sc = hx*hy*user%lambda 267*e7a95102SMartin Diehl hxdhy = hx/hy 268*e7a95102SMartin Diehl hydhx = hy/hx 269*e7a95102SMartin Diehl 270*e7a95102SMartin Diehl! Compute function over the locally owned part of the grid 271*e7a95102SMartin Diehl 272*e7a95102SMartin Diehl do j = user%ys, user%ye 273*e7a95102SMartin Diehl do i = user%xs, user%xe 274*e7a95102SMartin Diehl if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 275*e7a95102SMartin Diehl f(i, j) = x(i, j) 276*e7a95102SMartin Diehl else 277*e7a95102SMartin Diehl u = x(i, j) 278*e7a95102SMartin Diehl uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j)) 279*e7a95102SMartin Diehl uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1)) 280*e7a95102SMartin Diehl f(i, j) = uxx + uyy - sc*exp(u) 281*e7a95102SMartin Diehl end if 282*e7a95102SMartin Diehl end do 283*e7a95102SMartin Diehl end do 284*e7a95102SMartin Diehl 285*e7a95102SMartin Diehl end 286*e7a95102SMartin Diehl 287*e7a95102SMartin Diehl! --------------------------------------------------------------------- 288*e7a95102SMartin Diehl! 289*e7a95102SMartin Diehl! FormJacobian - Evaluates Jacobian matrix. 290*e7a95102SMartin Diehl! 291*e7a95102SMartin Diehl! Input Parameters: 292*e7a95102SMartin Diehl! snes - the SNES context 293*e7a95102SMartin Diehl! x - input vector 294*e7a95102SMartin Diehl! dummy - optional user-defined context, as set by SNESSetJacobian() 295*e7a95102SMartin Diehl! (not used here) 296*e7a95102SMartin Diehl! 297*e7a95102SMartin Diehl! Output Parameters: 298*e7a95102SMartin Diehl! jac - Jacobian matrix 299*e7a95102SMartin Diehl! jac_prec - optionally different matrix used to construct the preconditioner (not used here) 300*e7a95102SMartin Diehl! 301*e7a95102SMartin Diehl! Notes: 302*e7a95102SMartin Diehl! This routine serves as a wrapper for the lower-level routine 303*e7a95102SMartin Diehl! "FormJacobianLocal", where the actual computations are 304*e7a95102SMartin Diehl! done using the standard Fortran style of treating the local 305*e7a95102SMartin Diehl! vector data as a multidimensional array over the local mesh. 306*e7a95102SMartin Diehl! This routine merely accesses the local vector data via 307*e7a95102SMartin Diehl! VecGetArray() and VecRestoreArray(). 308*e7a95102SMartin Diehl! 309*e7a95102SMartin Diehl! Notes: 310*e7a95102SMartin Diehl! Due to grid point reordering with DMDAs, we must always work 311*e7a95102SMartin Diehl! with the local grid points, and then transform them to the new 312*e7a95102SMartin Diehl! global numbering with the "ltog" mapping 313*e7a95102SMartin Diehl! We cannot work directly with the global numbers for the original 314*e7a95102SMartin Diehl! uniprocessor grid! 315*e7a95102SMartin Diehl! 316*e7a95102SMartin Diehl! Two methods are available for imposing this transformation 317*e7a95102SMartin Diehl! when setting matrix entries: 318*e7a95102SMartin Diehl! (A) MatSetValuesLocal(), using the local ordering (including 319*e7a95102SMartin Diehl! ghost points!) 320*e7a95102SMartin Diehl! - Set matrix entries using the local ordering 321*e7a95102SMartin Diehl! by calling MatSetValuesLocal() 322*e7a95102SMartin Diehl! (B) MatSetValues(), using the global ordering 323*e7a95102SMartin Diehl 324*e7a95102SMartin Diehl! - Set matrix entries using the global ordering by calling 325*e7a95102SMartin Diehl! MatSetValues() 326*e7a95102SMartin Diehl! Option (A) seems cleaner/easier in many cases, and is the procedure 327*e7a95102SMartin Diehl! used in this example. 328*e7a95102SMartin Diehl! 329*e7a95102SMartin Diehl subroutine FormJacobian(snes, X, jac, jac_prec, user, ierr) 330*e7a95102SMartin Diehl! Input/output variables: 331*e7a95102SMartin Diehl SNES snes 332*e7a95102SMartin Diehl Vec X 333*e7a95102SMartin Diehl Mat jac, jac_prec 334*e7a95102SMartin Diehl type(userctx) user 335*e7a95102SMartin Diehl PetscErrorCode ierr 336*e7a95102SMartin Diehl DM da 337*e7a95102SMartin Diehl 338*e7a95102SMartin Diehl! Declarations for use with local arrays: 339*e7a95102SMartin Diehl PetscScalar, pointer :: lx_v(:) 340*e7a95102SMartin Diehl Vec localX 341*e7a95102SMartin Diehl 342*e7a95102SMartin Diehl! Scatter ghost points to local vector, using the 2-step process 343*e7a95102SMartin Diehl! DMGlobalToLocalBegin(), DMGlobalToLocalEnd() 344*e7a95102SMartin Diehl! Computations can be done while messages are in transition, 345*e7a95102SMartin Diehl! by placing code between these two statements. 346*e7a95102SMartin Diehl 347*e7a95102SMartin Diehl PetscCallA(SNESGetDM(snes, da, ierr)) 348*e7a95102SMartin Diehl PetscCallA(DMGetLocalVector(da, localX, ierr)) 349*e7a95102SMartin Diehl PetscCallA(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr)) 350*e7a95102SMartin Diehl PetscCallA(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr)) 351*e7a95102SMartin Diehl 352*e7a95102SMartin Diehl! Get a pointer to vector data 353*e7a95102SMartin Diehl PetscCallA(VecGetArray(localX, lx_v, ierr)) 354*e7a95102SMartin Diehl 355*e7a95102SMartin Diehl! Compute entries for the locally owned part of the Jacobian preconditioner. 356*e7a95102SMartin Diehl PetscCallA(FormJacobianLocal(lx_v, jac_prec, user, ierr)) 357*e7a95102SMartin Diehl 358*e7a95102SMartin Diehl! Assemble matrix, using the 2-step process: 359*e7a95102SMartin Diehl! MatAssemblyBegin(), MatAssemblyEnd() 360*e7a95102SMartin Diehl! Computations can be done while messages are in transition, 361*e7a95102SMartin Diehl! by placing code between these two statements. 362*e7a95102SMartin Diehl 363*e7a95102SMartin Diehl PetscCallA(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)) 364*e7a95102SMartin Diehl if (jac /= jac_prec) then 365*e7a95102SMartin Diehl PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 366*e7a95102SMartin Diehl end if 367*e7a95102SMartin Diehl PetscCallA(VecRestoreArray(localX, lx_v, ierr)) 368*e7a95102SMartin Diehl PetscCallA(DMRestoreLocalVector(da, localX, ierr)) 369*e7a95102SMartin Diehl PetscCallA(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)) 370*e7a95102SMartin Diehl if (jac /= jac_prec) then 371*e7a95102SMartin Diehl PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 372*e7a95102SMartin Diehl end if 373*e7a95102SMartin Diehl 374*e7a95102SMartin Diehl! Tell the matrix we will never add a new nonzero location to the 375*e7a95102SMartin Diehl! matrix. If we do it will generate an error. 376*e7a95102SMartin Diehl 377*e7a95102SMartin Diehl PetscCallA(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr)) 378*e7a95102SMartin Diehl 379*e7a95102SMartin Diehl end 380*e7a95102SMartin Diehl 381*e7a95102SMartin Diehl! --------------------------------------------------------------------- 382*e7a95102SMartin Diehl! 383*e7a95102SMartin Diehl! FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner, 384*e7a95102SMartin Diehl! called by the higher level routine FormJacobian(). 385*e7a95102SMartin Diehl! 386*e7a95102SMartin Diehl! Input Parameters: 387*e7a95102SMartin Diehl! x - local vector data 388*e7a95102SMartin Diehl! 389*e7a95102SMartin Diehl! Output Parameters: 390*e7a95102SMartin Diehl! jac_prec - Jacobian matrix used to compute the preconditioner 391*e7a95102SMartin Diehl! ierr - error code 392*e7a95102SMartin Diehl! 393*e7a95102SMartin Diehl! Notes: 394*e7a95102SMartin Diehl! This routine uses standard Fortran-style computations over a 2-dim array. 395*e7a95102SMartin Diehl! 396*e7a95102SMartin Diehl! Notes: 397*e7a95102SMartin Diehl! Due to grid point reordering with DMDAs, we must always work 398*e7a95102SMartin Diehl! with the local grid points, and then transform them to the new 399*e7a95102SMartin Diehl! global numbering with the "ltog" mapping 400*e7a95102SMartin Diehl! We cannot work directly with the global numbers for the original 401*e7a95102SMartin Diehl! uniprocessor grid! 402*e7a95102SMartin Diehl! 403*e7a95102SMartin Diehl! Two methods are available for imposing this transformation 404*e7a95102SMartin Diehl! when setting matrix entries: 405*e7a95102SMartin Diehl! (A) MatSetValuesLocal(), using the local ordering (including 406*e7a95102SMartin Diehl! ghost points!) 407*e7a95102SMartin Diehl! - Set matrix entries using the local ordering 408*e7a95102SMartin Diehl! by calling MatSetValuesLocal() 409*e7a95102SMartin Diehl! (B) MatSetValues(), using the global ordering 410*e7a95102SMartin Diehl! - Then apply this map explicitly yourself 411*e7a95102SMartin Diehl! - Set matrix entries using the global ordering by calling 412*e7a95102SMartin Diehl! MatSetValues() 413*e7a95102SMartin Diehl! Option (A) seems cleaner/easier in many cases, and is the procedure 414*e7a95102SMartin Diehl! used in this example. 415*e7a95102SMartin Diehl! 416*e7a95102SMartin Diehl subroutine FormJacobianLocal(x, jac_prec, user, ierr) 417*e7a95102SMartin Diehl! Input/output variables: 418*e7a95102SMartin Diehl type(userctx) user 419*e7a95102SMartin Diehl PetscScalar x(user%gxs:user%gxe, user%gys:user%gye) 420*e7a95102SMartin Diehl Mat jac_prec 421*e7a95102SMartin Diehl PetscErrorCode ierr 422*e7a95102SMartin Diehl 423*e7a95102SMartin Diehl! Local variables: 424*e7a95102SMartin Diehl PetscInt row, col(5), i, j 425*e7a95102SMartin Diehl PetscInt ione, ifive 426*e7a95102SMartin Diehl PetscScalar two, one, hx, hy, hxdhy 427*e7a95102SMartin Diehl PetscScalar hydhx, sc, v(5) 428*e7a95102SMartin Diehl 429*e7a95102SMartin Diehl! Set parameters 430*e7a95102SMartin Diehl ione = 1 431*e7a95102SMartin Diehl ifive = 5 432*e7a95102SMartin Diehl one = 1.0 433*e7a95102SMartin Diehl two = 2.0 434*e7a95102SMartin Diehl hx = one/(user%mx - 1) 435*e7a95102SMartin Diehl hy = one/(user%my - 1) 436*e7a95102SMartin Diehl sc = hx*hy 437*e7a95102SMartin Diehl hxdhy = hx/hy 438*e7a95102SMartin Diehl hydhx = hy/hx 439*e7a95102SMartin Diehl 440*e7a95102SMartin Diehl! Compute entries for the locally owned part of the Jacobian. 441*e7a95102SMartin Diehl! - Currently, all PETSc parallel matrix formats are partitioned by 442*e7a95102SMartin Diehl! contiguous chunks of rows across the processors. 443*e7a95102SMartin Diehl! - Each processor needs to insert only elements that it owns 444*e7a95102SMartin Diehl! locally (but any non-local elements will be sent to the 445*e7a95102SMartin Diehl! appropriate processor during matrix assembly). 446*e7a95102SMartin Diehl! - Here, we set all entries for a particular row at once. 447*e7a95102SMartin Diehl! - We can set matrix entries either using either 448*e7a95102SMartin Diehl! MatSetValuesLocal() or MatSetValues(), as discussed above. 449*e7a95102SMartin Diehl! - Note that MatSetValues() uses 0-based row and column numbers 450*e7a95102SMartin Diehl! in Fortran as well as in C. 451*e7a95102SMartin Diehl 452*e7a95102SMartin Diehl do j = user%ys, user%ye 453*e7a95102SMartin Diehl row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1 454*e7a95102SMartin Diehl do i = user%xs, user%xe 455*e7a95102SMartin Diehl row = row + 1 456*e7a95102SMartin Diehl! boundary points 457*e7a95102SMartin Diehl if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 458*e7a95102SMartin Diehl col(1) = row 459*e7a95102SMartin Diehl v(1) = one 460*e7a95102SMartin Diehl PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ione, col, v, INSERT_VALUES, ierr)) 461*e7a95102SMartin Diehl! interior grid points 462*e7a95102SMartin Diehl else 463*e7a95102SMartin Diehl v(1) = -hxdhy 464*e7a95102SMartin Diehl v(2) = -hydhx 465*e7a95102SMartin Diehl v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i, j)) 466*e7a95102SMartin Diehl v(4) = -hydhx 467*e7a95102SMartin Diehl v(5) = -hxdhy 468*e7a95102SMartin Diehl col(1) = row - user%gxm 469*e7a95102SMartin Diehl col(2) = row - 1 470*e7a95102SMartin Diehl col(3) = row 471*e7a95102SMartin Diehl col(4) = row + 1 472*e7a95102SMartin Diehl col(5) = row + user%gxm 473*e7a95102SMartin Diehl PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ifive, col, v, INSERT_VALUES, ierr)) 474*e7a95102SMartin Diehl end if 475*e7a95102SMartin Diehl end do 476*e7a95102SMartin Diehl end do 477*e7a95102SMartin Diehl 478*e7a95102SMartin Diehl end 479*e7a95102SMartin Diehl 480dfbbaf82SBarry Smithend module ex5f90module 481c4762a1bSJed Brown 482c4762a1bSJed Brownprogram main 483dfbbaf82SBarry Smith use ex5f90module 484c4762a1bSJed Brown implicit none 485c4762a1bSJed Brown! 486c4762a1bSJed Brown 487c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 488c4762a1bSJed Brown! Variable declarations 489c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 490c4762a1bSJed Brown! 491c4762a1bSJed Brown! Variables: 492c4762a1bSJed Brown! snes - nonlinear solver 493c4762a1bSJed Brown! x, r - solution, residual vectors 494c4762a1bSJed Brown! J - Jacobian matrix 495c4762a1bSJed Brown! its - iterations for convergence 496c4762a1bSJed Brown! Nx, Ny - number of preocessors in x- and y- directions 497c4762a1bSJed Brown! matrix_free - flag - 1 indicates matrix-free version 498c4762a1bSJed Brown! 499c4762a1bSJed Brown SNES snes 500c4762a1bSJed Brown Vec x, r 501c4762a1bSJed Brown Mat J 502c4762a1bSJed Brown PetscErrorCode ierr 503c4762a1bSJed Brown PetscInt its 504c4762a1bSJed Brown PetscBool flg, matrix_free 505c4762a1bSJed Brown PetscInt ione, nfour 506c4762a1bSJed Brown PetscReal lambda_max, lambda_min 507c4762a1bSJed Brown type(userctx) user 508c4762a1bSJed Brown DM da 509c4762a1bSJed Brown 510c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 511c4762a1bSJed Brown! Initialize program 512c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 513d8606c27SBarry Smith PetscCallA(PetscInitialize(ierr)) 514d8606c27SBarry Smith PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, user%rank, ierr)) 515c4762a1bSJed Brown 516c4762a1bSJed Brown! Initialize problem parameters 517c4762a1bSJed Brown lambda_max = 6.81 518c4762a1bSJed Brown lambda_min = 0.0 519c4762a1bSJed Brown user%lambda = 6.0 520c4762a1bSJed Brown ione = 1 521c4762a1bSJed Brown nfour = 4 522d8606c27SBarry Smith PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', user%lambda, flg, ierr)) 5234820e4eaSBarry Smith PetscCheckA(user%lambda < lambda_max .and. user%lambda > lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range') 524c4762a1bSJed Brown 525c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 526c4762a1bSJed Brown! Create nonlinear solver context 527c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 528d8606c27SBarry Smith PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr)) 529c4762a1bSJed Brown 530c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 531c4762a1bSJed Brown! Create vector data structures; set function evaluation routine 532c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 533c4762a1bSJed Brown 534c4762a1bSJed Brown! Create distributed array (DMDA) to manage parallel grid and vectors 535c4762a1bSJed Brown 536c4762a1bSJed Brown! This really needs only the star-type stencil, but we use the box 537c4762a1bSJed Brown! stencil temporarily. 5385d83a8b1SBarry Smith 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, da, ierr)) 539d8606c27SBarry Smith PetscCallA(DMSetFromOptions(da, ierr)) 540d8606c27SBarry Smith PetscCallA(DMSetUp(da, ierr)) 541c4762a1bSJed Brown 542ce78bad3SBarry Smith PetscCallA(DMDAGetInfo(da, PETSC_NULL_INTEGER, user%mx, user%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)) 543c4762a1bSJed Brown 544c4762a1bSJed Brown! 545c4762a1bSJed Brown! Visualize the distribution of the array across the processors 546c4762a1bSJed Brown! 547d8606c27SBarry Smith! PetscCallA(DMView(da,PETSC_VIEWER_DRAW_WORLD,ierr)) 548c4762a1bSJed Brown 549c4762a1bSJed Brown! Extract global and local vectors from DMDA; then duplicate for remaining 550c4762a1bSJed Brown! vectors that are the same types 551d8606c27SBarry Smith PetscCallA(DMCreateGlobalVector(da, x, ierr)) 552d8606c27SBarry Smith PetscCallA(VecDuplicate(x, r, ierr)) 553c4762a1bSJed Brown 554c4762a1bSJed Brown! Get local grid boundaries (for 2-dimensional DMDA) 555d8606c27SBarry Smith PetscCallA(DMDAGetCorners(da, user%xs, user%ys, PETSC_NULL_INTEGER, user%xm, user%ym, PETSC_NULL_INTEGER, ierr)) 556d8606c27SBarry Smith PetscCallA(DMDAGetGhostCorners(da, user%gxs, user%gys, PETSC_NULL_INTEGER, user%gxm, user%gym, PETSC_NULL_INTEGER, ierr)) 557c4762a1bSJed Brown 558c4762a1bSJed Brown! Here we shift the starting indices up by one so that we can easily 559c4762a1bSJed Brown! use the Fortran convention of 1-based indices (rather 0-based indices). 560c4762a1bSJed Brown user%xs = user%xs + 1 561c4762a1bSJed Brown user%ys = user%ys + 1 562c4762a1bSJed Brown user%gxs = user%gxs + 1 563c4762a1bSJed Brown user%gys = user%gys + 1 564c4762a1bSJed Brown 565c4762a1bSJed Brown user%ye = user%ys + user%ym - 1 566c4762a1bSJed Brown user%xe = user%xs + user%xm - 1 567c4762a1bSJed Brown user%gye = user%gys + user%gym - 1 568c4762a1bSJed Brown user%gxe = user%gxs + user%gxm - 1 569c4762a1bSJed Brown 570d8606c27SBarry Smith PetscCallA(SNESSetApplicationContext(snes, user, ierr)) 571c4762a1bSJed Brown 572c4762a1bSJed Brown! Set function evaluation routine and vector 573d8606c27SBarry Smith PetscCallA(SNESSetFunction(snes, r, FormFunction, user, ierr)) 574c4762a1bSJed Brown 575c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 576c4762a1bSJed Brown! Create matrix data structure; set Jacobian evaluation routine 577c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 578c4762a1bSJed Brown 579c4762a1bSJed Brown! Set Jacobian matrix data structure and default Jacobian evaluation 580c4762a1bSJed Brown! routine. User can override with: 581c4762a1bSJed Brown! -snes_fd : default finite differencing approximation of Jacobian 582c4762a1bSJed Brown! -snes_mf : matrix-free Newton-Krylov method with no preconditioning 583c4762a1bSJed Brown! (unless user explicitly sets preconditioner) 5847addb90fSBarry Smith! -snes_mf_operator : form matrix used to construct the preconditioner as set by the user, 585c4762a1bSJed Brown! but use matrix-free approx for Jacobian-vector 586c4762a1bSJed Brown! products within Newton-Krylov method 587c4762a1bSJed Brown! 588c4762a1bSJed Brown! Note: For the parallel case, vectors and matrices MUST be partitioned 589c4762a1bSJed Brown! accordingly. When using distributed arrays (DMDAs) to create vectors, 590c4762a1bSJed Brown! the DMDAs determine the problem partitioning. We must explicitly 591c4762a1bSJed Brown! specify the local matrix dimensions upon its creation for compatibility 592c4762a1bSJed Brown! with the vector distribution. Thus, the generic MatCreate() routine 593c4762a1bSJed Brown! is NOT sufficient when working with distributed arrays. 594c4762a1bSJed Brown! 595c4762a1bSJed Brown! Note: Here we only approximately preallocate storage space for the 596c4762a1bSJed Brown! Jacobian. See the users manual for a discussion of better techniques 597c4762a1bSJed Brown! for preallocating matrix memory. 598c4762a1bSJed Brown 599d8606c27SBarry Smith PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-snes_mf', matrix_free, ierr)) 600c4762a1bSJed Brown if (.not. matrix_free) then 601d8606c27SBarry Smith PetscCallA(DMSetMatType(da, MATAIJ, ierr)) 602d8606c27SBarry Smith PetscCallA(DMCreateMatrix(da, J, ierr)) 603d8606c27SBarry Smith PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, user, ierr)) 604c4762a1bSJed Brown end if 605c4762a1bSJed Brown 606c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 607c4762a1bSJed Brown! Customize nonlinear solver; set runtime options 608c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 609c4762a1bSJed Brown! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 610d8606c27SBarry Smith PetscCallA(SNESSetDM(snes, da, ierr)) 611d8606c27SBarry Smith PetscCallA(SNESSetFromOptions(snes, ierr)) 612c4762a1bSJed Brown 613c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 614c4762a1bSJed Brown! Evaluate initial guess; then solve nonlinear system. 615c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 616c4762a1bSJed Brown! Note: The user should initialize the vector, x, with the initial guess 617c4762a1bSJed Brown! for the nonlinear solver prior to calling SNESSolve(). In particular, 618c4762a1bSJed Brown! to employ an initial guess of zero, the user should explicitly set 619c4762a1bSJed Brown! this vector to zero by calling VecSet(). 620c4762a1bSJed Brown 621d8606c27SBarry Smith PetscCallA(FormInitialGuess(snes, x, ierr)) 622d8606c27SBarry Smith PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr)) 623d8606c27SBarry Smith PetscCallA(SNESGetIterationNumber(snes, its, ierr)) 6244820e4eaSBarry Smith if (user%rank == 0) then 625c4762a1bSJed Brown write (6, 100) its 626c4762a1bSJed Brown end if 627c4762a1bSJed Brown100 format('Number of SNES iterations = ', i5) 628c4762a1bSJed Brown 629c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 630c4762a1bSJed Brown! Free work space. All PETSc objects should be destroyed when they 631c4762a1bSJed Brown! are no longer needed. 632c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 633d8606c27SBarry Smith if (.not. matrix_free) PetscCallA(MatDestroy(J, ierr)) 634d8606c27SBarry Smith PetscCallA(VecDestroy(x, ierr)) 635d8606c27SBarry Smith PetscCallA(VecDestroy(r, ierr)) 636d8606c27SBarry Smith PetscCallA(SNESDestroy(snes, ierr)) 637d8606c27SBarry Smith PetscCallA(DMDestroy(da, ierr)) 638c4762a1bSJed Brown 639d8606c27SBarry Smith PetscCallA(PetscFinalize(ierr)) 640c4762a1bSJed Brownend 641c4762a1bSJed Brown! 642c4762a1bSJed Brown!/*TEST 643c4762a1bSJed Brown! 644c4762a1bSJed Brown! test: 645c4762a1bSJed Brown! nsize: 4 646c4762a1bSJed Brown! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 647c4762a1bSJed Brown! requires: !single 648c4762a1bSJed Brown! 649c4762a1bSJed Brown! test: 650c4762a1bSJed Brown! suffix: 2 651c4762a1bSJed Brown! nsize: 4 652c4762a1bSJed Brown! args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 653c4762a1bSJed Brown! requires: !single 654c4762a1bSJed Brown! 655c4762a1bSJed Brown! test: 656c4762a1bSJed Brown! suffix: 3 657c4762a1bSJed Brown! nsize: 3 658c4762a1bSJed Brown! args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 659c4762a1bSJed Brown! requires: !single 660c4762a1bSJed Brown! 661c4762a1bSJed Brown! test: 662c4762a1bSJed Brown! suffix: 4 663c4762a1bSJed Brown! nsize: 3 664c4762a1bSJed Brown! args: -snes_mf_operator -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 665c4762a1bSJed Brown! requires: !single 666c4762a1bSJed Brown! 667c4762a1bSJed Brown! test: 668c4762a1bSJed Brown! suffix: 5 669c4762a1bSJed Brown! requires: !single 670c4762a1bSJed Brown! 671c4762a1bSJed Brown!TEST*/ 672