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