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