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