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 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 PetscScalar pfive = .50; 176 PetscCall(VecSet(x, pfive)); 177 return 0; 178 } 179 /* ------------------------------------------------------------------- */ 180 /* 181 FormFunction - Evaluates nonlinear function, F(x). 182 183 Input Parameters: 184 . snes - the SNES context 185 . x - input vector 186 . ctx - optional user-defined context, as set by SNESSetFunction() 187 188 Output Parameter: 189 . f - function vector 190 191 Note: 192 The user-defined context can contain any application-specific data 193 needed for the function evaluation (such as various parameters, work 194 vectors, and grid information). In this program the context is just 195 a vector containing the right-hand-side of the discretized PDE. 196 */ 197 198 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) { 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); 220 d = d * d; 221 ff[0] = xx[0]; 222 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]; 223 ff[n - 1] = xx[n - 1] - 1.0; 224 225 /* 226 Restore vectors 227 */ 228 PetscCall(VecRestoreArrayRead(x, &xx)); 229 PetscCall(VecRestoreArray(f, &ff)); 230 PetscCall(VecRestoreArrayRead(g, &gg)); 231 return 0; 232 } 233 /* ------------------------------------------------------------------- */ 234 /* 235 FormJacobian - Evaluates Jacobian matrix. 236 237 Input Parameters: 238 . snes - the SNES context 239 . x - input vector 240 . dummy - optional user-defined context (not used here) 241 242 Output Parameters: 243 . jac - Jacobian matrix 244 . B - optionally different preconditioning matrix 245 246 */ 247 248 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) { 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); 265 d = d * d; 266 267 /* 268 Interior grid points 269 */ 270 for (i = 1; i < n - 1; i++) { 271 j[0] = i - 1; 272 j[1] = i; 273 j[2] = i + 1; 274 A[0] = A[2] = d; 275 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; 283 A[0] = 1.0; 284 285 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 286 287 i = n - 1; 288 A[0] = 1.0; 289 290 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 291 292 /* 293 Restore vector 294 */ 295 PetscCall(VecRestoreArrayRead(x, &xx)); 296 297 /* 298 Assemble matrix 299 */ 300 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 301 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 302 if (jac != B) { 303 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 304 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 305 } 306 return 0; 307 } 308 309 /*TEST 310 311 testset: 312 args: -snes_monitor_short -snes_view -ksp_monitor 313 output_file: output/ex5_1.out 314 filter: grep -v "type: seqaij" 315 316 test: 317 suffix: 1 318 319 test: 320 suffix: cuda 321 requires: cuda 322 args: -mat_type aijcusparse -vec_type cuda 323 324 test: 325 suffix: kok 326 requires: kokkos_kernels 327 args: -mat_type aijkokkos -vec_type kokkos 328 329 # this is just a test for SNESKSPTRASPOSEONLY and KSPSolveTranspose to behave properly 330 # the solution is wrong on purpose 331 test: 332 requires: !single !complex 333 suffix: transpose_only 334 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 335 336 TEST*/ 337