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