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