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 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 49 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 50 PetscCheckFalse(size > 1,PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs"); 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create nonlinear solver context 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 CHKERRQ(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 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&x)); 64 CHKERRQ(VecSetSizes(x,PETSC_DECIDE,2)); 65 CHKERRQ(VecSetFromOptions(x)); 66 CHKERRQ(VecDuplicate(x,&r)); 67 68 /* 69 Create Jacobian matrix data structure 70 */ 71 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&J)); 72 CHKERRQ(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2)); 73 CHKERRQ(MatSetFromOptions(J)); 74 CHKERRQ(MatSetUp(J)); 75 76 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-hard",&flg)); 77 if (!flg) { 78 /* 79 Set function evaluation routine and vector. 80 */ 81 CHKERRQ(SNESSetFunction(snes,r,FormFunction1,NULL)); 82 83 /* 84 Set Jacobian matrix data structure and Jacobian evaluation routine 85 */ 86 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian1,NULL)); 87 } else { 88 CHKERRQ(SNESSetFunction(snes,r,FormFunction2,NULL)); 89 CHKERRQ(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 CHKERRQ(SNESGetKSP(snes,&ksp)); 101 CHKERRQ(KSPGetPC(ksp,&pc)); 102 CHKERRQ(PCSetType(pc,PCNONE)); 103 CHKERRQ(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 CHKERRQ(SNESSetFromOptions(snes)); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Evaluate initial guess; then solve nonlinear system 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 if (!flg) { 118 CHKERRQ(VecSet(x,pfive)); 119 } else { 120 CHKERRQ(VecGetArray(x,&xx)); 121 xx[0] = 2.0; xx[1] = 3.0; 122 CHKERRQ(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 CHKERRQ(SNESSolve(snes,NULL,x)); 132 if (flg) { 133 Vec f; 134 CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); 135 CHKERRQ(SNESGetFunction(snes,&f,0,0)); 136 CHKERRQ(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 CHKERRQ(VecDestroy(&x)); CHKERRQ(VecDestroy(&r)); 145 CHKERRQ(MatDestroy(&J)); CHKERRQ(SNESDestroy(&snes)); 146 CHKERRQ(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 CHKERRQ(VecGetArrayRead(x,&xx)); 174 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(x,&xx)); 182 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 218 219 /* 220 Restore vector 221 */ 222 CHKERRQ(VecRestoreArrayRead(x,&xx)); 223 224 /* 225 Assemble matrix 226 */ 227 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 228 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 229 if (jac != B) { 230 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 231 CHKERRQ(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 CHKERRQ(VecGetArrayRead(x,&xx)); 250 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(x,&xx)); 262 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES)); 285 286 /* 287 Restore vector 288 */ 289 CHKERRQ(VecRestoreArrayRead(x,&xx)); 290 291 /* 292 Assemble matrix 293 */ 294 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 295 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 296 if (jac != B) { 297 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 298 CHKERRQ(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