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