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