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