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