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> 35*b06eb4cdSBarry Smithmodule ex5module 36c4762a1bSJed Brown use petscsnes 37dfbbaf82SBarry Smith use petscdmda 38e7a95102SMartin 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 47e7a95102SMartin Diehl interface 48e7a95102SMartin Diehl subroutine SNESSetApplicationContext(snes, ctx, ierr) 49e7a95102SMartin Diehl use petscsnes 50e7a95102SMartin Diehl import userctx 51e7a95102SMartin Diehl implicit none 52e7a95102SMartin Diehl SNES snes 53e7a95102SMartin Diehl type(userctx) ctx 54e7a95102SMartin Diehl PetscErrorCode ierr 55e7a95102SMartin Diehl end subroutine 56e7a95102SMartin Diehl subroutine SNESGetApplicationContext(snes, ctx, ierr) 57e7a95102SMartin Diehl use petscsnes 58e7a95102SMartin Diehl import userctx 59e7a95102SMartin Diehl implicit none 60e7a95102SMartin Diehl SNES snes 61e7a95102SMartin Diehl type(userctx), pointer :: ctx 62e7a95102SMartin Diehl PetscErrorCode ierr 63e7a95102SMartin Diehl end subroutine 64e7a95102SMartin Diehl end interface 65e7a95102SMartin 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 136e7a95102SMartin Diehl 137e7a95102SMartin Diehl! --------------------------------------------------------------------- 138e7a95102SMartin Diehl! 139e7a95102SMartin Diehl! FormInitialGuess - Forms initial approximation. 140e7a95102SMartin Diehl! 141e7a95102SMartin Diehl! Input Parameters: 142e7a95102SMartin Diehl! X - vector 143e7a95102SMartin Diehl! 144e7a95102SMartin Diehl! Output Parameter: 145e7a95102SMartin Diehl! X - vector 146e7a95102SMartin Diehl! 147e7a95102SMartin Diehl! Notes: 148e7a95102SMartin Diehl! This routine serves as a wrapper for the lower-level routine 149e7a95102SMartin Diehl! "InitialGuessLocal", where the actual computations are 150e7a95102SMartin Diehl! done using the standard Fortran style of treating the local 151e7a95102SMartin Diehl! vector data as a multidimensional array over the local mesh. 152e7a95102SMartin Diehl! This routine merely handles ghost point scatters and accesses 153e7a95102SMartin Diehl! the local vector data via VecGetArray() and VecRestoreArray(). 154e7a95102SMartin Diehl! 155e7a95102SMartin Diehl subroutine FormInitialGuess(snes, X, ierr) 156e7a95102SMartin Diehl! Input/output variables: 157e7a95102SMartin Diehl SNES snes 158e7a95102SMartin Diehl type(userctx), pointer:: puser 159e7a95102SMartin Diehl Vec X 160e7a95102SMartin Diehl PetscErrorCode ierr 161e7a95102SMartin Diehl DM da 162e7a95102SMartin Diehl 163e7a95102SMartin Diehl! Declarations for use with local arrays: 164e7a95102SMartin Diehl PetscScalar, pointer :: lx_v(:) 165e7a95102SMartin Diehl 166e7a95102SMartin Diehl ierr = 0 167e7a95102SMartin Diehl PetscCallA(SNESGetDM(snes, da, ierr)) 168e7a95102SMartin Diehl PetscCallA(SNESGetApplicationContext(snes, puser, ierr)) 169e7a95102SMartin Diehl! Get a pointer to vector data. 170e7a95102SMartin Diehl! - For default PETSc vectors, VecGetArray() returns a pointer to 171e7a95102SMartin Diehl! the data array. Otherwise, the routine is implementation dependent. 172e7a95102SMartin Diehl! - You MUST call VecRestoreArray() when you no longer need access to 173e7a95102SMartin Diehl! the array. 174e7a95102SMartin Diehl! - Note that the interface to VecGetArray() differs from VecGetArray(). 175e7a95102SMartin Diehl 176e7a95102SMartin Diehl PetscCallA(VecGetArray(X, lx_v, ierr)) 177e7a95102SMartin Diehl 178e7a95102SMartin Diehl! Compute initial guess over the locally owned part of the grid 179e7a95102SMartin Diehl PetscCallA(InitialGuessLocal(puser, lx_v, ierr)) 180e7a95102SMartin Diehl 181e7a95102SMartin Diehl! Restore vector 182e7a95102SMartin Diehl PetscCallA(VecRestoreArray(X, lx_v, ierr)) 183e7a95102SMartin Diehl 184e7a95102SMartin Diehl! Insert values into global vector 185e7a95102SMartin Diehl 186e7a95102SMartin Diehl end 187e7a95102SMartin Diehl 188e7a95102SMartin Diehl! --------------------------------------------------------------------- 189e7a95102SMartin Diehl! 190e7a95102SMartin Diehl! InitialGuessLocal - Computes initial approximation, called by 191e7a95102SMartin Diehl! the higher level routine FormInitialGuess(). 192e7a95102SMartin Diehl! 193e7a95102SMartin Diehl! Input Parameter: 194e7a95102SMartin Diehl! x - local vector data 195e7a95102SMartin Diehl! 196e7a95102SMartin Diehl! Output Parameters: 197e7a95102SMartin Diehl! x - local vector data 198e7a95102SMartin Diehl! ierr - error code 199e7a95102SMartin Diehl! 200e7a95102SMartin Diehl! Notes: 201e7a95102SMartin Diehl! This routine uses standard Fortran-style computations over a 2-dim array. 202e7a95102SMartin Diehl! 203e7a95102SMartin Diehl subroutine InitialGuessLocal(user, x, ierr) 204e7a95102SMartin Diehl! Input/output variables: 205e7a95102SMartin Diehl type(userctx) user 206e7a95102SMartin Diehl PetscScalar x(user%xs:user%xe, user%ys:user%ye) 207e7a95102SMartin Diehl PetscErrorCode ierr 208e7a95102SMartin Diehl 209e7a95102SMartin Diehl! Local variables: 210e7a95102SMartin Diehl PetscInt i, j 211e7a95102SMartin Diehl PetscReal temp1, temp, hx, hy 212e7a95102SMartin Diehl PetscReal one 213e7a95102SMartin Diehl 214e7a95102SMartin Diehl! Set parameters 215e7a95102SMartin Diehl 216e7a95102SMartin Diehl ierr = 0 217e7a95102SMartin Diehl one = 1.0 218e7a95102SMartin Diehl hx = one/(user%mx - 1) 219e7a95102SMartin Diehl hy = one/(user%my - 1) 220e7a95102SMartin Diehl temp1 = user%lambda/(user%lambda + one) 221e7a95102SMartin Diehl 222e7a95102SMartin Diehl do j = user%ys, user%ye 223e7a95102SMartin Diehl temp = min(j - 1, user%my - j)*hy 224e7a95102SMartin Diehl do i = user%xs, user%xe 225e7a95102SMartin Diehl if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 226e7a95102SMartin Diehl x(i, j) = 0.0 227e7a95102SMartin Diehl else 228e7a95102SMartin Diehl x(i, j) = temp1*sqrt(min(hx*min(i - 1, user%mx - i), temp)) 229e7a95102SMartin Diehl end if 230e7a95102SMartin Diehl end do 231e7a95102SMartin Diehl end do 232e7a95102SMartin Diehl 233e7a95102SMartin Diehl end 234e7a95102SMartin Diehl 235e7a95102SMartin Diehl! --------------------------------------------------------------------- 236e7a95102SMartin Diehl! 237e7a95102SMartin Diehl! FormFunctionLocal - Computes nonlinear function, called by 238e7a95102SMartin Diehl! the higher level routine FormFunction(). 239e7a95102SMartin Diehl! 240e7a95102SMartin Diehl! Input Parameter: 241e7a95102SMartin Diehl! x - local vector data 242e7a95102SMartin Diehl! 243e7a95102SMartin Diehl! Output Parameters: 244e7a95102SMartin Diehl! f - local vector data, f(x) 245e7a95102SMartin Diehl! ierr - error code 246e7a95102SMartin Diehl! 247e7a95102SMartin Diehl! Notes: 248e7a95102SMartin Diehl! This routine uses standard Fortran-style computations over a 2-dim array. 249e7a95102SMartin Diehl! 250e7a95102SMartin Diehl subroutine FormFunctionLocal(x, f, user, ierr) 251e7a95102SMartin Diehl! Input/output variables: 252e7a95102SMartin Diehl type(userctx) user 253e7a95102SMartin Diehl PetscScalar x(user%gxs:user%gxe, user%gys:user%gye) 254e7a95102SMartin Diehl PetscScalar f(user%xs:user%xe, user%ys:user%ye) 255e7a95102SMartin Diehl PetscErrorCode ierr 256e7a95102SMartin Diehl 257e7a95102SMartin Diehl! Local variables: 258e7a95102SMartin Diehl PetscScalar two, one, hx, hy, hxdhy, hydhx, sc 259e7a95102SMartin Diehl PetscScalar u, uxx, uyy 260e7a95102SMartin Diehl PetscInt i, j 261e7a95102SMartin Diehl 262e7a95102SMartin Diehl one = 1.0 263e7a95102SMartin Diehl two = 2.0 264e7a95102SMartin Diehl hx = one/(user%mx - 1) 265e7a95102SMartin Diehl hy = one/(user%my - 1) 266e7a95102SMartin Diehl sc = hx*hy*user%lambda 267e7a95102SMartin Diehl hxdhy = hx/hy 268e7a95102SMartin Diehl hydhx = hy/hx 269e7a95102SMartin Diehl 270e7a95102SMartin Diehl! Compute function over the locally owned part of the grid 271e7a95102SMartin Diehl 272e7a95102SMartin Diehl do j = user%ys, user%ye 273e7a95102SMartin Diehl do i = user%xs, user%xe 274e7a95102SMartin Diehl if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 275e7a95102SMartin Diehl f(i, j) = x(i, j) 276e7a95102SMartin Diehl else 277e7a95102SMartin Diehl u = x(i, j) 278e7a95102SMartin Diehl uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j)) 279e7a95102SMartin Diehl uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1)) 280e7a95102SMartin Diehl f(i, j) = uxx + uyy - sc*exp(u) 281e7a95102SMartin Diehl end if 282e7a95102SMartin Diehl end do 283e7a95102SMartin Diehl end do 284e7a95102SMartin Diehl 285e7a95102SMartin Diehl end 286e7a95102SMartin Diehl 287e7a95102SMartin Diehl! --------------------------------------------------------------------- 288e7a95102SMartin Diehl! 289e7a95102SMartin Diehl! FormJacobian - Evaluates Jacobian matrix. 290e7a95102SMartin Diehl! 291e7a95102SMartin Diehl! Input Parameters: 292e7a95102SMartin Diehl! snes - the SNES context 293e7a95102SMartin Diehl! x - input vector 294e7a95102SMartin Diehl! dummy - optional user-defined context, as set by SNESSetJacobian() 295e7a95102SMartin Diehl! (not used here) 296e7a95102SMartin Diehl! 297e7a95102SMartin Diehl! Output Parameters: 298e7a95102SMartin Diehl! jac - Jacobian matrix 299e7a95102SMartin Diehl! jac_prec - optionally different matrix used to construct the preconditioner (not used here) 300e7a95102SMartin Diehl! 301e7a95102SMartin Diehl! Notes: 302e7a95102SMartin Diehl! This routine serves as a wrapper for the lower-level routine 303e7a95102SMartin Diehl! "FormJacobianLocal", where the actual computations are 304e7a95102SMartin Diehl! done using the standard Fortran style of treating the local 305e7a95102SMartin Diehl! vector data as a multidimensional array over the local mesh. 306e7a95102SMartin Diehl! This routine merely accesses the local vector data via 307e7a95102SMartin Diehl! VecGetArray() and VecRestoreArray(). 308e7a95102SMartin Diehl! 309e7a95102SMartin Diehl! Notes: 310e7a95102SMartin Diehl! Due to grid point reordering with DMDAs, we must always work 311e7a95102SMartin Diehl! with the local grid points, and then transform them to the new 312e7a95102SMartin Diehl! global numbering with the "ltog" mapping 313e7a95102SMartin Diehl! We cannot work directly with the global numbers for the original 314e7a95102SMartin Diehl! uniprocessor grid! 315e7a95102SMartin Diehl! 316e7a95102SMartin Diehl! Two methods are available for imposing this transformation 317e7a95102SMartin Diehl! when setting matrix entries: 318e7a95102SMartin Diehl! (A) MatSetValuesLocal(), using the local ordering (including 319e7a95102SMartin Diehl! ghost points!) 320e7a95102SMartin Diehl! - Set matrix entries using the local ordering 321e7a95102SMartin Diehl! by calling MatSetValuesLocal() 322e7a95102SMartin Diehl! (B) MatSetValues(), using the global ordering 323e7a95102SMartin Diehl 324e7a95102SMartin Diehl! - Set matrix entries using the global ordering by calling 325e7a95102SMartin Diehl! MatSetValues() 326e7a95102SMartin Diehl! Option (A) seems cleaner/easier in many cases, and is the procedure 327e7a95102SMartin Diehl! used in this example. 328e7a95102SMartin Diehl! 329e7a95102SMartin Diehl subroutine FormJacobian(snes, X, jac, jac_prec, user, ierr) 330e7a95102SMartin Diehl! Input/output variables: 331e7a95102SMartin Diehl SNES snes 332e7a95102SMartin Diehl Vec X 333e7a95102SMartin Diehl Mat jac, jac_prec 334e7a95102SMartin Diehl type(userctx) user 335e7a95102SMartin Diehl PetscErrorCode ierr 336e7a95102SMartin Diehl DM da 337e7a95102SMartin Diehl 338e7a95102SMartin Diehl! Declarations for use with local arrays: 339e7a95102SMartin Diehl PetscScalar, pointer :: lx_v(:) 340e7a95102SMartin Diehl Vec localX 341e7a95102SMartin Diehl 342e7a95102SMartin Diehl! Scatter ghost points to local vector, using the 2-step process 343e7a95102SMartin Diehl! DMGlobalToLocalBegin(), DMGlobalToLocalEnd() 344e7a95102SMartin Diehl! Computations can be done while messages are in transition, 345e7a95102SMartin Diehl! by placing code between these two statements. 346e7a95102SMartin Diehl 347e7a95102SMartin Diehl PetscCallA(SNESGetDM(snes, da, ierr)) 348e7a95102SMartin Diehl PetscCallA(DMGetLocalVector(da, localX, ierr)) 349e7a95102SMartin Diehl PetscCallA(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr)) 350e7a95102SMartin Diehl PetscCallA(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr)) 351e7a95102SMartin Diehl 352e7a95102SMartin Diehl! Get a pointer to vector data 353e7a95102SMartin Diehl PetscCallA(VecGetArray(localX, lx_v, ierr)) 354e7a95102SMartin Diehl 355e7a95102SMartin Diehl! Compute entries for the locally owned part of the Jacobian preconditioner. 356e7a95102SMartin Diehl PetscCallA(FormJacobianLocal(lx_v, jac_prec, user, ierr)) 357e7a95102SMartin Diehl 358e7a95102SMartin Diehl! Assemble matrix, using the 2-step process: 359e7a95102SMartin Diehl! MatAssemblyBegin(), MatAssemblyEnd() 360e7a95102SMartin Diehl! Computations can be done while messages are in transition, 361e7a95102SMartin Diehl! by placing code between these two statements. 362e7a95102SMartin Diehl 363e7a95102SMartin Diehl PetscCallA(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)) 364e7a95102SMartin Diehl if (jac /= jac_prec) then 365e7a95102SMartin Diehl PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 366e7a95102SMartin Diehl end if 367e7a95102SMartin Diehl PetscCallA(VecRestoreArray(localX, lx_v, ierr)) 368e7a95102SMartin Diehl PetscCallA(DMRestoreLocalVector(da, localX, ierr)) 369e7a95102SMartin Diehl PetscCallA(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)) 370e7a95102SMartin Diehl if (jac /= jac_prec) then 371e7a95102SMartin Diehl PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 372e7a95102SMartin Diehl end if 373e7a95102SMartin Diehl 374e7a95102SMartin Diehl! Tell the matrix we will never add a new nonzero location to the 375e7a95102SMartin Diehl! matrix. If we do it will generate an error. 376e7a95102SMartin Diehl 377e7a95102SMartin Diehl PetscCallA(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr)) 378e7a95102SMartin Diehl 379e7a95102SMartin Diehl end 380e7a95102SMartin Diehl 381e7a95102SMartin Diehl! --------------------------------------------------------------------- 382e7a95102SMartin Diehl! 383e7a95102SMartin Diehl! FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner, 384e7a95102SMartin Diehl! called by the higher level routine FormJacobian(). 385e7a95102SMartin Diehl! 386e7a95102SMartin Diehl! Input Parameters: 387e7a95102SMartin Diehl! x - local vector data 388e7a95102SMartin Diehl! 389e7a95102SMartin Diehl! Output Parameters: 390e7a95102SMartin Diehl! jac_prec - Jacobian matrix used to compute the preconditioner 391e7a95102SMartin Diehl! ierr - error code 392e7a95102SMartin Diehl! 393e7a95102SMartin Diehl! Notes: 394e7a95102SMartin Diehl! This routine uses standard Fortran-style computations over a 2-dim array. 395e7a95102SMartin Diehl! 396e7a95102SMartin Diehl! Notes: 397e7a95102SMartin Diehl! Due to grid point reordering with DMDAs, we must always work 398e7a95102SMartin Diehl! with the local grid points, and then transform them to the new 399e7a95102SMartin Diehl! global numbering with the "ltog" mapping 400e7a95102SMartin Diehl! We cannot work directly with the global numbers for the original 401e7a95102SMartin Diehl! uniprocessor grid! 402e7a95102SMartin Diehl! 403e7a95102SMartin Diehl! Two methods are available for imposing this transformation 404e7a95102SMartin Diehl! when setting matrix entries: 405e7a95102SMartin Diehl! (A) MatSetValuesLocal(), using the local ordering (including 406e7a95102SMartin Diehl! ghost points!) 407e7a95102SMartin Diehl! - Set matrix entries using the local ordering 408e7a95102SMartin Diehl! by calling MatSetValuesLocal() 409e7a95102SMartin Diehl! (B) MatSetValues(), using the global ordering 410e7a95102SMartin Diehl! - Then apply this map explicitly yourself 411e7a95102SMartin Diehl! - Set matrix entries using the global ordering by calling 412e7a95102SMartin Diehl! MatSetValues() 413e7a95102SMartin Diehl! Option (A) seems cleaner/easier in many cases, and is the procedure 414e7a95102SMartin Diehl! used in this example. 415e7a95102SMartin Diehl! 416e7a95102SMartin Diehl subroutine FormJacobianLocal(x, jac_prec, user, ierr) 417e7a95102SMartin Diehl! Input/output variables: 418e7a95102SMartin Diehl type(userctx) user 419e7a95102SMartin Diehl PetscScalar x(user%gxs:user%gxe, user%gys:user%gye) 420e7a95102SMartin Diehl Mat jac_prec 421e7a95102SMartin Diehl PetscErrorCode ierr 422e7a95102SMartin Diehl 423e7a95102SMartin Diehl! Local variables: 424e7a95102SMartin Diehl PetscInt row, col(5), i, j 425e7a95102SMartin Diehl PetscInt ione, ifive 426e7a95102SMartin Diehl PetscScalar two, one, hx, hy, hxdhy 427e7a95102SMartin Diehl PetscScalar hydhx, sc, v(5) 428e7a95102SMartin Diehl 429e7a95102SMartin Diehl! Set parameters 430e7a95102SMartin Diehl ione = 1 431e7a95102SMartin Diehl ifive = 5 432e7a95102SMartin Diehl one = 1.0 433e7a95102SMartin Diehl two = 2.0 434e7a95102SMartin Diehl hx = one/(user%mx - 1) 435e7a95102SMartin Diehl hy = one/(user%my - 1) 436e7a95102SMartin Diehl sc = hx*hy 437e7a95102SMartin Diehl hxdhy = hx/hy 438e7a95102SMartin Diehl hydhx = hy/hx 439e7a95102SMartin Diehl 440e7a95102SMartin Diehl! Compute entries for the locally owned part of the Jacobian. 441e7a95102SMartin Diehl! - Currently, all PETSc parallel matrix formats are partitioned by 442e7a95102SMartin Diehl! contiguous chunks of rows across the processors. 443e7a95102SMartin Diehl! - Each processor needs to insert only elements that it owns 444e7a95102SMartin Diehl! locally (but any non-local elements will be sent to the 445e7a95102SMartin Diehl! appropriate processor during matrix assembly). 446e7a95102SMartin Diehl! - Here, we set all entries for a particular row at once. 447e7a95102SMartin Diehl! - We can set matrix entries either using either 448e7a95102SMartin Diehl! MatSetValuesLocal() or MatSetValues(), as discussed above. 449e7a95102SMartin Diehl! - Note that MatSetValues() uses 0-based row and column numbers 450e7a95102SMartin Diehl! in Fortran as well as in C. 451e7a95102SMartin Diehl 452e7a95102SMartin Diehl do j = user%ys, user%ye 453e7a95102SMartin Diehl row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1 454e7a95102SMartin Diehl do i = user%xs, user%xe 455e7a95102SMartin Diehl row = row + 1 456e7a95102SMartin Diehl! boundary points 457e7a95102SMartin Diehl if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then 458e7a95102SMartin Diehl col(1) = row 459e7a95102SMartin Diehl v(1) = one 460e7a95102SMartin Diehl PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ione, col, v, INSERT_VALUES, ierr)) 461e7a95102SMartin Diehl! interior grid points 462e7a95102SMartin Diehl else 463e7a95102SMartin Diehl v(1) = -hxdhy 464e7a95102SMartin Diehl v(2) = -hydhx 465e7a95102SMartin Diehl v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i, j)) 466e7a95102SMartin Diehl v(4) = -hydhx 467e7a95102SMartin Diehl v(5) = -hxdhy 468e7a95102SMartin Diehl col(1) = row - user%gxm 469e7a95102SMartin Diehl col(2) = row - 1 470e7a95102SMartin Diehl col(3) = row 471e7a95102SMartin Diehl col(4) = row + 1 472e7a95102SMartin Diehl col(5) = row + user%gxm 473e7a95102SMartin Diehl PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ifive, col, v, INSERT_VALUES, ierr)) 474e7a95102SMartin Diehl end if 475e7a95102SMartin Diehl end do 476e7a95102SMartin Diehl end do 477e7a95102SMartin Diehl 478e7a95102SMartin Diehl end 479e7a95102SMartin Diehl 480*b06eb4cdSBarry Smithend module ex5module 481c4762a1bSJed Brown 482c4762a1bSJed Brownprogram main 483*b06eb4cdSBarry Smith use ex5module 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