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 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 51 PetscCheckFalse(size > 1,PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs"); 52 53 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 54 Create nonlinear solver context 55 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 56 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 57 58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 59 Create matrix and vector data structures; set corresponding routines 60 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 61 /* 62 Create vectors for solution and nonlinear function 63 */ 64 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&x)); 65 CHKERRQ(VecSetSizes(x,PETSC_DECIDE,2)); 66 CHKERRQ(VecSetFromOptions(x)); 67 CHKERRQ(VecDuplicate(x,&r)); 68 69 /* 70 Create Jacobian matrix data structure 71 */ 72 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&J)); 73 CHKERRQ(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2)); 74 CHKERRQ(MatSetFromOptions(J)); 75 CHKERRQ(MatSetUp(J)); 76 77 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-hard",&flg)); 78 if (!flg) { 79 /* 80 Set function evaluation routine and vector. 81 */ 82 CHKERRQ(SNESSetFunction(snes,r,FormFunction1,NULL)); 83 84 /* 85 Set Jacobian matrix data structure and Jacobian evaluation routine 86 */ 87 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian1,NULL)); 88 } else { 89 CHKERRQ(SNESSetFunction(snes,r,FormFunction2,NULL)); 90 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian2,NULL)); 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 CHKERRQ(SNESGetKSP(snes,&ksp)); 102 CHKERRQ(KSPGetPC(ksp,&pc)); 103 CHKERRQ(PCSetType(pc,PCNONE)); 104 CHKERRQ(KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20)); 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 CHKERRQ(SNESSetFromOptions(snes)); 114 115 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 116 Evaluate initial guess; then solve nonlinear system 117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 118 if (!flg) { 119 CHKERRQ(VecSet(x,pfive)); 120 } else { 121 CHKERRQ(VecGetArray(x,&xx)); 122 xx[0] = 2.0; xx[1] = 3.0; 123 CHKERRQ(VecRestoreArray(x,&xx)); 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 CHKERRQ(SNESSolve(snes,NULL,x)); 133 if (flg) { 134 Vec f; 135 CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 136 CHKERRQ(SNESGetFunction(snes,&f,0,0)); 137 CHKERRQ(VecView(r,PETSC_VIEWER_STDOUT_WORLD)); 138 } 139 140 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 Free work space. All PETSc objects should be destroyed when they 142 are no longer needed. 143 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 144 145 CHKERRQ(VecDestroy(&x)); CHKERRQ(VecDestroy(&r)); 146 CHKERRQ(MatDestroy(&J)); CHKERRQ(SNESDestroy(&snes)); 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 const PetscScalar *xx; 165 PetscScalar *ff; 166 167 /* 168 Get pointers to vector data. 169 - For default PETSc vectors, VecGetArray() returns a pointer to 170 the data array. Otherwise, the routine is implementation dependent. 171 - You MUST call VecRestoreArray() when you no longer need access to 172 the array. 173 */ 174 CHKERRQ(VecGetArrayRead(x,&xx)); 175 CHKERRQ(VecGetArray(f,&ff)); 176 177 /* Compute function */ 178 ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0; 179 ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0; 180 181 /* Restore vectors */ 182 CHKERRQ(VecRestoreArrayRead(x,&xx)); 183 CHKERRQ(VecRestoreArray(f,&ff)); 184 return 0; 185 } 186 /* ------------------------------------------------------------------- */ 187 /* 188 FormJacobian1 - Evaluates Jacobian matrix. 189 190 Input Parameters: 191 . snes - the SNES context 192 . x - input vector 193 . dummy - optional user-defined context (not used here) 194 195 Output Parameters: 196 . jac - Jacobian matrix 197 . B - optionally different preconditioning matrix 198 . flag - flag indicating matrix structure 199 */ 200 PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 201 { 202 const PetscScalar *xx; 203 PetscScalar A[4]; 204 PetscInt idx[2] = {0,1}; 205 206 /* 207 Get pointer to vector data 208 */ 209 CHKERRQ(VecGetArrayRead(x,&xx)); 210 211 /* 212 Compute Jacobian entries and insert into matrix. 213 - Since this is such a small problem, we set all entries for 214 the matrix at once. 215 */ 216 A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0]; 217 A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1]; 218 CHKERRQ(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 219 220 /* 221 Restore vector 222 */ 223 CHKERRQ(VecRestoreArrayRead(x,&xx)); 224 225 /* 226 Assemble matrix 227 */ 228 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 229 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 230 if (jac != B) { 231 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 232 CHKERRQ(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 233 } 234 return 0; 235 } 236 237 /* ------------------------------------------------------------------- */ 238 PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy) 239 { 240 const PetscScalar *xx; 241 PetscScalar *ff; 242 243 /* 244 Get pointers to vector data. 245 - For default PETSc vectors, VecGetArray() returns a pointer to 246 the data array. Otherwise, the routine is implementation dependent. 247 - You MUST call VecRestoreArray() when you no longer need access to 248 the array. 249 */ 250 CHKERRQ(VecGetArrayRead(x,&xx)); 251 CHKERRQ(VecGetArray(f,&ff)); 252 253 /* 254 Compute function 255 */ 256 ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0]; 257 ff[1] = xx[1]; 258 259 /* 260 Restore vectors 261 */ 262 CHKERRQ(VecRestoreArrayRead(x,&xx)); 263 CHKERRQ(VecRestoreArray(f,&ff)); 264 return 0; 265 } 266 /* ------------------------------------------------------------------- */ 267 PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 268 { 269 const PetscScalar *xx; 270 PetscScalar A[4]; 271 PetscInt idx[2] = {0,1}; 272 273 /* 274 Get pointer to vector data 275 */ 276 CHKERRQ(VecGetArrayRead(x,&xx)); 277 278 /* 279 Compute Jacobian entries and insert into matrix. 280 - Since this is such a small problem, we set all entries for 281 the matrix at once. 282 */ 283 A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0; 284 A[2] = 0.0; A[3] = 1.0; 285 CHKERRQ(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 286 287 /* 288 Restore vector 289 */ 290 CHKERRQ(VecRestoreArrayRead(x,&xx)); 291 292 /* 293 Assemble matrix 294 */ 295 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 296 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 297 if (jac != B) { 298 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 299 CHKERRQ(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 300 } 301 return 0; 302 } 303 304 /*TEST 305 306 test: 307 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 308 requires: !single 309 310 test: 311 suffix: 2 312 requires: !single 313 args: -snes_monitor_short 314 output_file: output/ex1_1.out 315 316 test: 317 suffix: 3 318 args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short 319 requires: !single 320 output_file: output/ex1_1.out 321 322 test: 323 suffix: 4 324 args: -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short 325 requires: !single 326 output_file: output/ex1_1.out 327 328 test: 329 suffix: 5 330 args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short 331 requires: !single 332 output_file: output/ex1_1.out 333 334 test: 335 suffix: 6 336 args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short 337 requires: !single 338 output_file: output/ex1_1.out 339 340 test: 341 suffix: X 342 args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4 343 requires: !single x 344 345 TEST*/ 346