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 PetscFunctionBeginUser; 169 /* 170 Get pointers to vector data. 171 - For default PETSc vectors, VecGetArray() returns a pointer to 172 the data array. Otherwise, the routine is implementation dependent. 173 - You MUST call VecRestoreArray() when you no longer need access to 174 the array. 175 */ 176 PetscCall(VecGetArrayRead(x, &xx)); 177 PetscCall(VecGetArray(f, &ff)); 178 179 /* Compute function */ 180 ff[0] = xx[0] * xx[0] + xx[0] * xx[1] - 3.0; 181 ff[1] = xx[0] * xx[1] + xx[1] * xx[1] - 6.0; 182 183 /* Restore vectors */ 184 PetscCall(VecRestoreArrayRead(x, &xx)); 185 PetscCall(VecRestoreArray(f, &ff)); 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 /* ------------------------------------------------------------------- */ 189 /* 190 FormJacobian1 - Evaluates Jacobian matrix. 191 192 Input Parameters: 193 . snes - the SNES context 194 . x - input vector 195 . dummy - optional user-defined context (not used here) 196 197 Output Parameters: 198 . jac - Jacobian matrix 199 . B - optionally different preconditioning matrix 200 . flag - flag indicating matrix structure 201 */ 202 PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 203 { 204 const PetscScalar *xx; 205 PetscScalar A[4]; 206 PetscInt idx[2] = {0, 1}; 207 208 PetscFunctionBeginUser; 209 /* 210 Get pointer to vector data 211 */ 212 PetscCall(VecGetArrayRead(x, &xx)); 213 214 /* 215 Compute Jacobian entries and insert into matrix. 216 - Since this is such a small problem, we set all entries for 217 the matrix at once. 218 */ 219 A[0] = 2.0 * xx[0] + xx[1]; 220 A[1] = xx[0]; 221 A[2] = xx[1]; 222 A[3] = xx[0] + 2.0 * xx[1]; 223 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 224 225 /* 226 Restore vector 227 */ 228 PetscCall(VecRestoreArrayRead(x, &xx)); 229 230 /* 231 Assemble matrix 232 */ 233 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 234 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 235 if (jac != B) { 236 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 237 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 238 } 239 PetscFunctionReturn(PETSC_SUCCESS); 240 } 241 242 /* ------------------------------------------------------------------- */ 243 PetscErrorCode FormFunction2(SNES snes, Vec x, Vec f, void *dummy) 244 { 245 const PetscScalar *xx; 246 PetscScalar *ff; 247 248 PetscFunctionBeginUser; 249 /* 250 Get pointers to vector data. 251 - For default PETSc vectors, VecGetArray() returns a pointer to 252 the data array. Otherwise, the routine is implementation dependent. 253 - You MUST call VecRestoreArray() when you no longer need access to 254 the array. 255 */ 256 PetscCall(VecGetArrayRead(x, &xx)); 257 PetscCall(VecGetArray(f, &ff)); 258 259 /* 260 Compute function 261 */ 262 ff[0] = PetscSinScalar(3.0 * xx[0]) + xx[0]; 263 ff[1] = xx[1]; 264 265 /* 266 Restore vectors 267 */ 268 PetscCall(VecRestoreArrayRead(x, &xx)); 269 PetscCall(VecRestoreArray(f, &ff)); 270 PetscFunctionReturn(PETSC_SUCCESS); 271 } 272 /* ------------------------------------------------------------------- */ 273 PetscErrorCode FormJacobian2(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 274 { 275 const PetscScalar *xx; 276 PetscScalar A[4]; 277 PetscInt idx[2] = {0, 1}; 278 279 PetscFunctionBeginUser; 280 /* 281 Get pointer to vector data 282 */ 283 PetscCall(VecGetArrayRead(x, &xx)); 284 285 /* 286 Compute Jacobian entries and insert into matrix. 287 - Since this is such a small problem, we set all entries for 288 the matrix at once. 289 */ 290 A[0] = 3.0 * PetscCosScalar(3.0 * xx[0]) + 1.0; 291 A[1] = 0.0; 292 A[2] = 0.0; 293 A[3] = 1.0; 294 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES)); 295 296 /* 297 Restore vector 298 */ 299 PetscCall(VecRestoreArrayRead(x, &xx)); 300 301 /* 302 Assemble matrix 303 */ 304 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 305 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 306 if (jac != B) { 307 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 308 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 309 } 310 PetscFunctionReturn(PETSC_SUCCESS); 311 } 312 313 /*TEST 314 315 test: 316 args: -prefix_push mysolver_ -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short -prefix_pop 317 requires: !single 318 319 # test harness puts {{ }} options always at the end, need to specify the prefix explicitly 320 test: 321 suffix: 2 322 requires: !single 323 args: -prefix_push mysolver_ -snes_monitor_short -prefix_pop -mysolver_snes_ksp_ew {{0 1}} 324 output_file: output/ex1_1.out 325 326 test: 327 suffix: 3 328 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short -prefix_pop 329 requires: !single 330 output_file: output/ex1_1.out 331 332 test: 333 suffix: 4 334 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short -prefix_pop 335 requires: !single 336 output_file: output/ex1_1.out 337 338 test: 339 suffix: 5 340 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short -prefix_pop 341 requires: !single 342 output_file: output/ex1_1.out 343 344 test: 345 suffix: 6 346 args: -prefix_push mysolver_ -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short -prefix_pop 347 requires: !single 348 output_file: output/ex1_1.out 349 350 test: 351 suffix: X 352 args: -prefix_push mysolver_ -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 -prefix_pop 353 requires: !single x 354 355 TEST*/ 356