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