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