1 2 static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\ 3 This example employs a user-defined monitoring routine.\n\n"; 4 5 /* 6 Include "petscdraw.h" so that we can use PETSc drawing routines. 7 Include "petscsnes.h" so that we can use SNES solvers. Note that this 8 file automatically includes: 9 petscsys.h - base PETSc routines petscvec.h - vectors 10 petscmat.h - matrices 11 petscis.h - index sets petscksp.h - Krylov subspace methods 12 petscviewer.h - viewers petscpc.h - preconditioners 13 petscksp.h - linear solvers 14 */ 15 16 #include <petscsnes.h> 17 18 /* 19 User-defined routines 20 */ 21 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 22 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 23 extern PetscErrorCode FormInitialGuess(Vec); 24 extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*); 25 26 /* 27 User-defined context for monitoring 28 */ 29 typedef struct { 30 PetscViewer viewer; 31 } MonitorCtx; 32 33 int main(int argc,char **argv) 34 { 35 SNES snes; /* SNES context */ 36 Vec x,r,F,U; /* vectors */ 37 Mat J; /* Jacobian matrix */ 38 MonitorCtx monP; /* monitoring context */ 39 PetscInt its,n = 5,i,maxit,maxf; 40 PetscMPIInt size; 41 PetscScalar h,xp,v,none = -1.0; 42 PetscReal abstol,rtol,stol,norm; 43 44 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 45 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 46 PetscCheck(size == 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 47 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 48 h = 1.0/(n-1); 49 50 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 51 Create nonlinear solver context 52 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 53 54 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 55 56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 57 Create vector data structures; set function evaluation routine 58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 59 /* 60 Note that we form 1 vector from scratch and then duplicate as needed. 61 */ 62 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 63 PetscCall(VecSetSizes(x,PETSC_DECIDE,n)); 64 PetscCall(VecSetFromOptions(x)); 65 PetscCall(VecDuplicate(x,&r)); 66 PetscCall(VecDuplicate(x,&F)); 67 PetscCall(VecDuplicate(x,&U)); 68 69 /* 70 Set function evaluation routine and vector 71 */ 72 PetscCall(SNESSetFunction(snes,r,FormFunction,(void*)F)); 73 74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 75 Create matrix data structure; set Jacobian evaluation routine 76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 77 78 PetscCall(MatCreate(PETSC_COMM_WORLD,&J)); 79 PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n)); 80 PetscCall(MatSetFromOptions(J)); 81 PetscCall(MatSeqAIJSetPreallocation(J,3,NULL)); 82 83 /* 84 Set Jacobian matrix data structure and default Jacobian evaluation 85 routine. User can override with: 86 -snes_fd : default finite differencing approximation of Jacobian 87 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 88 (unless user explicitly sets preconditioner) 89 -snes_mf_operator : form preconditioning matrix as set by the user, 90 but use matrix-free approx for Jacobian-vector 91 products within Newton-Krylov method 92 */ 93 94 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian,NULL)); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Customize nonlinear solver; set runtime options 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 100 /* 101 Set an optional user-defined monitoring routine 102 */ 103 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer)); 104 PetscCall(SNESMonitorSet(snes,Monitor,&monP,0)); 105 106 /* 107 Set names for some vectors to facilitate monitoring (optional) 108 */ 109 PetscCall(PetscObjectSetName((PetscObject)x,"Approximate Solution")); 110 PetscCall(PetscObjectSetName((PetscObject)U,"Exact Solution")); 111 112 /* 113 Set SNES/KSP/KSP/PC runtime options, e.g., 114 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 115 */ 116 PetscCall(SNESSetFromOptions(snes)); 117 118 /* 119 Print parameters used for convergence testing (optional) ... just 120 to demonstrate this routine; this information is also printed with 121 the option -snes_view 122 */ 123 PetscCall(SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf)); 124 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf)); 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Initialize application: 128 Store right-hand-side of PDE and exact solution 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 131 xp = 0.0; 132 for (i=0; i<n; i++) { 133 v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */ 134 PetscCall(VecSetValues(F,1,&i,&v,INSERT_VALUES)); 135 v = xp*xp*xp; 136 PetscCall(VecSetValues(U,1,&i,&v,INSERT_VALUES)); 137 xp += h; 138 } 139 140 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 Evaluate initial guess; then solve nonlinear system 142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 143 /* 144 Note: The user should initialize the vector, x, with the initial guess 145 for the nonlinear solver prior to calling SNESSolve(). In particular, 146 to employ an initial guess of zero, the user should explicitly set 147 this vector to zero by calling VecSet(). 148 */ 149 PetscCall(FormInitialGuess(x)); 150 PetscCall(SNESSolve(snes,NULL,x)); 151 PetscCall(SNESGetIterationNumber(snes,&its)); 152 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %" PetscInt_FMT "\n\n",its)); 153 154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155 Check solution and clean up 156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 157 158 /* 159 Check the error 160 */ 161 PetscCall(VecAXPY(x,none,U)); 162 PetscCall(VecNorm(x,NORM_2,&norm)); 163 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %" PetscInt_FMT "\n",(double)norm,its)); 164 165 /* 166 Free work space. All PETSc objects should be destroyed when they 167 are no longer needed. 168 */ 169 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 170 PetscCall(VecDestroy(&U)); PetscCall(VecDestroy(&F)); 171 PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes)); 172 PetscCall(PetscViewerDestroy(&monP.viewer)); 173 PetscCall(PetscFinalize()); 174 return 0; 175 } 176 /* ------------------------------------------------------------------- */ 177 /* 178 FormInitialGuess - Computes initial guess. 179 180 Input/Output Parameter: 181 . x - the solution vector 182 */ 183 PetscErrorCode FormInitialGuess(Vec x) 184 { 185 PetscScalar pfive = .50; 186 PetscCall(VecSet(x,pfive)); 187 return 0; 188 } 189 /* ------------------------------------------------------------------- */ 190 /* 191 FormFunction - Evaluates nonlinear function, F(x). 192 193 Input Parameters: 194 . snes - the SNES context 195 . x - input vector 196 . ctx - optional user-defined context, as set by SNESSetFunction() 197 198 Output Parameter: 199 . f - function vector 200 201 Note: 202 The user-defined context can contain any application-specific data 203 needed for the function evaluation (such as various parameters, work 204 vectors, and grid information). In this program the context is just 205 a vector containing the right-hand-side of the discretized PDE. 206 */ 207 208 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx) 209 { 210 Vec g = (Vec)ctx; 211 const PetscScalar *xx,*gg; 212 PetscScalar *ff,d; 213 PetscInt i,n; 214 215 /* 216 Get pointers to vector data. 217 - For default PETSc vectors, VecGetArray() returns a pointer to 218 the data array. Otherwise, the routine is implementation dependent. 219 - You MUST call VecRestoreArray() when you no longer need access to 220 the array. 221 */ 222 PetscCall(VecGetArrayRead(x,&xx)); 223 PetscCall(VecGetArray(f,&ff)); 224 PetscCall(VecGetArrayRead(g,&gg)); 225 226 /* 227 Compute function 228 */ 229 PetscCall(VecGetSize(x,&n)); 230 d = (PetscReal)(n - 1); d = d*d; 231 ff[0] = xx[0]; 232 for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i]; 233 ff[n-1] = xx[n-1] - 1.0; 234 235 /* 236 Restore vectors 237 */ 238 PetscCall(VecRestoreArrayRead(x,&xx)); 239 PetscCall(VecRestoreArray(f,&ff)); 240 PetscCall(VecRestoreArrayRead(g,&gg)); 241 return 0; 242 } 243 /* ------------------------------------------------------------------- */ 244 /* 245 FormJacobian - Evaluates Jacobian matrix. 246 247 Input Parameters: 248 . snes - the SNES context 249 . x - input vector 250 . dummy - optional user-defined context (not used here) 251 252 Output Parameters: 253 . jac - Jacobian matrix 254 . B - optionally different preconditioning matrix 255 256 */ 257 258 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 259 { 260 const PetscScalar *xx; 261 PetscScalar A[3],d; 262 PetscInt i,n,j[3]; 263 264 /* 265 Get pointer to vector data 266 */ 267 PetscCall(VecGetArrayRead(x,&xx)); 268 269 /* 270 Compute Jacobian entries and insert into matrix. 271 - Note that in this case we set all elements for a particular 272 row at once. 273 */ 274 PetscCall(VecGetSize(x,&n)); 275 d = (PetscReal)(n - 1); d = d*d; 276 277 /* 278 Interior grid points 279 */ 280 for (i=1; i<n-1; i++) { 281 j[0] = i - 1; j[1] = i; j[2] = i + 1; 282 A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i]; 283 PetscCall(MatSetValues(B,1,&i,3,j,A,INSERT_VALUES)); 284 } 285 286 /* 287 Boundary points 288 */ 289 i = 0; A[0] = 1.0; 290 291 PetscCall(MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES)); 292 293 i = n-1; A[0] = 1.0; 294 295 PetscCall(MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES)); 296 297 /* 298 Restore vector 299 */ 300 PetscCall(VecRestoreArrayRead(x,&xx)); 301 302 /* 303 Assemble matrix 304 */ 305 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 306 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 307 if (jac != B) { 308 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 309 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 310 } 311 return 0; 312 } 313 /* ------------------------------------------------------------------- */ 314 /* 315 Monitor - User-defined monitoring routine that views the 316 current iterate with an x-window plot. 317 318 Input Parameters: 319 snes - the SNES context 320 its - iteration number 321 norm - 2-norm function value (may be estimated) 322 ctx - optional user-defined context for private data for the 323 monitor routine, as set by SNESMonitorSet() 324 325 Note: 326 See the manpage for PetscViewerDrawOpen() for useful runtime options, 327 such as -nox to deactivate all x-window output. 328 */ 329 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx) 330 { 331 MonitorCtx *monP = (MonitorCtx*) ctx; 332 Vec x; 333 334 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"iter = %" PetscInt_FMT ", SNES Function norm %g\n",its,(double)fnorm)); 335 PetscCall(SNESGetSolution(snes,&x)); 336 PetscCall(VecView(x,monP->viewer)); 337 return 0; 338 } 339 340 /*TEST 341 342 test: 343 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always 344 345 test: 346 suffix: 2 347 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view 348 requires: !single 349 350 test: 351 suffix: 3 352 args: -nox -malloc no -options_left no -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always 353 354 test: 355 suffix: 4 356 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view 357 requires: !single 358 359 TEST*/ 360