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