1! 2! 3! This example demonstrates basic use of the SNES Fortran interface. 4! 5! 6#include <petsc/finclude/petscsnes.h> 7module ex12fmodule 8 use petscsnes 9 implicit none 10 type User 11 DM da 12 Vec F 13 Vec xl 14 MPIU_Comm comm 15 PetscInt N 16 end type User 17 type monctx 18 PetscInt :: its, lag 19 end type monctx 20 21contains 22! -------------------- Evaluate Function F(x) --------------------- 23 24 subroutine FormFunction(snes, x, f, ctx, ierr) 25 SNES snes 26 Vec x, f 27 type(User) ctx 28 PetscMPIInt rank, size, zero 29 PetscInt i, s, n 30 PetscErrorCode ierr 31 PetscScalar h, d 32 PetscScalar, pointer :: vf2(:), vxx(:), vff(:) 33 34 zero = 0 35 PetscCallMPI(MPI_Comm_rank(ctx%comm, rank, ierr)) 36 PetscCallMPI(MPI_Comm_size(ctx%comm, size, ierr)) 37 h = 1.0/(real(ctx%N) - 1.0) 38 PetscCall(DMGlobalToLocalBegin(ctx%da, x, INSERT_VALUES, ctx%xl, ierr)) 39 PetscCall(DMGlobalToLocalEnd(ctx%da, x, INSERT_VALUES, ctx%xl, ierr)) 40 41 PetscCall(VecGetLocalSize(ctx%xl, n, ierr)) 42 if (n > 1000) then 43 print *, 'Local work array not big enough' 44 call MPI_Abort(PETSC_COMM_WORLD, zero, ierr) 45 end if 46 47 PetscCall(VecGetArrayRead(ctx%xl, vxx, ierr)) 48 PetscCall(VecGetArray(f, vff, ierr)) 49 PetscCall(VecGetArray(ctx%F, vF2, ierr)) 50 51 d = h*h 52 53! 54! Note that the array vxx() was obtained from a ghosted local vector 55! ctx%xl while the array vff() was obtained from the non-ghosted parallel 56! vector F. This is why there is a need for shift variable s. Since vff() 57! does not have locations for the ghost variables we need to index in it 58! slightly different then indexing into vxx(). For example on processor 59! 1 (the second processor) 60! 61! xx(1) xx(2) xx(3) ..... 62! ^^^^^^^ ^^^^^ ^^^^^ 63! ghost value 1st local value 2nd local value 64! 65! ff(1) ff(2) 66! ^^^^^^^ ^^^^^^^ 67! 1st local value 2nd local value 68! 69 if (rank == 0) then 70 s = 0 71 vff(1) = vxx(1) 72 else 73 s = 1 74 end if 75 76 do i = 1, n - 2 77 vff(i - s + 1) = d*(vxx(i) - 2.0*vxx(i + 1) + vxx(i + 2)) + vxx(i + 1)*vxx(i + 1) - vF2(i - s + 1) 78 end do 79 80 if (rank == size - 1) then 81 vff(n - s) = vxx(n) - 1.0 82 end if 83 84 PetscCall(VecRestoreArray(f, vff, ierr)) 85 PetscCall(VecRestoreArrayRead(ctx%xl, vxx, ierr)) 86 PetscCall(VecRestoreArray(ctx%F, vF2, ierr)) 87 end 88 89! --------------------------------------------------------------------- 90! Subroutine FormMonitor 91! This function lets up keep track of the SNES progress at each step 92! In this routine, we determine when the Jacobian is rebuilt with the parameter 'jag' 93! 94! Input Parameters: 95! snes - SNES nonlinear solver context 96! its - current nonlinear iteration, starting from a call of SNESSolve() 97! norm - 2-norm of current residual (may be approximate) 98! snesm - monctx designed module (included in Snesmmod) 99! --------------------------------------------------------------------- 100 subroutine FormMonitor(snes, its, norm, snesm, ierr) 101 102 SNES :: snes 103 PetscInt :: its, one, mone 104 PetscScalar :: norm 105 type(monctx) :: snesm 106 PetscErrorCode :: ierr 107 108! write(6,*) ' ' 109! write(6,*) ' its ',its,snesm%its,'lag', 110! & snesm%lag 111! call flush(6) 112 if (mod(snesm%its, snesm%lag) == 0) then 113 one = 1 114 PetscCall(SNESSetLagJacobian(snes, one, ierr)) ! build jacobian 115 else 116 mone = -1 117 PetscCall(SNESSetLagJacobian(snes, mone, ierr)) ! do NOT build jacobian 118 end if 119 snesm%its = snesm%its + 1 120 end subroutine FormMonitor 121 122! -------------------- Form initial approximation ----------------- 123 124 subroutine FormInitialGuess(snes, x, ierr) 125 126 PetscErrorCode ierr 127 Vec x 128 SNES snes 129 PetscScalar five 130 131 five = .5 132 PetscCall(VecSet(x, five, ierr)) 133 end 134 135! -------------------- Evaluate Jacobian -------------------- 136 137 subroutine FormJacobian(snes, x, jac, B, ctx, ierr) 138 139 SNES snes 140 Vec x 141 Mat jac, B 142 type(User) ctx 143 PetscInt ii, istart, iend 144 PetscInt i, j, n, end, start, i1 145 PetscErrorCode ierr 146 PetscMPIInt rank, size 147 PetscScalar d, A, h 148 PetscScalar, pointer :: vxx(:) 149 150 i1 = 1 151 h = 1.0/(real(ctx%N) - 1.0) 152 d = h*h 153 PetscCallMPI(MPI_Comm_rank(ctx%comm, rank, ierr)) 154 PetscCallMPI(MPI_Comm_size(ctx%comm, size, ierr)) 155 156 PetscCall(VecGetArrayRead(x, vxx, ierr)) 157 PetscCall(VecGetOwnershipRange(x, start, end, ierr)) 158 n = end - start 159 160 if (rank == 0) then 161 A = 1.0 162 PetscCall(MatSetValues(jac, i1, [start], i1, [start], [A], INSERT_VALUES, ierr)) 163 istart = 1 164 else 165 istart = 0 166 end if 167 if (rank == size - 1) then 168 i = INT(ctx%N - 1) 169 A = 1.0 170 PetscCall(MatSetValues(jac, i1, [i], i1, [i], [A], INSERT_VALUES, ierr)) 171 iend = n - 1 172 else 173 iend = n 174 end if 175 do i = istart, iend - 1 176 ii = i + start 177 j = start + i - 1 178 PetscCall(MatSetValues(jac, i1, [ii], i1, [j], [d], INSERT_VALUES, ierr)) 179 j = start + i + 1 180 PetscCall(MatSetValues(jac, i1, [ii], i1, [j], [d], INSERT_VALUES, ierr)) 181 A = -2.0*d + 2.0*vxx(i + 1) 182 PetscCall(MatSetValues(jac, i1, [ii], i1, [ii], [A], INSERT_VALUES, ierr)) 183 end do 184 PetscCall(VecRestoreArrayRead(x, vxx, ierr)) 185 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)) 186 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)) 187 end 188 189end module 190 191program main 192 use ex12fmodule 193 implicit none 194 type(User) ctx 195 PetscMPIInt rank, size 196 PetscErrorCode ierr 197 PetscInt N, start, end, nn, i 198 PetscInt ii, its, i1, i0, i3 199 PetscBool flg 200 SNES snes 201 Mat J 202 Vec x, r, u 203 PetscScalar xp, FF, UU, h 204 character*(10) matrixname 205 type(monctx) :: snesm 206 207 PetscCallA(PetscInitialize(ierr)) 208 i1 = 1 209 i0 = 0 210 i3 = 3 211 N = 10 212 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', N, flg, ierr)) 213 h = 1.0/real(N - 1) 214 ctx%N = N 215 ctx%comm = PETSC_COMM_WORLD 216 217 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr)) 218 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)) 219 220! Set up data structures 221 PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, i1, i1, PETSC_NULL_INTEGER_ARRAY, ctx%da, ierr)) 222 PetscCallA(DMSetFromOptions(ctx%da, ierr)) 223 PetscCallA(DMSetUp(ctx%da, ierr)) 224 PetscCallA(DMCreateGlobalVector(ctx%da, x, ierr)) 225 PetscCallA(DMCreateLocalVector(ctx%da, ctx%xl, ierr)) 226 227 PetscCallA(PetscObjectSetName(x, 'Approximate Solution', ierr)) 228 PetscCallA(VecDuplicate(x, r, ierr)) 229 PetscCallA(VecDuplicate(x, ctx%F, ierr)) 230 PetscCallA(VecDuplicate(x, U, ierr)) 231 PetscCallA(PetscObjectSetName(U, 'Exact Solution', ierr)) 232 233 PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, i3, PETSC_NULL_INTEGER_ARRAY, i0, PETSC_NULL_INTEGER_ARRAY, J, ierr)) 234 PetscCallA(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE, ierr)) 235 PetscCallA(MatGetType(J, matrixname, ierr)) 236 237! Store right-hand side of PDE and exact solution 238 PetscCallA(VecGetOwnershipRange(x, start, end, ierr)) 239 xp = h*start 240 nn = end - start 241 ii = start 242 do i = 0, nn - 1 243 FF = 6.0*xp + (xp + 1.e-12)**6.e0 244 UU = xp*xp*xp 245 PetscCallA(VecSetValues(ctx%F, i1, [ii], [FF], INSERT_VALUES, ierr)) 246 PetscCallA(VecSetValues(U, i1, [ii], [UU], INSERT_VALUES, ierr)) 247 xp = xp + h 248 ii = ii + 1 249 end do 250 PetscCallA(VecAssemblyBegin(ctx%F, ierr)) 251 PetscCallA(VecAssemblyEnd(ctx%F, ierr)) 252 PetscCallA(VecAssemblyBegin(U, ierr)) 253 PetscCallA(VecAssemblyEnd(U, ierr)) 254 255! Create nonlinear solver 256 PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr)) 257 258! Set various routines and options 259 PetscCallA(SNESSetFunction(snes, r, FormFunction, ctx, ierr)) 260 PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, ctx, ierr)) 261 262 snesm%its = 0 263 PetscCallA(SNESGetLagJacobian(snes, snesm%lag, ierr)) 264 PetscCallA(SNESMonitorSet(snes, FormMonitor, snesm, PETSC_NULL_FUNCTION, ierr)) 265 PetscCallA(SNESSetFromOptions(snes, ierr)) 266 267! Solve nonlinear system 268 PetscCallA(FormInitialGuess(snes, x, ierr)) 269 PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr)) 270 PetscCallA(SNESGetIterationNumber(snes, its, ierr)) 271 272! Free work space. All PETSc objects should be destroyed when they 273! are no longer needed. 274 PetscCallA(VecDestroy(x, ierr)) 275 PetscCallA(VecDestroy(ctx%xl, ierr)) 276 PetscCallA(VecDestroy(r, ierr)) 277 PetscCallA(VecDestroy(U, ierr)) 278 PetscCallA(VecDestroy(ctx%F, ierr)) 279 PetscCallA(MatDestroy(J, ierr)) 280 PetscCallA(SNESDestroy(snes, ierr)) 281 PetscCallA(DMDestroy(ctx%da, ierr)) 282 PetscCallA(PetscFinalize(ierr)) 283end 284 285!/*TEST 286! 287! test: 288! nsize: 2 289! args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short 290! output_file: output/ex12_1.out 291! 292!TEST*/ 293