1 static char help[] = "Solves the nonlinear system, the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\ 2 This example also illustrates the use of matrix coloring. Runtime options include:\n\ 3 -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\ 4 problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\ 5 -mx <xg>, where <xg> = number of grid points in the x-direction\n\ 6 -my <yg>, where <yg> = number of grid points in the y-direction\n\n"; 7 8 /* 9 10 Solid Fuel Ignition (SFI) problem. This problem is modeled by 11 the partial differential equation 12 13 -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 14 15 with boundary conditions 16 17 u = 0 for x = 0, x = 1, y = 0, y = 1. 18 19 A finite difference approximation with the usual 5-point stencil 20 is used to discretize the boundary value problem to obtain a nonlinear 21 system of equations. 22 23 The parallel version of this code is snes/tutorials/ex5.c 24 25 */ 26 27 /* 28 Include "petscsnes.h" so that we can use SNES solvers. Note that 29 this file automatically includes: 30 petscsys.h - base PETSc routines petscvec.h - vectors 31 petscmat.h - matrices 32 petscis.h - index sets petscksp.h - Krylov subspace methods 33 petscviewer.h - viewers petscpc.h - preconditioners 34 petscksp.h - linear solvers 35 */ 36 37 #include <petscsnes.h> 38 39 /* 40 User-defined application context - contains data needed by the 41 application-provided call-back routines, FormJacobian() and 42 FormFunction(). 43 */ 44 typedef struct { 45 PetscReal param; /* test problem parameter */ 46 PetscInt mx; /* Discretization in x-direction */ 47 PetscInt my; /* Discretization in y-direction */ 48 } AppCtx; 49 50 /* 51 User-defined routines 52 */ 53 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 54 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 55 extern PetscErrorCode FormInitialGuess(AppCtx *, Vec); 56 extern PetscErrorCode ConvergenceTest(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *); 57 extern PetscErrorCode ConvergenceDestroy(void *); 58 extern PetscErrorCode postcheck(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *); 59 60 int main(int argc, char **argv) 61 { 62 SNES snes; /* nonlinear solver context */ 63 Vec x, r; /* solution, residual vectors */ 64 Mat J; /* Jacobian matrix */ 65 AppCtx user; /* user-defined application context */ 66 PetscInt i, its, N, hist_its[50]; 67 PetscMPIInt size; 68 PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0., history[50]; 69 MatFDColoring fdcoloring; 70 PetscBool matrix_free = PETSC_FALSE, flg, fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE, pc = PETSC_FALSE, prunejacobian = PETSC_FALSE, null_appctx = PETSC_TRUE; 71 KSP ksp; 72 PetscInt *testarray; 73 74 PetscFunctionBeginUser; 75 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 76 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 77 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 78 79 /* 80 Initialize problem parameters 81 */ 82 user.mx = 4; 83 user.my = 4; 84 user.param = 6.0; 85 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL)); 86 PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL)); 87 PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL)); 88 PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc", &pc, NULL)); 89 PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lambda is out of range"); 90 N = user.mx * user.my; 91 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL)); 92 PetscCall(PetscOptionsGetBool(NULL, NULL, "-prune_jacobian", &prunejacobian, NULL)); 93 PetscCall(PetscOptionsGetBool(NULL, NULL, "-null_appctx", &null_appctx, NULL)); 94 95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96 Create nonlinear solver context 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 99 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 100 101 if (pc) { 102 PetscCall(SNESSetType(snes, SNESNEWTONTR)); 103 PetscCall(SNESNewtonTRSetPostCheck(snes, postcheck, NULL)); 104 } 105 106 /* Test application context handling from Python */ 107 if (!null_appctx) { PetscCall(SNESSetApplicationContext(snes, (void *)&user)); } 108 109 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110 Create vector data structures; set function evaluation routine 111 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112 113 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 114 PetscCall(VecSetSizes(x, PETSC_DECIDE, N)); 115 PetscCall(VecSetFromOptions(x)); 116 PetscCall(VecDuplicate(x, &r)); 117 118 /* 119 Set function evaluation routine and vector. Whenever the nonlinear 120 solver needs to evaluate the nonlinear function, it will call this 121 routine. 122 - Note that the final routine argument is the user-defined 123 context that provides application-specific data for the 124 function evaluation routine. 125 */ 126 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user)); 127 128 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 129 Create matrix data structure; set Jacobian evaluation routine 130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 131 132 /* 133 Create matrix. Here we only approximately preallocate storage space 134 for the Jacobian. See the users manual for a discussion of better 135 techniques for preallocating matrix memory. 136 */ 137 PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL)); 138 if (!matrix_free) { 139 PetscBool matrix_free_operator = PETSC_FALSE; 140 PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &matrix_free_operator, NULL)); 141 if (matrix_free_operator) matrix_free = PETSC_FALSE; 142 } 143 if (!matrix_free) PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, 5, NULL, &J)); 144 145 /* 146 This option will cause the Jacobian to be computed via finite differences 147 efficiently using a coloring of the columns of the matrix. 148 */ 149 PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_fd_coloring", &fd_coloring, NULL)); 150 PetscCheck(!matrix_free || !fd_coloring, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Use only one of -snes_mf, -snes_fd_coloring options! You can do -snes_mf_operator -snes_fd_coloring"); 151 152 if (fd_coloring) { 153 ISColoring iscoloring; 154 MatColoring mc; 155 if (prunejacobian) { 156 /* Initialize x with random nonzero values so that the nonzeros in the Jacobian 157 can better reflect the sparsity structure of the Jacobian. */ 158 PetscRandom rctx; 159 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 160 PetscCall(PetscRandomSetInterval(rctx, 1.0, 2.0)); 161 PetscCall(VecSetRandom(x, rctx)); 162 PetscCall(PetscRandomDestroy(&rctx)); 163 } 164 165 /* 166 This initializes the nonzero structure of the Jacobian. This is artificial 167 because clearly if we had a routine to compute the Jacobian we won't need 168 to use finite differences. 169 */ 170 PetscCall(FormJacobian(snes, x, J, J, &user)); 171 172 /* 173 Color the matrix, i.e. determine groups of columns that share no common 174 rows. These columns in the Jacobian can all be computed simultaneously. 175 */ 176 PetscCall(MatColoringCreate(J, &mc)); 177 PetscCall(MatColoringSetType(mc, MATCOLORINGSL)); 178 PetscCall(MatColoringSetFromOptions(mc)); 179 PetscCall(MatColoringApply(mc, &iscoloring)); 180 PetscCall(MatColoringDestroy(&mc)); 181 /* 182 Create the data structure that SNESComputeJacobianDefaultColor() uses 183 to compute the actual Jacobians via finite differences. 184 */ 185 PetscCall(MatFDColoringCreate(J, iscoloring, &fdcoloring)); 186 PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))FormFunction, &user)); 187 PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 188 PetscCall(MatFDColoringSetUp(J, iscoloring, fdcoloring)); 189 /* 190 Tell SNES to use the routine SNESComputeJacobianDefaultColor() 191 to compute Jacobians. 192 */ 193 PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring)); 194 PetscCall(ISColoringDestroy(&iscoloring)); 195 if (prunejacobian) PetscCall(SNESPruneJacobianColor(snes, J, J)); 196 } 197 /* 198 Set Jacobian matrix data structure and default Jacobian evaluation 199 routine. Whenever the nonlinear solver needs to compute the 200 Jacobian matrix, it will call this routine. 201 - Note that the final routine argument is the user-defined 202 context that provides application-specific data for the 203 Jacobian evaluation routine. 204 - The user can override with: 205 -snes_fd : default finite differencing approximation of Jacobian 206 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 207 (unless user explicitly sets preconditioner) 208 -snes_mf_operator : form preconditioning matrix as set by the user, 209 but use matrix-free approx for Jacobian-vector 210 products within Newton-Krylov method 211 */ 212 else if (!matrix_free) { 213 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, (void *)&user)); 214 } 215 216 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 217 Customize nonlinear solver; set runtime options 218 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 219 220 /* 221 Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 222 */ 223 PetscCall(SNESSetFromOptions(snes)); 224 225 /* 226 Set array that saves the function norms. This array is intended 227 when the user wants to save the convergence history for later use 228 rather than just to view the function norms via -snes_monitor. 229 */ 230 PetscCall(SNESSetConvergenceHistory(snes, history, hist_its, 50, PETSC_TRUE)); 231 232 /* 233 Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the 234 user provided test before the specialized test. The convergence context is just an array to 235 test that it gets properly freed at the end 236 */ 237 if (use_convergence_test) { 238 PetscCall(SNESGetKSP(snes, &ksp)); 239 PetscCall(PetscMalloc1(5, &testarray)); 240 PetscCall(KSPSetConvergenceTest(ksp, ConvergenceTest, testarray, ConvergenceDestroy)); 241 } 242 243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244 Evaluate initial guess; then solve nonlinear system 245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 246 /* 247 Note: The user should initialize the vector, x, with the initial guess 248 for the nonlinear solver prior to calling SNESSolve(). In particular, 249 to employ an initial guess of zero, the user should explicitly set 250 this vector to zero by calling VecSet(). 251 */ 252 PetscCall(FormInitialGuess(&user, x)); 253 PetscCall(SNESSolve(snes, NULL, x)); 254 PetscCall(SNESGetIterationNumber(snes, &its)); 255 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its)); 256 257 /* 258 Print the convergence history. This is intended just to demonstrate 259 use of the data attained via SNESSetConvergenceHistory(). 260 */ 261 PetscCall(PetscOptionsHasName(NULL, NULL, "-print_history", &flg)); 262 if (flg) { 263 for (i = 0; i < its + 1; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iteration %" PetscInt_FMT ": Linear iterations %" PetscInt_FMT " Function norm = %g\n", i, hist_its[i], (double)history[i])); 264 } 265 266 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 267 Free work space. All PETSc objects should be destroyed when they 268 are no longer needed. 269 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 270 271 if (!matrix_free) PetscCall(MatDestroy(&J)); 272 if (fd_coloring) PetscCall(MatFDColoringDestroy(&fdcoloring)); 273 PetscCall(VecDestroy(&x)); 274 PetscCall(VecDestroy(&r)); 275 PetscCall(SNESDestroy(&snes)); 276 PetscCall(PetscFinalize()); 277 return 0; 278 } 279 280 /* 281 FormInitialGuess - Forms initial approximation. 282 283 Input Parameters: 284 user - user-defined application context 285 X - vector 286 287 Output Parameter: 288 X - vector 289 */ 290 PetscErrorCode FormInitialGuess(AppCtx *user, Vec X) 291 { 292 PetscInt i, j, row, mx, my; 293 PetscReal lambda, temp1, temp, hx, hy; 294 PetscScalar *x; 295 296 PetscFunctionBeginUser; 297 mx = user->mx; 298 my = user->my; 299 lambda = user->param; 300 301 hx = 1.0 / (PetscReal)(mx - 1); 302 hy = 1.0 / (PetscReal)(my - 1); 303 304 /* 305 Get a pointer to vector data. 306 - For default PETSc vectors, VecGetArray() returns a pointer to 307 the data array. Otherwise, the routine is implementation dependent. 308 - You MUST call VecRestoreArray() when you no longer need access to 309 the array. 310 */ 311 PetscCall(VecGetArray(X, &x)); 312 temp1 = lambda / (lambda + 1.0); 313 for (j = 0; j < my; j++) { 314 temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy; 315 for (i = 0; i < mx; i++) { 316 row = i + j * mx; 317 if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 318 x[row] = 0.0; 319 continue; 320 } 321 x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp)); 322 } 323 } 324 325 /* 326 Restore vector 327 */ 328 PetscCall(VecRestoreArray(X, &x)); 329 PetscFunctionReturn(PETSC_SUCCESS); 330 } 331 332 /* 333 FormFunction - Evaluates nonlinear function, F(x). 334 335 Input Parameters: 336 . snes - the SNES context 337 . X - input vector 338 . ptr - optional user-defined context, as set by SNESSetFunction() 339 340 Output Parameter: 341 . F - function vector 342 */ 343 PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr) 344 { 345 AppCtx *user = (AppCtx *)ptr; 346 PetscInt i, j, row, mx, my; 347 PetscReal two = 2.0, one = 1.0, lambda, hx, hy, hxdhy, hydhx; 348 PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f; 349 const PetscScalar *x; 350 351 PetscFunctionBeginUser; 352 mx = user->mx; 353 my = user->my; 354 lambda = user->param; 355 hx = one / (PetscReal)(mx - 1); 356 hy = one / (PetscReal)(my - 1); 357 sc = hx * hy; 358 hxdhy = hx / hy; 359 hydhx = hy / hx; 360 361 /* 362 Get pointers to vector data 363 */ 364 PetscCall(VecGetArrayRead(X, &x)); 365 PetscCall(VecGetArray(F, &f)); 366 367 /* 368 Compute function 369 */ 370 for (j = 0; j < my; j++) { 371 for (i = 0; i < mx; i++) { 372 row = i + j * mx; 373 if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 374 f[row] = x[row]; 375 continue; 376 } 377 u = x[row]; 378 ub = x[row - mx]; 379 ul = x[row - 1]; 380 ut = x[row + mx]; 381 ur = x[row + 1]; 382 uxx = (-ur + two * u - ul) * hydhx; 383 uyy = (-ut + two * u - ub) * hxdhy; 384 f[row] = uxx + uyy - sc * lambda * PetscExpScalar(u); 385 } 386 } 387 388 /* 389 Restore vectors 390 */ 391 PetscCall(VecRestoreArrayRead(X, &x)); 392 PetscCall(VecRestoreArray(F, &f)); 393 PetscFunctionReturn(PETSC_SUCCESS); 394 } 395 396 /* 397 FormJacobian - Evaluates Jacobian matrix. 398 399 Input Parameters: 400 . snes - the SNES context 401 . x - input vector 402 . ptr - optional user-defined context, as set by SNESSetJacobian() 403 404 Output Parameters: 405 . A - Jacobian matrix 406 . B - optionally different preconditioning matrix 407 . flag - flag indicating matrix structure 408 */ 409 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr) 410 { 411 AppCtx *user = (AppCtx *)ptr; /* user-defined application context */ 412 PetscInt i, j, row, mx, my, col[5]; 413 PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc; 414 const PetscScalar *x; 415 PetscReal hx, hy, hxdhy, hydhx; 416 417 PetscFunctionBeginUser; 418 mx = user->mx; 419 my = user->my; 420 lambda = user->param; 421 hx = 1.0 / (PetscReal)(mx - 1); 422 hy = 1.0 / (PetscReal)(my - 1); 423 sc = hx * hy; 424 hxdhy = hx / hy; 425 hydhx = hy / hx; 426 427 /* 428 Get pointer to vector data 429 */ 430 PetscCall(VecGetArrayRead(X, &x)); 431 432 /* 433 Compute entries of the Jacobian 434 */ 435 for (j = 0; j < my; j++) { 436 for (i = 0; i < mx; i++) { 437 row = i + j * mx; 438 if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 439 PetscCall(MatSetValues(jac, 1, &row, 1, &row, &one, INSERT_VALUES)); 440 continue; 441 } 442 v[0] = -hxdhy; 443 col[0] = row - mx; 444 v[1] = -hydhx; 445 col[1] = row - 1; 446 v[2] = two * (hydhx + hxdhy) - sc * lambda * PetscExpScalar(x[row]); 447 col[2] = row; 448 v[3] = -hydhx; 449 col[3] = row + 1; 450 v[4] = -hxdhy; 451 col[4] = row + mx; 452 PetscCall(MatSetValues(jac, 1, &row, 5, col, v, INSERT_VALUES)); 453 } 454 } 455 456 /* 457 Restore vector 458 */ 459 PetscCall(VecRestoreArrayRead(X, &x)); 460 461 /* 462 Assemble matrix 463 */ 464 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 465 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 466 467 if (jac != J) { 468 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 469 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 470 } 471 PetscFunctionReturn(PETSC_SUCCESS); 472 } 473 474 PetscErrorCode ConvergenceTest(KSP ksp, PetscInt it, PetscReal nrm, KSPConvergedReason *reason, void *ctx) 475 { 476 PetscFunctionBeginUser; 477 *reason = KSP_CONVERGED_ITERATING; 478 if (it > 1) { 479 *reason = KSP_CONVERGED_ITS; 480 PetscCall(PetscInfo(NULL, "User provided convergence test returning after 2 iterations\n")); 481 } 482 PetscFunctionReturn(PETSC_SUCCESS); 483 } 484 485 PetscErrorCode ConvergenceDestroy(void *ctx) 486 { 487 PetscFunctionBeginUser; 488 PetscCall(PetscInfo(NULL, "User provided convergence destroy called\n")); 489 PetscCall(PetscFree(ctx)); 490 PetscFunctionReturn(PETSC_SUCCESS); 491 } 492 493 PetscErrorCode postcheck(SNES snes, Vec x, Vec y, Vec w, PetscBool *changed_y, PetscBool *changed_w, void *ctx) 494 { 495 PetscReal norm; 496 Vec tmp; 497 498 PetscFunctionBeginUser; 499 PetscCall(VecDuplicate(x, &tmp)); 500 PetscCall(VecWAXPY(tmp, -1.0, x, w)); 501 PetscCall(VecNorm(tmp, NORM_2, &norm)); 502 PetscCall(VecDestroy(&tmp)); 503 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of search step %g\n", (double)norm)); 504 PetscFunctionReturn(PETSC_SUCCESS); 505 } 506 507 /*TEST 508 509 build: 510 requires: !single 511 512 test: 513 args: -ksp_gmres_cgs_refinement_type refine_always 514 515 test: 516 suffix: 2 517 args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always 518 519 test: 520 suffix: 2a 521 filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault" 522 args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info 523 requires: defined(PETSC_USE_INFO) 524 525 test: 526 suffix: 2b 527 filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test" 528 args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info 529 requires: defined(PETSC_USE_INFO) 530 531 test: 532 suffix: 2c 533 args: -snes_converged_reason -snes_type newtontr -snes_tr_qn {{same different}separate output} -pc_type mat -snes_view -snes_tr_qn_mat_type lmvmdfp -snes_tr_norm_type infinity 534 535 test: 536 suffix: 3 537 args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always 538 539 test: 540 suffix: 4 541 args: -pc -par 6.807 -snes_monitor -snes_converged_reason 542 543 test: 544 suffix: 5 545 args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always -prune_jacobian 546 output_file: output/ex1_3.out 547 548 test: 549 suffix: 6 550 args: -snes_monitor draw:image:testfile -viewer_view 551 552 test: 553 suffix: python 554 requires: petsc4py 555 args: -python -snes_type python -snes_python_type ex1.py:MySNES -snes_view -null_appctx {{0 1}separate output} 556 localrunfiles: ex1.py 557 558 TEST*/ 559