1 2 static char help[] = "Newton's method for a two-variable system, sequential.\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 /*F 14 This examples solves either 15 \begin{equation} 16 F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{x^2_0 + x_0 x_1 - 3}{x_0 x_1 + x^2_1 - 6} 17 \end{equation} 18 or if the {\tt -hard} options is given 19 \begin{equation} 20 F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1} 21 \end{equation} 22 F*/ 23 #include <petscsnes.h> 24 25 /* 26 User-defined routines 27 */ 28 extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *); 29 extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *); 30 extern PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *); 31 extern PetscErrorCode FormFunction2(SNES, Vec, Vec, void *); 32 33 int main(int argc, char **argv) 34 { 35 SNES snes; /* nonlinear solver context */ 36 KSP ksp; /* linear solver context */ 37 PC pc; /* preconditioner context */ 38 Vec x, r; /* solution, residual vectors */ 39 Mat J; /* Jacobian matrix */ 40 PetscMPIInt size; 41 PetscScalar pfive = .5, *xx; 42 PetscBool flg; 43 44 PetscFunctionBeginUser; 45 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 46 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 47 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Example is only for sequential runs"); 48 49 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 50 Create nonlinear solver context 51 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 52 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 53 PetscCall(SNESSetType(snes, SNESNEWTONLS)); 54 PetscCall(SNESSetOptionsPrefix(snes, "mysolver_")); 55 56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 57 Create matrix and vector data structures; set corresponding routines 58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 59 /* 60 Create vectors for solution and nonlinear function 61 */ 62 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 63 PetscCall(VecSetSizes(x, PETSC_DECIDE, 2)); 64 PetscCall(VecSetFromOptions(x)); 65 PetscCall(VecDuplicate(x, &r)); 66 67 /* 68 Create Jacobian matrix data structure 69 */ 70 PetscCall(MatCreate(PETSC_COMM_WORLD, &J)); 71 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 72 PetscCall(MatSetFromOptions(J)); 73 PetscCall(MatSetUp(J)); 74 75 PetscCall(PetscOptionsHasName(NULL, NULL, "-hard", &flg)); 76 if (!flg) { 77 /* 78 Set function evaluation routine and vector. 79 */ 80 PetscCall(SNESSetFunction(snes, r, FormFunction1, NULL)); 81 82 /* 83 Set Jacobian matrix data structure and Jacobian evaluation routine 84 */ 85 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, NULL)); 86 } else { 87 PetscCall(SNESSetFunction(snes, r, FormFunction2, NULL)); 88 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian2, NULL)); 89 } 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Customize nonlinear solver; set runtime options 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 /* 95 Set linear solver defaults for this problem. By extracting the 96 KSP and PC contexts from the SNES context, we can then 97 directly call any KSP and PC routines to set various options. 98 */ 99 PetscCall(SNESGetKSP(snes, &ksp)); 100 PetscCall(KSPGetPC(ksp, &pc)); 101 PetscCall(PCSetType(pc, PCNONE)); 102 PetscCall(KSPSetTolerances(ksp, 1.e-4, PETSC_DEFAULT, PETSC_DEFAULT, 20)); 103 104 /* 105 Set SNES/KSP/KSP/PC runtime options, e.g., 106 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 107 These options will override those specified above as long as 108 SNESSetFromOptions() is called _after_ any other customization 109 routines. 110 */ 111 PetscCall(SNESSetFromOptions(snes)); 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Evaluate initial guess; then solve nonlinear system 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 if (!flg) { 117 PetscCall(VecSet(x, pfive)); 118 } else { 119 PetscCall(VecGetArray(x, &xx)); 120 xx[0] = 2.0; 121 xx[1] = 3.0; 122 PetscCall(VecRestoreArray(x, &xx)); 123 } 124 /* 125 Note: The user should initialize the vector, x, with the initial guess 126 for the nonlinear solver prior to calling SNESSolve(). In particular, 127 to employ an initial guess of zero, the user should explicitly set 128 this vector to zero by calling VecSet(). 129 */ 130 131 PetscCall(SNESSolve(snes, NULL, x)); 132 if (flg) { 133 Vec f; 134 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 135 PetscCall(SNESGetFunction(snes, &f, 0, 0)); 136 PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD)); 137 } 138 139 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 140 Free work space. All PETSc objects should be destroyed when they 141 are no longer needed. 142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 143 144 PetscCall(VecDestroy(&x)); 145 PetscCall(VecDestroy(&r)); 146 PetscCall(MatDestroy(&J)); 147 PetscCall(SNESDestroy(&snes)); 148 PetscCall(PetscFinalize()); 149 return 0; 150 } 151 /* ------------------------------------------------------------------- */ 152 /* 153 FormFunction1 - Evaluates nonlinear function, F(x). 154 155 Input Parameters: 156 . snes - the SNES context 157 . x - input vector 158 . ctx - optional user-defined context 159 160 Output Parameter: 161 . f - function vector 162 */ 163 PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx) 164 { 165 const PetscScalar *xx; 166 PetscScalar *ff; 167 168 /* 169 Get pointers to vector data. 170 - For default PETSc vectors, VecGetArray() returns a pointer to 171 the data array. Otherwise, the routine is implementation dependent. 172 - You MUST call VecRestoreArray() when you no longer need access to 173 the array. 174 */ 175 PetscCall(VecGetArrayRead(x, &xx)); 176 PetscCall(VecGetArray(f, &ff)); 177 178 /* Compute function */ 179 ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0; 180 ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0; 181 182 /* Restore vectors */ 183 PetscCall(VecRestoreArrayRead(x, &xx)); 184 PetscCall(VecRestoreArray(f, &ff)); 185 return 0; 186 } 187 /* ------------------------------------------------------------------- */ 188 /* 189 FormJacobian1 - Evaluates Jacobian matrix. 190 191 Input Parameters: 192 . snes - the SNES context 193 . x - input vector 194 . dummy - optional user-defined context (not used here) 195 196 Output Parameters: 197 . jac - Jacobian matrix 198 . B - optionally different preconditioning matrix 199 . flag - flag indicating matrix structure 200 */ 201 PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 202 { 203 const PetscScalar *xx; 204 PetscScalar A[4]; 205 PetscInt idx[2] = {0, 1}; 206 207 /* 208 Get pointer to vector data 209 */ 210 PetscCall(VecGetArrayRead(x, &xx)); 211 212 /* 213 Compute Jacobian entries and insert into matrix. 214 - Since this is such a small problem, we set all entries for 215 the matrix at once. 216 */ 217 A[0] = 2.0 * xx[0] + xx[1]; 218 A[1] = xx[0]; 219 A[2] = xx[1]; 220 A[3] = xx[0] + 2.0 * xx[1]; 221 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 222 223 /* 224 Restore vector 225 */ 226 PetscCall(VecRestoreArrayRead(x, &xx)); 227 228 /* 229 Assemble matrix 230 */ 231 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 232 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 233 if (jac != B) { 234 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 235 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 236 } 237 return 0; 238 } 239 240 /* ------------------------------------------------------------------- */ 241 PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy) 242 { 243 const PetscScalar *xx; 244 PetscScalar *ff; 245 246 /* 247 Get pointers to vector data. 248 - For default PETSc vectors, VecGetArray() returns a pointer to 249 the data array. Otherwise, the routine is implementation dependent. 250 - You MUST call VecRestoreArray() when you no longer need access to 251 the array. 252 */ 253 PetscCall(VecGetArrayRead(x, &xx)); 254 PetscCall(VecGetArray(f, &ff)); 255 256 /* 257 Compute function 258 */ 259 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0]; 260 ff[1] = xx[1]; 261 262 /* 263 Restore vectors 264 */ 265 PetscCall(VecRestoreArrayRead(x, &xx)); 266 PetscCall(VecRestoreArray(f, &ff)); 267 return 0; 268 } 269 /* ------------------------------------------------------------------- */ 270 PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 271 { 272 const PetscScalar *xx; 273 PetscScalar A[4]; 274 PetscInt idx[2] = {0, 1}; 275 276 /* 277 Get pointer to vector data 278 */ 279 PetscCall(VecGetArrayRead(x, &xx)); 280 281 /* 282 Compute Jacobian entries and insert into matrix. 283 - Since this is such a small problem, we set all entries for 284 the matrix at once. 285 */ 286 A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0; 287 A[1] = 0.0; 288 A[2] = 0.0; 289 A[3] = 1.0; 290 PetscCall(MatSetValues(B, 2, idx, 2, idx, 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 test: 312 args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop 313 requires: !single 314 315 # test harness puts {{ }} options always at the end, need to specify the prefix explicitly 316 test: 317 suffix: 2 318 requires: !single 319 args: -prefix_push mysolver_ -snes_monitor_short -prefix_pop -mysolver_snes_ksp_ew {{0 1}} 320 output_file: output/ex1_1.out 321 322 test: 323 suffix: 3 324 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short -prefix_pop 325 requires: !single 326 output_file: output/ex1_1.out 327 328 test: 329 suffix: 4 330 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short -prefix_pop 331 requires: !single 332 output_file: output/ex1_1.out 333 334 test: 335 suffix: 5 336 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short -prefix_pop 337 requires: !single 338 output_file: output/ex1_1.out 339 340 test: 341 suffix: 6 342 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short -prefix_pop 343 requires: !single 344 output_file: output/ex1_1.out 345 346 test: 347 suffix: X 348 args: -prefix_push mysolver_ -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 -prefix_pop 349 requires: !single x 350 351 TEST*/ 352