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