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