static char help[] = "Solves the nonlinear system, the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\ This example also illustrates the use of matrix coloring. Runtime options include:\n\ -par , where indicates the problem's nonlinearity\n\ problem SFI: = Bratu parameter (0 <= par <= 6.81)\n\ -mx , where = number of grid points in the x-direction\n\ -my , where = number of grid points in the y-direction\n\n"; /* Solid Fuel Ignition (SFI) problem. This problem is modeled by the partial differential equation -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, with boundary conditions u = 0 for x = 0, x = 1, y = 0, y = 1. A finite difference approximation with the usual 5-point stencil is used to discretize the boundary value problem to obtain a nonlinear system of equations. The parallel version of this code is snes/tutorials/ex5.c */ /* Include "petscsnes.h" so that we can use SNES solvers. Note that this file automatically includes: petscsys.h - base PETSc routines petscvec.h - vectors petscmat.h - matrices petscis.h - index sets petscksp.h - Krylov subspace methods petscviewer.h - viewers petscpc.h - preconditioners petscksp.h - linear solvers */ #include /* User-defined application context - contains data needed by the application-provided call-back routines, FormJacobian() and FormFunction(). */ typedef struct { PetscReal param; /* test problem parameter */ PetscInt mx; /* Discretization in x-direction */ PetscInt my; /* Discretization in y-direction */ } AppCtx; /* User-defined routines */ static PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); static PetscErrorCode FormFunction(SNES, Vec, Vec, void *); static PetscErrorCode FormInitialGuess(AppCtx *, Vec); static PetscErrorCode ConvergenceTest(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *); static PetscErrorCode ConvergenceDestroy(PetscCtxRt); static PetscErrorCode postcheck(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *); static PetscErrorCode monitor_change_deltamax(SNES, PetscInt, PetscReal, void *); int main(int argc, char **argv) { SNES snes; /* nonlinear solver context */ Vec x, r; /* solution, residual vectors */ Mat J; /* Jacobian matrix */ AppCtx user; /* user-defined application context */ PetscInt i, its, N, hist_its[50]; PetscMPIInt size; PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0., history[50]; MatFDColoring fdcoloring; 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, test_tr_deltamax = PETSC_FALSE; KSP ksp; PetscInt *testarray; PetscFunctionBeginUser; PetscCall(PetscInitialize(&argc, &argv, NULL, help)); PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); /* Initialize problem parameters */ user.mx = 4; user.my = 4; user.param = 6.0; PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL)); PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL)); PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL)); PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc", &pc, NULL)); PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lambda is out of range"); N = user.mx * user.my; PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL)); PetscCall(PetscOptionsGetBool(NULL, NULL, "-prune_jacobian", &prunejacobian, NULL)); PetscCall(PetscOptionsGetBool(NULL, NULL, "-null_appctx", &null_appctx, NULL)); PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_tr_deltamax", &test_tr_deltamax, NULL)); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create nonlinear solver context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); if (pc) { PetscCall(SNESSetType(snes, SNESNEWTONTR)); PetscCall(SNESNewtonTRSetPostCheck(snes, postcheck, NULL)); } if (test_tr_deltamax) { PetscCall(SNESSetType(snes, SNESNEWTONTR)); PetscCall(SNESMonitorSet(snes, monitor_change_deltamax, NULL, NULL)); } /* Test application context handling from Python */ if (!null_appctx) PetscCall(SNESSetApplicationContext(snes, (void *)&user)); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create vector data structures; set function evaluation routine - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); PetscCall(VecSetSizes(x, PETSC_DECIDE, N)); PetscCall(VecSetFromOptions(x)); PetscCall(VecDuplicate(x, &r)); /* Set function evaluation routine and vector. Whenever the nonlinear solver needs to evaluate the nonlinear function, it will call this routine. - Note that the final routine argument is the user-defined context that provides application-specific data for the function evaluation routine. */ PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user)); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create matrix data structure; set Jacobian evaluation routine - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Create matrix. Here we only approximately preallocate storage space for the Jacobian. See the users manual for a discussion of better techniques for preallocating matrix memory. */ PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL)); if (!matrix_free) { PetscBool matrix_free_operator = PETSC_FALSE; PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &matrix_free_operator, NULL)); if (matrix_free_operator) matrix_free = PETSC_FALSE; } if (!matrix_free) PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, 5, NULL, &J)); /* This option will cause the Jacobian to be computed via finite differences efficiently using a coloring of the columns of the matrix. */ PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_fd_coloring", &fd_coloring, NULL)); 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"); if (fd_coloring) { ISColoring iscoloring; MatColoring mc; if (prunejacobian) { /* Initialize x with random nonzero values so that the nonzeros in the Jacobian can better reflect the sparsity structure of the Jacobian. */ PetscRandom rctx; PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); PetscCall(PetscRandomSetInterval(rctx, 1.0, 2.0)); PetscCall(VecSetRandom(x, rctx)); PetscCall(PetscRandomDestroy(&rctx)); } /* This initializes the nonzero structure of the Jacobian. This is artificial because clearly if we had a routine to compute the Jacobian we won't need to use finite differences. */ PetscCall(FormJacobian(snes, x, J, J, &user)); /* Color the matrix, i.e. determine groups of columns that share no common rows. These columns in the Jacobian can all be computed simultaneously. */ PetscCall(MatColoringCreate(J, &mc)); PetscCall(MatColoringSetType(mc, MATCOLORINGSL)); PetscCall(MatColoringSetFromOptions(mc)); PetscCall(MatColoringApply(mc, &iscoloring)); PetscCall(MatColoringDestroy(&mc)); /* Create the data structure that SNESComputeJacobianDefaultColor() uses to compute the actual Jacobians via finite differences. */ PetscCall(MatFDColoringCreate(J, iscoloring, &fdcoloring)); PetscCall(MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn *)FormFunction, &user)); PetscCall(MatFDColoringSetFromOptions(fdcoloring)); PetscCall(MatFDColoringSetUp(J, iscoloring, fdcoloring)); /* Tell SNES to use the routine SNESComputeJacobianDefaultColor() to compute Jacobians. */ PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring)); PetscCall(ISColoringDestroy(&iscoloring)); if (prunejacobian) PetscCall(SNESPruneJacobianColor(snes, J, J)); } /* Set Jacobian matrix data structure and default Jacobian evaluation routine. Whenever the nonlinear solver needs to compute the Jacobian matrix, it will call this routine. - Note that the final routine argument is the user-defined context that provides application-specific data for the Jacobian evaluation routine. - The user can override with: -snes_fd : default finite differencing approximation of Jacobian -snes_mf : matrix-free Newton-Krylov method with no preconditioning (unless user explicitly sets preconditioner) -snes_mf_operator : form matrix used to construct the preconditioner as set by the user, but use matrix-free approx for Jacobian-vector products within Newton-Krylov method */ else if (!matrix_free) { PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, (void *)&user)); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Customize nonlinear solver; set runtime options - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Set runtime options (e.g., -snes_monitor -snes_rtol -ksp_type ) */ PetscCall(SNESSetFromOptions(snes)); /* Set array that saves the function norms. This array is intended when the user wants to save the convergence history for later use rather than just to view the function norms via -snes_monitor. */ PetscCall(SNESSetConvergenceHistory(snes, history, hist_its, 50, PETSC_TRUE)); /* Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the user provided test before the specialized test. The convergence context is just an array to test that it gets properly freed at the end */ if (use_convergence_test) { PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(PetscMalloc1(5, &testarray)); PetscCall(KSPSetConvergenceTest(ksp, ConvergenceTest, testarray, ConvergenceDestroy)); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Evaluate initial guess; then solve nonlinear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Note: The user should initialize the vector, x, with the initial guess for the nonlinear solver prior to calling SNESSolve(). In particular, to employ an initial guess of zero, the user should explicitly set this vector to zero by calling VecSet(). */ PetscCall(FormInitialGuess(&user, x)); PetscCall(SNESSolve(snes, NULL, x)); PetscCall(SNESGetIterationNumber(snes, &its)); PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its)); /* Print the convergence history. This is intended just to demonstrate use of the data attained via SNESSetConvergenceHistory(). */ PetscCall(PetscOptionsHasName(NULL, NULL, "-print_history", &flg)); if (flg) { 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])); } /* Test NewtonTR API */ PetscCall(SNESNewtonTRSetTolerances(snes, 1.0, 2.0, 3.0)); PetscCall(SNESNewtonTRSetUpdateParameters(snes, 4.0, 5.0, 6.0, 7.0, 8.0)); PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg)); if (flg) { PetscReal tmp[5]; PetscCall(SNESNewtonTRGetTolerances(snes, &tmp[0], &tmp[1], &tmp[2])); PetscCheck(tmp[0] == 1.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); PetscCheck(tmp[1] == 2.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); PetscCheck(tmp[2] == 3.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); PetscCall(SNESNewtonTRGetUpdateParameters(snes, &tmp[0], &tmp[1], &tmp[2], &tmp[3], &tmp[4])); PetscCheck(tmp[0] == 4.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); PetscCheck(tmp[1] == 5.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); PetscCheck(tmp[2] == 6.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); PetscCheck(tmp[3] == 7.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); PetscCheck(tmp[4] == 8.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Free work space. All PETSc objects should be destroyed when they are no longer needed. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (!matrix_free) PetscCall(MatDestroy(&J)); if (fd_coloring) PetscCall(MatFDColoringDestroy(&fdcoloring)); PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); PetscCall(SNESDestroy(&snes)); PetscCall(PetscFinalize()); return 0; } /* FormInitialGuess - Forms initial approximation. Input Parameters: user - user-defined application context X - vector Output Parameter: X - vector */ PetscErrorCode FormInitialGuess(AppCtx *user, Vec X) { PetscInt i, j, row, mx, my; PetscReal lambda, temp1, temp, hx, hy; PetscScalar *x; PetscFunctionBeginUser; mx = user->mx; my = user->my; lambda = user->param; hx = 1.0 / (PetscReal)(mx - 1); hy = 1.0 / (PetscReal)(my - 1); /* Get a pointer to vector data. - For default PETSc vectors, VecGetArray() returns a pointer to the data array. Otherwise, the routine is implementation dependent. - You MUST call VecRestoreArray() when you no longer need access to the array. */ PetscCall(VecGetArray(X, &x)); temp1 = lambda / (lambda + 1.0); for (j = 0; j < my; j++) { temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy; for (i = 0; i < mx; i++) { row = i + j * mx; if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { x[row] = 0.0; continue; } x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp)); } } /* Restore vector */ PetscCall(VecRestoreArray(X, &x)); PetscFunctionReturn(PETSC_SUCCESS); } /* FormFunction - Evaluates nonlinear function, F(x). Input Parameters: . snes - the SNES context . X - input vector . ptr - optional user-defined context, as set by SNESSetFunction() Output Parameter: . F - function vector */ PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr) { AppCtx *user = (AppCtx *)ptr; PetscInt i, j, row, mx, my; PetscReal two = 2.0, one = 1.0, lambda, hx, hy, hxdhy, hydhx; PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f; const PetscScalar *x; PetscFunctionBeginUser; mx = user->mx; my = user->my; lambda = user->param; hx = one / (PetscReal)(mx - 1); hy = one / (PetscReal)(my - 1); sc = hx * hy; hxdhy = hx / hy; hydhx = hy / hx; /* Get pointers to vector data */ PetscCall(VecGetArrayRead(X, &x)); PetscCall(VecGetArray(F, &f)); /* Compute function */ for (j = 0; j < my; j++) { for (i = 0; i < mx; i++) { row = i + j * mx; if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { f[row] = x[row]; continue; } u = x[row]; ub = x[row - mx]; ul = x[row - 1]; ut = x[row + mx]; ur = x[row + 1]; uxx = (-ur + two * u - ul) * hydhx; uyy = (-ut + two * u - ub) * hxdhy; f[row] = uxx + uyy - sc * lambda * PetscExpScalar(u); } } /* Restore vectors */ PetscCall(VecRestoreArrayRead(X, &x)); PetscCall(VecRestoreArray(F, &f)); PetscFunctionReturn(PETSC_SUCCESS); } /* FormJacobian - Evaluates Jacobian matrix. Input Parameters: . snes - the SNES context . x - input vector . ptr - optional user-defined context, as set by SNESSetJacobian() Output Parameters: . A - Jacobian matrix . B - optionally different matrix used to construct the preconditioner */ PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr) { AppCtx *user = (AppCtx *)ptr; /* user-defined application context */ PetscInt i, j, row, mx, my, col[5]; PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc; const PetscScalar *x; PetscReal hx, hy, hxdhy, hydhx; PetscFunctionBeginUser; mx = user->mx; my = user->my; lambda = user->param; hx = 1.0 / (PetscReal)(mx - 1); hy = 1.0 / (PetscReal)(my - 1); sc = hx * hy; hxdhy = hx / hy; hydhx = hy / hx; /* Get pointer to vector data */ PetscCall(VecGetArrayRead(X, &x)); /* Compute entries of the Jacobian */ for (j = 0; j < my; j++) { for (i = 0; i < mx; i++) { row = i + j * mx; if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { PetscCall(MatSetValues(jac, 1, &row, 1, &row, &one, INSERT_VALUES)); continue; } v[0] = -hxdhy; col[0] = row - mx; v[1] = -hydhx; col[1] = row - 1; v[2] = two * (hydhx + hxdhy) - sc * lambda * PetscExpScalar(x[row]); col[2] = row; v[3] = -hydhx; col[3] = row + 1; v[4] = -hxdhy; col[4] = row + mx; PetscCall(MatSetValues(jac, 1, &row, 5, col, v, INSERT_VALUES)); } } /* Restore vector */ PetscCall(VecRestoreArrayRead(X, &x)); /* Assemble matrix */ PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); if (jac != J) { PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode ConvergenceTest(KSP ksp, PetscInt it, PetscReal nrm, KSPConvergedReason *reason, PetscCtx ctx) { PetscFunctionBeginUser; *reason = KSP_CONVERGED_ITERATING; if (it > 1) { *reason = KSP_CONVERGED_ITS; PetscCall(PetscInfo(NULL, "User provided convergence test returning after 2 iterations\n")); } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode ConvergenceDestroy(PetscCtxRt ctx) { PetscFunctionBeginUser; PetscCall(PetscInfo(NULL, "User provided convergence destroy called\n")); PetscCall(PetscFree(*(void **)ctx)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode postcheck(SNES snes, Vec x, Vec y, Vec w, PetscBool *changed_y, PetscBool *changed_w, PetscCtx ctx) { PetscReal norm; Vec tmp; PetscFunctionBeginUser; PetscCall(VecDuplicate(x, &tmp)); PetscCall(VecWAXPY(tmp, -1.0, x, w)); PetscCall(VecNorm(tmp, NORM_2, &norm)); PetscCall(VecDestroy(&tmp)); PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of search step %g\n", (double)norm)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode monitor_change_deltamax(SNES snes, PetscInt it, PetscReal fnorm, PetscCtx ctx) { PetscFunctionBeginUser; if (it == 0) PetscCall(SNESNewtonTRSetTolerances(snes, PETSC_CURRENT, 0.01, PETSC_CURRENT)); PetscFunctionReturn(PETSC_SUCCESS); } /*TEST build: requires: !single test: args: -ksp_gmres_cgs_refinement_type refine_always test: suffix: 2 args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always test: suffix: 2_trdeltamax_change args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -pc -test_tr_deltamax test: suffix: 2a filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault" args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info requires: defined(PETSC_USE_INFO) test: suffix: 2b filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test" args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info requires: defined(PETSC_USE_INFO) test: suffix: 2c 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 test: suffix: 3 args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always test: suffix: 4 args: -pc -par 6.807 -snes_monitor -snes_converged_reason test: suffix: 5 args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always -prune_jacobian output_file: output/ex1_3.out test: suffix: 6 args: -snes_monitor draw:image:testfile -viewer_view test: suffix: python requires: petsc4py args: -python -snes_type python -snes_python_type ex1.py:MySNES -snes_view -null_appctx {{0 1}separate output} localrunfiles: ex1.py TEST*/