1 static char help[] = "Newton's method for a two-variable system, sequential.\n\n"; 2 3 /* 4 Include "petscsnes.h" so that we can use SNES solvers. Note that this 5 file automatically includes: 6 petscsys.h - base PETSc routines petscvec.h - vectors 7 petscmat.h - matrices 8 petscis.h - index sets petscksp.h - Krylov subspace methods 9 petscviewer.h - viewers petscpc.h - preconditioners 10 petscksp.h - linear solvers 11 */ 12 /*F 13 This examples solves either 14 \begin{equation} 15 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} 16 \end{equation} 17 or if the {\tt -hard} options is given 18 \begin{equation} 19 F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1} 20 \end{equation} 21 F*/ 22 #include <petscsnes.h> 23 24 /* 25 User-defined routines 26 */ 27 extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *); 28 extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *); 29 extern PetscErrorCode FormJacobian2(SNES, Vec, Mat, Mat, void *); 30 extern PetscErrorCode FormFunction2(SNES, Vec, Vec, void *); 31 32 int main(int argc, char **argv) 33 { 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 { 164 const PetscScalar *xx; 165 PetscScalar *ff; 166 167 PetscFunctionBeginUser; 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 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 208 /* 209 Get pointer to vector data 210 */ 211 PetscCall(VecGetArrayRead(x, &xx)); 212 213 /* 214 Compute Jacobian entries and insert into matrix. 215 - Since this is such a small problem, we set all entries for 216 the matrix at once. 217 */ 218 A[0] = 2.0 * xx[0] + xx[1]; 219 A[1] = xx[0]; 220 A[2] = xx[1]; 221 A[3] = xx[0] + 2.0 * xx[1]; 222 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 223 224 /* 225 Restore vector 226 */ 227 PetscCall(VecRestoreArrayRead(x, &xx)); 228 229 /* 230 Assemble matrix 231 */ 232 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 233 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 234 if (jac != B) { 235 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 236 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 237 } 238 PetscFunctionReturn(PETSC_SUCCESS); 239 } 240 241 /* ------------------------------------------------------------------- */ 242 PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy) 243 { 244 const PetscScalar *xx; 245 PetscScalar *ff; 246 247 PetscFunctionBeginUser; 248 /* 249 Get pointers to vector data. 250 - For default PETSc vectors, VecGetArray() returns a pointer to 251 the data array. Otherwise, the routine is implementation dependent. 252 - You MUST call VecRestoreArray() when you no longer need access to 253 the array. 254 */ 255 PetscCall(VecGetArrayRead(x, &xx)); 256 PetscCall(VecGetArray(f, &ff)); 257 258 /* 259 Compute function 260 */ 261 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0]; 262 ff[1] = xx[1]; 263 264 /* 265 Restore vectors 266 */ 267 PetscCall(VecRestoreArrayRead(x, &xx)); 268 PetscCall(VecRestoreArray(f, &ff)); 269 PetscFunctionReturn(PETSC_SUCCESS); 270 } 271 /* ------------------------------------------------------------------- */ 272 PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 273 { 274 const PetscScalar *xx; 275 PetscScalar A[4]; 276 PetscInt idx[2] = {0, 1}; 277 278 PetscFunctionBeginUser; 279 /* 280 Get pointer to vector data 281 */ 282 PetscCall(VecGetArrayRead(x, &xx)); 283 284 /* 285 Compute Jacobian entries and insert into matrix. 286 - Since this is such a small problem, we set all entries for 287 the matrix at once. 288 */ 289 A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0; 290 A[1] = 0.0; 291 A[2] = 0.0; 292 A[3] = 1.0; 293 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 294 295 /* 296 Restore vector 297 */ 298 PetscCall(VecRestoreArrayRead(x, &xx)); 299 300 /* 301 Assemble matrix 302 */ 303 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 304 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 305 if (jac != B) { 306 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 307 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 308 } 309 PetscFunctionReturn(PETSC_SUCCESS); 310 } 311 312 /*TEST 313 314 test: 315 args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop 316 requires: !single 317 318 # test harness puts {{ }} options always at the end, need to specify the prefix explicitly 319 test: 320 suffix: 2 321 requires: !single 322 args: -prefix_push mysolver_ -snes_monitor_short -prefix_pop -mysolver_snes_ksp_ew {{0 1}} 323 output_file: output/ex1_1.out 324 325 test: 326 suffix: 3 327 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short -prefix_pop 328 requires: !single 329 output_file: output/ex1_1.out 330 331 test: 332 suffix: 4 333 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short -prefix_pop 334 requires: !single 335 output_file: output/ex1_1.out 336 337 test: 338 suffix: 5 339 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short -prefix_pop 340 requires: !single 341 output_file: output/ex1_1.out 342 343 test: 344 suffix: 6 345 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short -prefix_pop 346 requires: !single 347 output_file: output/ex1_1.out 348 349 test: 350 suffix: X 351 args: -prefix_push mysolver_ -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 -prefix_pop 352 requires: !single x 353 354 TEST*/ 355