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)); 160 PetscCall(VecDestroy(&r)); 161 PetscCall(VecDestroy(&U)); 162 PetscCall(VecDestroy(&F)); 163 PetscCall(MatDestroy(&J)); 164 PetscCall(SNESDestroy(&snes)); 165 PetscCall(PetscFinalize()); 166 return 0; 167 } 168 169 /* 170 FormInitialGuess - Computes initial guess. 171 172 Input/Output Parameter: 173 . x - the solution vector 174 */ 175 PetscErrorCode FormInitialGuess(Vec x) 176 { 177 PetscScalar pfive = .50; 178 PetscCall(VecSet(x, pfive)); 179 return 0; 180 } 181 182 /* 183 FormFunction - Evaluates nonlinear function, F(x). 184 185 Input Parameters: 186 . snes - the SNES context 187 . x - input vector 188 . ctx - optional user-defined context, as set by SNESSetFunction() 189 190 Output Parameter: 191 . f - function vector 192 193 Note: 194 The user-defined context can contain any application-specific data 195 needed for the function evaluation (such as various parameters, work 196 vectors, and grid information). In this program the context is just 197 a vector containing the right-hand-side of the discretized PDE. 198 */ 199 200 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) 201 { 202 Vec g = (Vec)ctx; 203 const PetscScalar *xx, *gg; 204 PetscScalar *ff, d; 205 PetscInt i, n; 206 207 /* 208 Get pointers to vector data. 209 - For default PETSc vectors, VecGetArray() returns a pointer to 210 the data array. Otherwise, the routine is implementation dependent. 211 - You MUST call VecRestoreArray() when you no longer need access to 212 the array. 213 */ 214 PetscCall(VecGetArrayRead(x, &xx)); 215 PetscCall(VecGetArray(f, &ff)); 216 PetscCall(VecGetArrayRead(g, &gg)); 217 218 /* 219 Compute function 220 */ 221 PetscCall(VecGetSize(x, &n)); 222 d = (PetscReal)(n - 1); 223 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); 269 d = d * d; 270 271 /* 272 Interior grid points 273 */ 274 for (i = 1; i < n - 1; i++) { 275 j[0] = i - 1; 276 j[1] = i; 277 j[2] = i + 1; 278 A[0] = A[2] = d; 279 A[1] = -2.0 * d + 2.0 * xx[i]; 280 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 281 } 282 283 /* 284 Boundary points 285 */ 286 i = 0; 287 A[0] = 1.0; 288 289 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 290 291 i = n - 1; 292 A[0] = 1.0; 293 294 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 295 296 /* 297 Restore vector 298 */ 299 PetscCall(VecRestoreArrayRead(x, &xx)); 300 301 /* 302 Assemble matrix 303 */ 304 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 305 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 306 if (jac != B) { 307 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 308 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 309 } 310 return 0; 311 } 312 313 /*TEST 314 315 testset: 316 args: -snes_monitor_short -snes_view -ksp_monitor 317 output_file: output/ex5_1.out 318 filter: grep -v "type: seqaij" 319 320 test: 321 suffix: 1 322 323 test: 324 suffix: cuda 325 requires: cuda 326 args: -mat_type aijcusparse -vec_type cuda 327 328 test: 329 suffix: kok 330 requires: kokkos_kernels 331 args: -mat_type aijkokkos -vec_type kokkos 332 333 # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly 334 # the solution is wrong on purpose 335 test: 336 requires: !single !complex 337 suffix: transpose_only 338 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 339 340 TEST*/ 341