1 static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\ 2 This example tests PCVPBJacobiSetBlocks().\n\n"; 3 4 /* 5 Include "petscsnes.h" so that we can use SNES solvers. Note that this 6 file automatically includes: 7 petscsys.h - base PETSc routines petscvec.h - vectors 8 petscmat.h - matrices 9 petscis.h - index sets petscksp.h - Krylov subspace methods 10 petscviewer.h - viewers petscpc.h - preconditioners 11 petscksp.h - linear solvers 12 */ 13 14 #include <petscsnes.h> 15 16 /* 17 User-defined routines 18 */ 19 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 20 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 21 extern PetscErrorCode FormInitialGuess(Vec); 22 23 int main(int argc, char **argv) 24 { 25 SNES snes; /* SNES context */ 26 Vec x, r, F, U; /* vectors */ 27 Mat J; /* Jacobian matrix */ 28 PetscInt its, n = 5, i, maxit, maxf, lens[3] = {1, 2, 2}; 29 PetscMPIInt size; 30 PetscScalar h, xp, v, none = -1.0; 31 PetscReal abstol, rtol, stol, norm; 32 KSP ksp; 33 PC pc; 34 35 PetscFunctionBeginUser; 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=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\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 = %" PetscInt_FMT "\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 %" PetscInt_FMT "\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)); 159 PetscCall(VecDestroy(&r)); 160 PetscCall(VecDestroy(&U)); 161 PetscCall(VecDestroy(&F)); 162 PetscCall(MatDestroy(&J)); 163 PetscCall(SNESDestroy(&snes)); 164 PetscCall(PetscFinalize()); 165 return 0; 166 } 167 168 /* 169 FormInitialGuess - Computes initial guess. 170 171 Input/Output Parameter: 172 . x - the solution vector 173 */ 174 PetscErrorCode FormInitialGuess(Vec x) 175 { 176 PetscScalar pfive = .50; 177 178 PetscFunctionBeginUser; 179 PetscCall(VecSet(x, pfive)); 180 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 209 /* 210 Get pointers to vector data. 211 - For default PETSc vectors, VecGetArray() returns a pointer to 212 the data array. Otherwise, the routine is implementation dependent. 213 - You MUST call VecRestoreArray() when you no longer need access to 214 the array. 215 */ 216 PetscCall(VecGetArrayRead(x, &xx)); 217 PetscCall(VecGetArray(f, &ff)); 218 PetscCall(VecGetArrayRead(g, &gg)); 219 220 /* 221 Compute function 222 */ 223 PetscCall(VecGetSize(x, &n)); 224 d = (PetscReal)(n - 1); 225 d = d * d; 226 ff[0] = xx[0]; 227 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]; 228 ff[n - 1] = xx[n - 1] - 1.0; 229 230 /* 231 Restore vectors 232 */ 233 PetscCall(VecRestoreArrayRead(x, &xx)); 234 PetscCall(VecRestoreArray(f, &ff)); 235 PetscCall(VecRestoreArrayRead(g, &gg)); 236 PetscFunctionReturn(PETSC_SUCCESS); 237 } 238 239 /* 240 FormJacobian - Evaluates Jacobian matrix. 241 242 Input Parameters: 243 . snes - the SNES context 244 . x - input vector 245 . dummy - optional user-defined context (not used here) 246 247 Output Parameters: 248 . jac - Jacobian matrix 249 . B - optionally different preconditioning matrix 250 251 */ 252 253 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 254 { 255 const PetscScalar *xx; 256 PetscScalar A[3], d; 257 PetscInt i, n, j[3]; 258 259 PetscFunctionBeginUser; 260 /* 261 Get pointer to vector data 262 */ 263 PetscCall(VecGetArrayRead(x, &xx)); 264 265 /* 266 Compute Jacobian entries and insert into matrix. 267 - Note that in this case we set all elements for a particular 268 row at once. 269 */ 270 PetscCall(VecGetSize(x, &n)); 271 d = (PetscReal)(n - 1); 272 d = d * d; 273 274 /* 275 Interior grid points 276 */ 277 for (i = 1; i < n - 1; i++) { 278 j[0] = i - 1; 279 j[1] = i; 280 j[2] = i + 1; 281 A[0] = d; 282 A[1] = -2.0 * d + 2.0 * xx[i]; 283 A[2] = d; 284 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 285 } 286 287 /* 288 Boundary points 289 */ 290 i = 0; 291 A[0] = 1.0; 292 293 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 294 295 i = n - 1; 296 A[0] = 1.0; 297 298 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 299 300 /* 301 Restore vector 302 */ 303 PetscCall(VecRestoreArrayRead(x, &xx)); 304 305 /* 306 Assemble matrix 307 */ 308 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 309 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 310 if (jac != B) { 311 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 312 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 313 } 314 PetscFunctionReturn(PETSC_SUCCESS); 315 } 316 317 /*TEST 318 319 testset: 320 args: -snes_monitor_short -snes_view -ksp_monitor 321 output_file: output/ex5_1.out 322 filter: grep -v "type: seqaij" 323 324 test: 325 suffix: 1 326 327 test: 328 suffix: cuda 329 requires: cuda 330 args: -mat_type aijcusparse -vec_type cuda 331 332 test: 333 suffix: kok 334 requires: kokkos_kernels 335 args: -mat_type aijkokkos -vec_type kokkos 336 337 # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly 338 # the solution is wrong on purpose 339 test: 340 requires: !single !complex 341 suffix: transpose_only 342 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 343 344 TEST*/ 345