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