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