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