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);CHKERRMPI(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 Create matrix data structure; set Jacobian evaluation routine 76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 77 78 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr); 79 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr); 80 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 81 ierr = MatSeqAIJSetPreallocation(J,3,NULL);CHKERRQ(ierr); 82 ierr = MatSetVariableBlockSizes(J,3,lens);CHKERRQ(ierr); 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 ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr); 96 97 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 98 Customize nonlinear solver; set runtime options 99 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 100 101 /* 102 Set names for some vectors to facilitate monitoring (optional) 103 */ 104 ierr = PetscObjectSetName((PetscObject)x,"Approximate Solution");CHKERRQ(ierr); 105 ierr = PetscObjectSetName((PetscObject)U,"Exact Solution");CHKERRQ(ierr); 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 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 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 ierr = SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);CHKERRQ(ierr); 119 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); 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 ierr = VecSetValues(F,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 130 v = xp*xp*xp; 131 ierr = VecSetValues(U,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = FormInitialGuess(x);CHKERRQ(ierr); 145 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); 146 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 147 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its);CHKERRQ(ierr); 148 149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 Check solution and clean up 151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152 153 /* 154 Check the error 155 */ 156 ierr = VecAXPY(x,none,U);CHKERRQ(ierr); 157 ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 158 ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr); 159 160 /* 161 Free work space. All PETSc objects should be destroyed when they 162 are no longer needed. 163 */ 164 ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); 165 ierr = VecDestroy(&U);CHKERRQ(ierr); ierr = VecDestroy(&F);CHKERRQ(ierr); 166 ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); 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 PetscErrorCode ierr; 180 PetscScalar pfive = .50; 181 ierr = VecSet(x,pfive);CHKERRQ(ierr); 182 return 0; 183 } 184 /* ------------------------------------------------------------------- */ 185 /* 186 FormFunction - Evaluates nonlinear function, F(x). 187 188 Input Parameters: 189 . snes - the SNES context 190 . x - input vector 191 . ctx - optional user-defined context, as set by SNESSetFunction() 192 193 Output Parameter: 194 . f - function vector 195 196 Note: 197 The user-defined context can contain any application-specific data 198 needed for the function evaluation (such as various parameters, work 199 vectors, and grid information). In this program the context is just 200 a vector containing the right-hand-side of the discretized PDE. 201 */ 202 203 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx) 204 { 205 Vec g = (Vec)ctx; 206 const PetscScalar *xx,*gg; 207 PetscScalar *ff,d; 208 PetscErrorCode ierr; 209 PetscInt i,n; 210 211 /* 212 Get pointers to vector data. 213 - For default PETSc vectors, VecGetArray() returns a pointer to 214 the data array. Otherwise, the routine is implementation dependent. 215 - You MUST call VecRestoreArray() when you no longer need access to 216 the array. 217 */ 218 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 219 ierr = VecGetArray(f,&ff);CHKERRQ(ierr); 220 ierr = VecGetArrayRead(g,&gg);CHKERRQ(ierr); 221 222 /* 223 Compute function 224 */ 225 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 226 d = (PetscReal)(n - 1); d = d*d; 227 ff[0] = xx[0]; 228 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]; 229 ff[n-1] = xx[n-1] - 1.0; 230 231 /* 232 Restore vectors 233 */ 234 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 235 ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr); 236 ierr = VecRestoreArrayRead(g,&gg);CHKERRQ(ierr); 237 return 0; 238 } 239 /* ------------------------------------------------------------------- */ 240 /* 241 FormJacobian - Evaluates Jacobian matrix. 242 243 Input Parameters: 244 . snes - the SNES context 245 . x - input vector 246 . dummy - optional user-defined context (not used here) 247 248 Output Parameters: 249 . jac - Jacobian matrix 250 . B - optionally different preconditioning matrix 251 252 */ 253 254 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 255 { 256 const PetscScalar *xx; 257 PetscScalar A[3],d; 258 PetscErrorCode ierr; 259 PetscInt i,n,j[3]; 260 261 /* 262 Get pointer to vector data 263 */ 264 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 265 266 /* 267 Compute Jacobian entries and insert into matrix. 268 - Note that in this case we set all elements for a particular 269 row at once. 270 */ 271 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 272 d = (PetscReal)(n - 1); d = d*d; 273 274 /* 275 Interior grid points 276 */ 277 for (i=1; i<n-1; i++) { 278 j[0] = i - 1; j[1] = i; j[2] = i + 1; 279 A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i]; 280 ierr = MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr); 281 } 282 283 /* 284 Boundary points 285 */ 286 i = 0; A[0] = 1.0; 287 288 ierr = MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr); 289 290 i = n-1; A[0] = 1.0; 291 292 ierr = MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr); 293 294 /* 295 Restore vector 296 */ 297 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 298 299 /* 300 Assemble matrix 301 */ 302 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 303 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 304 if (jac != B) { 305 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 306 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 307 } 308 return 0; 309 } 310 311 /*TEST 312 313 test: 314 args: -snes_monitor_short -snes_view -ksp_monitor 315 316 # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly 317 # the solution is wrong on purpose 318 test: 319 requires: !single !complex 320 suffix: transpose_only 321 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 322 323 TEST*/ 324