1 2 static char help[] = "Newton's method for a two-variable system, sequential.\n\n"; 3 4 /*T 5 Concepts: SNES^basic example 6 T*/ 7 8 /* 9 Include "petscsnes.h" so that we can use SNES solvers. Note that this 10 file automatically includes: 11 petscsys.h - base PETSc routines petscvec.h - vectors 12 petscmat.h - matrices 13 petscis.h - index sets petscksp.h - Krylov subspace methods 14 petscviewer.h - viewers petscpc.h - preconditioners 15 petscksp.h - linear solvers 16 */ 17 /*F 18 This examples solves either 19 \begin{equation} 20 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} 21 \end{equation} 22 or if the {\tt -hard} options is given 23 \begin{equation} 24 F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1} 25 \end{equation} 26 F*/ 27 #include <petscsnes.h> 28 29 /* 30 User-defined routines 31 */ 32 extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*); 33 extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*); 34 extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*); 35 extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*); 36 37 int main(int argc,char **argv) 38 { 39 SNES snes; /* nonlinear solver context */ 40 KSP ksp; /* linear solver context */ 41 PC pc; /* preconditioner context */ 42 Vec x,r; /* solution, residual vectors */ 43 Mat J; /* Jacobian matrix */ 44 PetscMPIInt size; 45 PetscScalar pfive = .5,*xx; 46 PetscBool flg; 47 48 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 49 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 50 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Example is only for sequential runs"); 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create nonlinear solver context 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 56 57 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 58 Create matrix and vector data structures; set corresponding routines 59 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 60 /* 61 Create vectors for solution and nonlinear function 62 */ 63 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 64 PetscCall(VecSetSizes(x,PETSC_DECIDE,2)); 65 PetscCall(VecSetFromOptions(x)); 66 PetscCall(VecDuplicate(x,&r)); 67 68 /* 69 Create Jacobian matrix data structure 70 */ 71 PetscCall(MatCreate(PETSC_COMM_WORLD,&J)); 72 PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2)); 73 PetscCall(MatSetFromOptions(J)); 74 PetscCall(MatSetUp(J)); 75 76 PetscCall(PetscOptionsHasName(NULL,NULL,"-hard",&flg)); 77 if (!flg) { 78 /* 79 Set function evaluation routine and vector. 80 */ 81 PetscCall(SNESSetFunction(snes,r,FormFunction1,NULL)); 82 83 /* 84 Set Jacobian matrix data structure and Jacobian evaluation routine 85 */ 86 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian1,NULL)); 87 } else { 88 PetscCall(SNESSetFunction(snes,r,FormFunction2,NULL)); 89 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian2,NULL)); 90 } 91 92 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 93 Customize nonlinear solver; set runtime options 94 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 95 /* 96 Set linear solver defaults for this problem. By extracting the 97 KSP and PC contexts from the SNES context, we can then 98 directly call any KSP and PC routines to set various options. 99 */ 100 PetscCall(SNESGetKSP(snes,&ksp)); 101 PetscCall(KSPGetPC(ksp,&pc)); 102 PetscCall(PCSetType(pc,PCNONE)); 103 PetscCall(KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20)); 104 105 /* 106 Set SNES/KSP/KSP/PC runtime options, e.g., 107 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 108 These options will override those specified above as long as 109 SNESSetFromOptions() is called _after_ any other customization 110 routines. 111 */ 112 PetscCall(SNESSetFromOptions(snes)); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Evaluate initial guess; then solve nonlinear system 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 if (!flg) { 118 PetscCall(VecSet(x,pfive)); 119 } else { 120 PetscCall(VecGetArray(x,&xx)); 121 xx[0] = 2.0; 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)); PetscCall(VecDestroy(&r)); 145 PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes)); 146 PetscCall(PetscFinalize()); 147 return 0; 148 } 149 /* ------------------------------------------------------------------- */ 150 /* 151 FormFunction1 - Evaluates nonlinear function, F(x). 152 153 Input Parameters: 154 . snes - the SNES context 155 . x - input vector 156 . ctx - optional user-defined context 157 158 Output Parameter: 159 . f - function vector 160 */ 161 PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx) 162 { 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 { 201 const PetscScalar *xx; 202 PetscScalar A[4]; 203 PetscInt idx[2] = {0,1}; 204 205 /* 206 Get pointer to vector data 207 */ 208 PetscCall(VecGetArrayRead(x,&xx)); 209 210 /* 211 Compute Jacobian entries and insert into matrix. 212 - Since this is such a small problem, we set all entries for 213 the matrix at once. 214 */ 215 A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0]; 216 A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1]; 217 PetscCall(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 218 219 /* 220 Restore vector 221 */ 222 PetscCall(VecRestoreArrayRead(x,&xx)); 223 224 /* 225 Assemble matrix 226 */ 227 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 228 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 229 if (jac != B) { 230 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 231 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 232 } 233 return 0; 234 } 235 236 /* ------------------------------------------------------------------- */ 237 PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy) 238 { 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 { 268 const PetscScalar *xx; 269 PetscScalar A[4]; 270 PetscInt idx[2] = {0,1}; 271 272 /* 273 Get pointer to vector data 274 */ 275 PetscCall(VecGetArrayRead(x,&xx)); 276 277 /* 278 Compute Jacobian entries and insert into matrix. 279 - Since this is such a small problem, we set all entries for 280 the matrix at once. 281 */ 282 A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0; 283 A[2] = 0.0; A[3] = 1.0; 284 PetscCall(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 285 286 /* 287 Restore vector 288 */ 289 PetscCall(VecRestoreArrayRead(x,&xx)); 290 291 /* 292 Assemble matrix 293 */ 294 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 295 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 296 if (jac != B) { 297 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 298 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 299 } 300 return 0; 301 } 302 303 /*TEST 304 305 test: 306 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 307 requires: !single 308 309 test: 310 suffix: 2 311 requires: !single 312 args: -snes_monitor_short 313 output_file: output/ex1_1.out 314 315 test: 316 suffix: 3 317 args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short 318 requires: !single 319 output_file: output/ex1_1.out 320 321 test: 322 suffix: 4 323 args: -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short 324 requires: !single 325 output_file: output/ex1_1.out 326 327 test: 328 suffix: 5 329 args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short 330 requires: !single 331 output_file: output/ex1_1.out 332 333 test: 334 suffix: 6 335 args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short 336 requires: !single 337 output_file: output/ex1_1.out 338 339 test: 340 suffix: X 341 args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 342 requires: !single x 343 344 TEST*/ 345