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 end 225 226! -------------------- Form initial approximation ----------------- 227 228 subroutine FormInitialGuess(snes,x,ierr) 229 use ex12fmodule 230 implicit none 231 232 PetscErrorCode ierr 233 Vec x 234 SNES snes 235 PetscScalar five 236 237 five = .5 238 PetscCall(VecSet(x,five,ierr)) 239 end 240 241! -------------------- Evaluate Jacobian -------------------- 242 243 subroutine FormJacobian(snes,x,jac,B,ctx,ierr) 244 use ex12fmodule 245 implicit none 246 247 SNES snes 248 Vec x 249 Mat jac,B 250 type(User) ctx 251 PetscInt ii,istart,iend 252 PetscInt i,j,n,end,start,i1 253 PetscErrorCode ierr 254 PetscMPIInt rank,size 255 PetscScalar d,A,h 256 PetscScalar,pointer :: vxx(:) 257 258 i1 = 1 259 h = 1.0/(real(ctx%N) - 1.0) 260 d = h*h 261 PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr)) 262 PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr)) 263 264 PetscCall(VecGetArrayReadF90(x,vxx,ierr)) 265 PetscCall(VecGetOwnershipRange(x,start,end,ierr)) 266 n = end - start 267 268 if (rank .eq. 0) then 269 A = 1.0 270 PetscCall(MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr)) 271 istart = 1 272 else 273 istart = 0 274 endif 275 if (rank .eq. size-1) then 276 i = INT(ctx%N-1) 277 A = 1.0 278 PetscCall(MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr)) 279 iend = n-1 280 else 281 iend = n 282 endif 283 do 10 i=istart,iend-1 284 ii = i + start 285 j = start + i - 1 286 PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)) 287 j = start + i + 1 288 PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)) 289 A = -2.0*d + 2.0*vxx(i+1) 290 PetscCall(MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr)) 291 10 continue 292 PetscCall(VecRestoreArrayReadF90(x,vxx,ierr)) 293 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 294 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 295 end 296 297!/*TEST 298! 299! test: 300! nsize: 2 301! args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short 302! output_file: output/ex12_1.out 303! 304!TEST*/ 305