1c4762a1bSJed Brown static char help[] = "Solves the nonlinear system, the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\ 2c4762a1bSJed Brown This example also illustrates the use of matrix coloring. Runtime options include:\n\ 3c4762a1bSJed Brown -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\ 4c4762a1bSJed Brown problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\ 5c4762a1bSJed Brown -mx <xg>, where <xg> = number of grid points in the x-direction\n\ 6c4762a1bSJed Brown -my <yg>, where <yg> = number of grid points in the y-direction\n\n"; 7c4762a1bSJed Brown 8f6dfbefdSBarry Smith /* 9c4762a1bSJed Brown 10c4762a1bSJed Brown Solid Fuel Ignition (SFI) problem. This problem is modeled by 11c4762a1bSJed Brown the partial differential equation 12c4762a1bSJed Brown 13c4762a1bSJed Brown -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 14c4762a1bSJed Brown 15c4762a1bSJed Brown with boundary conditions 16c4762a1bSJed Brown 17c4762a1bSJed Brown u = 0 for x = 0, x = 1, y = 0, y = 1. 18c4762a1bSJed Brown 19c4762a1bSJed Brown A finite difference approximation with the usual 5-point stencil 20c4762a1bSJed Brown is used to discretize the boundary value problem to obtain a nonlinear 21c4762a1bSJed Brown system of equations. 22c4762a1bSJed Brown 23c4762a1bSJed Brown The parallel version of this code is snes/tutorials/ex5.c 24c4762a1bSJed Brown 25f6dfbefdSBarry Smith */ 26c4762a1bSJed Brown 27c4762a1bSJed Brown /* 28c4762a1bSJed Brown Include "petscsnes.h" so that we can use SNES solvers. Note that 29c4762a1bSJed Brown this file automatically includes: 30c4762a1bSJed Brown petscsys.h - base PETSc routines petscvec.h - vectors 31c4762a1bSJed Brown petscmat.h - matrices 32c4762a1bSJed Brown petscis.h - index sets petscksp.h - Krylov subspace methods 33c4762a1bSJed Brown petscviewer.h - viewers petscpc.h - preconditioners 34c4762a1bSJed Brown petscksp.h - linear solvers 35c4762a1bSJed Brown */ 36c4762a1bSJed Brown 37c4762a1bSJed Brown #include <petscsnes.h> 38c4762a1bSJed Brown 39c4762a1bSJed Brown /* 40c4762a1bSJed Brown User-defined application context - contains data needed by the 41c4762a1bSJed Brown application-provided call-back routines, FormJacobian() and 42c4762a1bSJed Brown FormFunction(). 43c4762a1bSJed Brown */ 44c4762a1bSJed Brown typedef struct { 45c4762a1bSJed Brown PetscReal param; /* test problem parameter */ 46c4762a1bSJed Brown PetscInt mx; /* Discretization in x-direction */ 47c4762a1bSJed Brown PetscInt my; /* Discretization in y-direction */ 48c4762a1bSJed Brown } AppCtx; 49c4762a1bSJed Brown 50c4762a1bSJed Brown /* 51c4762a1bSJed Brown User-defined routines 52c4762a1bSJed Brown */ 53c4762a1bSJed Brown extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 54c4762a1bSJed Brown extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 55c4762a1bSJed Brown extern PetscErrorCode FormInitialGuess(AppCtx *, Vec); 56c4762a1bSJed Brown extern PetscErrorCode ConvergenceTest(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *); 57c4762a1bSJed Brown extern PetscErrorCode ConvergenceDestroy(void *); 58c4762a1bSJed Brown extern PetscErrorCode postcheck(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *); 59c4762a1bSJed Brown 60d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv) 61d71ae5a4SJacob Faibussowitsch { 62c4762a1bSJed Brown SNES snes; /* nonlinear solver context */ 63c4762a1bSJed Brown Vec x, r; /* solution, residual vectors */ 64c4762a1bSJed Brown Mat J; /* Jacobian matrix */ 65c4762a1bSJed Brown AppCtx user; /* user-defined application context */ 66c4762a1bSJed Brown PetscInt i, its, N, hist_its[50]; 67c4762a1bSJed Brown PetscMPIInt size; 68c4762a1bSJed Brown PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0., history[50]; 69c4762a1bSJed Brown MatFDColoring fdcoloring; 709be84c52SStefano Zampini 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; 71c4762a1bSJed Brown KSP ksp; 72c4762a1bSJed Brown PetscInt *testarray; 73c4762a1bSJed Brown 74327415f7SBarry Smith PetscFunctionBeginUser; 75c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 769566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 77be096a46SBarry Smith PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 78c4762a1bSJed Brown 79c4762a1bSJed Brown /* 80c4762a1bSJed Brown Initialize problem parameters 81c4762a1bSJed Brown */ 829371c9d4SSatish Balay user.mx = 4; 839371c9d4SSatish Balay user.my = 4; 849371c9d4SSatish Balay user.param = 6.0; 859566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL)); 869566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL)); 879566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL)); 889566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc", &pc, NULL)); 89e00437b9SBarry Smith PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lambda is out of range"); 90c4762a1bSJed Brown N = user.mx * user.my; 919566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL)); 9278e7fe0eSHong Zhang PetscCall(PetscOptionsGetBool(NULL, NULL, "-prune_jacobian", &prunejacobian, NULL)); 939be84c52SStefano Zampini PetscCall(PetscOptionsGetBool(NULL, NULL, "-null_appctx", &null_appctx, NULL)); 94c4762a1bSJed Brown 95c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96c4762a1bSJed Brown Create nonlinear solver context 97c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98c4762a1bSJed Brown 999566063dSJacob Faibussowitsch PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 100c4762a1bSJed Brown 101c4762a1bSJed Brown if (pc) { 1029566063dSJacob Faibussowitsch PetscCall(SNESSetType(snes, SNESNEWTONTR)); 1039566063dSJacob Faibussowitsch PetscCall(SNESNewtonTRSetPostCheck(snes, postcheck, NULL)); 104c4762a1bSJed Brown } 105c4762a1bSJed Brown 1069be84c52SStefano Zampini /* Test application context handling from Python */ 1079be84c52SStefano Zampini if (!null_appctx) { PetscCall(SNESSetApplicationContext(snes, (void *)&user)); } 1089be84c52SStefano Zampini 109c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 110c4762a1bSJed Brown Create vector data structures; set function evaluation routine 111c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 112c4762a1bSJed Brown 1139566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 1149566063dSJacob Faibussowitsch PetscCall(VecSetSizes(x, PETSC_DECIDE, N)); 1159566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(x)); 1169566063dSJacob Faibussowitsch PetscCall(VecDuplicate(x, &r)); 117c4762a1bSJed Brown 118c4762a1bSJed Brown /* 119c4762a1bSJed Brown Set function evaluation routine and vector. Whenever the nonlinear 120c4762a1bSJed Brown solver needs to evaluate the nonlinear function, it will call this 121c4762a1bSJed Brown routine. 122c4762a1bSJed Brown - Note that the final routine argument is the user-defined 123c4762a1bSJed Brown context that provides application-specific data for the 124c4762a1bSJed Brown function evaluation routine. 125c4762a1bSJed Brown */ 1269566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user)); 127c4762a1bSJed Brown 128c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 129c4762a1bSJed Brown Create matrix data structure; set Jacobian evaluation routine 130c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 131c4762a1bSJed Brown 132c4762a1bSJed Brown /* 133c4762a1bSJed Brown Create matrix. Here we only approximately preallocate storage space 134c4762a1bSJed Brown for the Jacobian. See the users manual for a discussion of better 135c4762a1bSJed Brown techniques for preallocating matrix memory. 136c4762a1bSJed Brown */ 1379566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL)); 138c4762a1bSJed Brown if (!matrix_free) { 139c4762a1bSJed Brown PetscBool matrix_free_operator = PETSC_FALSE; 1409566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &matrix_free_operator, NULL)); 141c4762a1bSJed Brown if (matrix_free_operator) matrix_free = PETSC_FALSE; 142c4762a1bSJed Brown } 14348a46eb9SPierre Jolivet if (!matrix_free) PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, 5, NULL, &J)); 144c4762a1bSJed Brown 145c4762a1bSJed Brown /* 146c4762a1bSJed Brown This option will cause the Jacobian to be computed via finite differences 147c4762a1bSJed Brown efficiently using a coloring of the columns of the matrix. 148c4762a1bSJed Brown */ 1499566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_fd_coloring", &fd_coloring, NULL)); 15000045ab3SPierre Jolivet 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"); 151c4762a1bSJed Brown 152c4762a1bSJed Brown if (fd_coloring) { 153c4762a1bSJed Brown ISColoring iscoloring; 154c4762a1bSJed Brown MatColoring mc; 15578e7fe0eSHong Zhang if (prunejacobian) { 15678e7fe0eSHong Zhang /* Initialize x with random nonzero values so that the nonzeros in the Jacobian 15778e7fe0eSHong Zhang can better reflect the sparsity structure of the Jacobian. */ 15878e7fe0eSHong Zhang PetscRandom rctx; 15978e7fe0eSHong Zhang PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 16078e7fe0eSHong Zhang PetscCall(PetscRandomSetInterval(rctx, 1.0, 2.0)); 16178e7fe0eSHong Zhang PetscCall(VecSetRandom(x, rctx)); 16278e7fe0eSHong Zhang PetscCall(PetscRandomDestroy(&rctx)); 16378e7fe0eSHong Zhang } 164c4762a1bSJed Brown 165c4762a1bSJed Brown /* 166c4762a1bSJed Brown This initializes the nonzero structure of the Jacobian. This is artificial 167c4762a1bSJed Brown because clearly if we had a routine to compute the Jacobian we won't need 168c4762a1bSJed Brown to use finite differences. 169c4762a1bSJed Brown */ 1709566063dSJacob Faibussowitsch PetscCall(FormJacobian(snes, x, J, J, &user)); 171c4762a1bSJed Brown 172c4762a1bSJed Brown /* 173c4762a1bSJed Brown Color the matrix, i.e. determine groups of columns that share no common 174a5b23f4aSJose E. Roman rows. These columns in the Jacobian can all be computed simultaneously. 175c4762a1bSJed Brown */ 1769566063dSJacob Faibussowitsch PetscCall(MatColoringCreate(J, &mc)); 1779566063dSJacob Faibussowitsch PetscCall(MatColoringSetType(mc, MATCOLORINGSL)); 1789566063dSJacob Faibussowitsch PetscCall(MatColoringSetFromOptions(mc)); 1799566063dSJacob Faibussowitsch PetscCall(MatColoringApply(mc, &iscoloring)); 1809566063dSJacob Faibussowitsch PetscCall(MatColoringDestroy(&mc)); 181c4762a1bSJed Brown /* 182c4762a1bSJed Brown Create the data structure that SNESComputeJacobianDefaultColor() uses 183c4762a1bSJed Brown to compute the actual Jacobians via finite differences. 184c4762a1bSJed Brown */ 1859566063dSJacob Faibussowitsch PetscCall(MatFDColoringCreate(J, iscoloring, &fdcoloring)); 1869566063dSJacob Faibussowitsch PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void))FormFunction, &user)); 1879566063dSJacob Faibussowitsch PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 1889566063dSJacob Faibussowitsch PetscCall(MatFDColoringSetUp(J, iscoloring, fdcoloring)); 189c4762a1bSJed Brown /* 190c4762a1bSJed Brown Tell SNES to use the routine SNESComputeJacobianDefaultColor() 191c4762a1bSJed Brown to compute Jacobians. 192c4762a1bSJed Brown */ 1939566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring)); 1949566063dSJacob Faibussowitsch PetscCall(ISColoringDestroy(&iscoloring)); 19578e7fe0eSHong Zhang if (prunejacobian) PetscCall(SNESPruneJacobianColor(snes, J, J)); 196c4762a1bSJed Brown } 197c4762a1bSJed Brown /* 198c4762a1bSJed Brown Set Jacobian matrix data structure and default Jacobian evaluation 199c4762a1bSJed Brown routine. Whenever the nonlinear solver needs to compute the 200c4762a1bSJed Brown Jacobian matrix, it will call this routine. 201c4762a1bSJed Brown - Note that the final routine argument is the user-defined 202c4762a1bSJed Brown context that provides application-specific data for the 203c4762a1bSJed Brown Jacobian evaluation routine. 204c4762a1bSJed Brown - The user can override with: 205c4762a1bSJed Brown -snes_fd : default finite differencing approximation of Jacobian 206c4762a1bSJed Brown -snes_mf : matrix-free Newton-Krylov method with no preconditioning 207c4762a1bSJed Brown (unless user explicitly sets preconditioner) 208c4762a1bSJed Brown -snes_mf_operator : form preconditioning matrix as set by the user, 209c4762a1bSJed Brown but use matrix-free approx for Jacobian-vector 210c4762a1bSJed Brown products within Newton-Krylov method 211c4762a1bSJed Brown */ 212c4762a1bSJed Brown else if (!matrix_free) { 2139566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, (void *)&user)); 214c4762a1bSJed Brown } 215c4762a1bSJed Brown 216c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 217c4762a1bSJed Brown Customize nonlinear solver; set runtime options 218c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 219c4762a1bSJed Brown 220c4762a1bSJed Brown /* 221c4762a1bSJed Brown Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 222c4762a1bSJed Brown */ 2239566063dSJacob Faibussowitsch PetscCall(SNESSetFromOptions(snes)); 224c4762a1bSJed Brown 225c4762a1bSJed Brown /* 226c4762a1bSJed Brown Set array that saves the function norms. This array is intended 227c4762a1bSJed Brown when the user wants to save the convergence history for later use 228c4762a1bSJed Brown rather than just to view the function norms via -snes_monitor. 229c4762a1bSJed Brown */ 2309566063dSJacob Faibussowitsch PetscCall(SNESSetConvergenceHistory(snes, history, hist_its, 50, PETSC_TRUE)); 231c4762a1bSJed Brown 232c4762a1bSJed Brown /* 233c4762a1bSJed Brown Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the 234c4762a1bSJed Brown user provided test before the specialized test. The convergence context is just an array to 235c4762a1bSJed Brown test that it gets properly freed at the end 236c4762a1bSJed Brown */ 237c4762a1bSJed Brown if (use_convergence_test) { 2389566063dSJacob Faibussowitsch PetscCall(SNESGetKSP(snes, &ksp)); 2399566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(5, &testarray)); 2409566063dSJacob Faibussowitsch PetscCall(KSPSetConvergenceTest(ksp, ConvergenceTest, testarray, ConvergenceDestroy)); 241c4762a1bSJed Brown } 242c4762a1bSJed Brown 243c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244c4762a1bSJed Brown Evaluate initial guess; then solve nonlinear system 245c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 246c4762a1bSJed Brown /* 247c4762a1bSJed Brown Note: The user should initialize the vector, x, with the initial guess 248c4762a1bSJed Brown for the nonlinear solver prior to calling SNESSolve(). In particular, 249c4762a1bSJed Brown to employ an initial guess of zero, the user should explicitly set 250c4762a1bSJed Brown this vector to zero by calling VecSet(). 251c4762a1bSJed Brown */ 2529566063dSJacob Faibussowitsch PetscCall(FormInitialGuess(&user, x)); 2539566063dSJacob Faibussowitsch PetscCall(SNESSolve(snes, NULL, x)); 2549566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(snes, &its)); 25563a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its)); 256c4762a1bSJed Brown 257c4762a1bSJed Brown /* 258c4762a1bSJed Brown Print the convergence history. This is intended just to demonstrate 259c4762a1bSJed Brown use of the data attained via SNESSetConvergenceHistory(). 260c4762a1bSJed Brown */ 2619566063dSJacob Faibussowitsch PetscCall(PetscOptionsHasName(NULL, NULL, "-print_history", &flg)); 262c4762a1bSJed Brown if (flg) { 26348a46eb9SPierre Jolivet 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])); 264c4762a1bSJed Brown } 265c4762a1bSJed Brown 2663201ab8dSStefano Zampini /* Test NewtonTR API */ 2673201ab8dSStefano Zampini PetscCall(SNESNewtonTRSetTolerances(snes, 1.0, 2.0, 3.0)); 2683201ab8dSStefano Zampini PetscCall(SNESNewtonTRSetUpdateParameters(snes, 4.0, 5.0, 6.0, 7.0, 8.0)); 2693201ab8dSStefano Zampini PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg)); 2703201ab8dSStefano Zampini if (flg) { 2713201ab8dSStefano Zampini PetscReal tmp[5]; 2723201ab8dSStefano Zampini 2733201ab8dSStefano Zampini PetscCall(SNESNewtonTRGetTolerances(snes, &tmp[0], &tmp[1], &tmp[2])); 2743201ab8dSStefano Zampini PetscCheck(tmp[0] == 1.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2753201ab8dSStefano Zampini PetscCheck(tmp[1] == 2.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2763201ab8dSStefano Zampini PetscCheck(tmp[2] == 3.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2773201ab8dSStefano Zampini PetscCall(SNESNewtonTRGetUpdateParameters(snes, &tmp[0], &tmp[1], &tmp[2], &tmp[3], &tmp[4])); 2783201ab8dSStefano Zampini PetscCheck(tmp[0] == 4.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2793201ab8dSStefano Zampini PetscCheck(tmp[1] == 5.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2803201ab8dSStefano Zampini PetscCheck(tmp[2] == 6.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2813201ab8dSStefano Zampini PetscCheck(tmp[3] == 7.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2823201ab8dSStefano Zampini PetscCheck(tmp[4] == 8.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value"); 2833201ab8dSStefano Zampini } 2843201ab8dSStefano Zampini 285c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 286c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they 287c4762a1bSJed Brown are no longer needed. 288c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 289c4762a1bSJed Brown 29048a46eb9SPierre Jolivet if (!matrix_free) PetscCall(MatDestroy(&J)); 29148a46eb9SPierre Jolivet if (fd_coloring) PetscCall(MatFDColoringDestroy(&fdcoloring)); 2929566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 2939566063dSJacob Faibussowitsch PetscCall(VecDestroy(&r)); 2949566063dSJacob Faibussowitsch PetscCall(SNESDestroy(&snes)); 2959566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 296b122ec5aSJacob Faibussowitsch return 0; 297c4762a1bSJed Brown } 298f6dfbefdSBarry Smith 299c4762a1bSJed Brown /* 300c4762a1bSJed Brown FormInitialGuess - Forms initial approximation. 301c4762a1bSJed Brown 302c4762a1bSJed Brown Input Parameters: 303c4762a1bSJed Brown user - user-defined application context 304c4762a1bSJed Brown X - vector 305c4762a1bSJed Brown 306c4762a1bSJed Brown Output Parameter: 307c4762a1bSJed Brown X - vector 308c4762a1bSJed Brown */ 309d71ae5a4SJacob Faibussowitsch PetscErrorCode FormInitialGuess(AppCtx *user, Vec X) 310d71ae5a4SJacob Faibussowitsch { 311c4762a1bSJed Brown PetscInt i, j, row, mx, my; 312c4762a1bSJed Brown PetscReal lambda, temp1, temp, hx, hy; 313c4762a1bSJed Brown PetscScalar *x; 314c4762a1bSJed Brown 3153ba16761SJacob Faibussowitsch PetscFunctionBeginUser; 316c4762a1bSJed Brown mx = user->mx; 317c4762a1bSJed Brown my = user->my; 318c4762a1bSJed Brown lambda = user->param; 319c4762a1bSJed Brown 320c4762a1bSJed Brown hx = 1.0 / (PetscReal)(mx - 1); 321c4762a1bSJed Brown hy = 1.0 / (PetscReal)(my - 1); 322c4762a1bSJed Brown 323c4762a1bSJed Brown /* 324c4762a1bSJed Brown Get a pointer to vector data. 325c4762a1bSJed Brown - For default PETSc vectors, VecGetArray() returns a pointer to 326c4762a1bSJed Brown the data array. Otherwise, the routine is implementation dependent. 327c4762a1bSJed Brown - You MUST call VecRestoreArray() when you no longer need access to 328c4762a1bSJed Brown the array. 329c4762a1bSJed Brown */ 3309566063dSJacob Faibussowitsch PetscCall(VecGetArray(X, &x)); 331c4762a1bSJed Brown temp1 = lambda / (lambda + 1.0); 332c4762a1bSJed Brown for (j = 0; j < my; j++) { 333c4762a1bSJed Brown temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy; 334c4762a1bSJed Brown for (i = 0; i < mx; i++) { 335c4762a1bSJed Brown row = i + j * mx; 336c4762a1bSJed Brown if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 337c4762a1bSJed Brown x[row] = 0.0; 338c4762a1bSJed Brown continue; 339c4762a1bSJed Brown } 340c4762a1bSJed Brown x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp)); 341c4762a1bSJed Brown } 342c4762a1bSJed Brown } 343c4762a1bSJed Brown 344c4762a1bSJed Brown /* 345c4762a1bSJed Brown Restore vector 346c4762a1bSJed Brown */ 3479566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(X, &x)); 3483ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 349c4762a1bSJed Brown } 350f6dfbefdSBarry Smith 351c4762a1bSJed Brown /* 352c4762a1bSJed Brown FormFunction - Evaluates nonlinear function, F(x). 353c4762a1bSJed Brown 354c4762a1bSJed Brown Input Parameters: 355c4762a1bSJed Brown . snes - the SNES context 356c4762a1bSJed Brown . X - input vector 357c4762a1bSJed Brown . ptr - optional user-defined context, as set by SNESSetFunction() 358c4762a1bSJed Brown 359c4762a1bSJed Brown Output Parameter: 360c4762a1bSJed Brown . F - function vector 361c4762a1bSJed Brown */ 362d71ae5a4SJacob Faibussowitsch PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr) 363d71ae5a4SJacob Faibussowitsch { 364c4762a1bSJed Brown AppCtx *user = (AppCtx *)ptr; 365c4762a1bSJed Brown PetscInt i, j, row, mx, my; 366c4762a1bSJed Brown PetscReal two = 2.0, one = 1.0, lambda, hx, hy, hxdhy, hydhx; 367c4762a1bSJed Brown PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f; 368c4762a1bSJed Brown const PetscScalar *x; 369c4762a1bSJed Brown 3703ba16761SJacob Faibussowitsch PetscFunctionBeginUser; 371c4762a1bSJed Brown mx = user->mx; 372c4762a1bSJed Brown my = user->my; 373c4762a1bSJed Brown lambda = user->param; 374c4762a1bSJed Brown hx = one / (PetscReal)(mx - 1); 375c4762a1bSJed Brown hy = one / (PetscReal)(my - 1); 376c4762a1bSJed Brown sc = hx * hy; 377c4762a1bSJed Brown hxdhy = hx / hy; 378c4762a1bSJed Brown hydhx = hy / hx; 379c4762a1bSJed Brown 380c4762a1bSJed Brown /* 381c4762a1bSJed Brown Get pointers to vector data 382c4762a1bSJed Brown */ 3839566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(X, &x)); 3849566063dSJacob Faibussowitsch PetscCall(VecGetArray(F, &f)); 385c4762a1bSJed Brown 386c4762a1bSJed Brown /* 387c4762a1bSJed Brown Compute function 388c4762a1bSJed Brown */ 389c4762a1bSJed Brown for (j = 0; j < my; j++) { 390c4762a1bSJed Brown for (i = 0; i < mx; i++) { 391c4762a1bSJed Brown row = i + j * mx; 392c4762a1bSJed Brown if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 393c4762a1bSJed Brown f[row] = x[row]; 394c4762a1bSJed Brown continue; 395c4762a1bSJed Brown } 396c4762a1bSJed Brown u = x[row]; 397c4762a1bSJed Brown ub = x[row - mx]; 398c4762a1bSJed Brown ul = x[row - 1]; 399c4762a1bSJed Brown ut = x[row + mx]; 400c4762a1bSJed Brown ur = x[row + 1]; 401c4762a1bSJed Brown uxx = (-ur + two * u - ul) * hydhx; 402c4762a1bSJed Brown uyy = (-ut + two * u - ub) * hxdhy; 403c4762a1bSJed Brown f[row] = uxx + uyy - sc * lambda * PetscExpScalar(u); 404c4762a1bSJed Brown } 405c4762a1bSJed Brown } 406c4762a1bSJed Brown 407c4762a1bSJed Brown /* 408c4762a1bSJed Brown Restore vectors 409c4762a1bSJed Brown */ 4109566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(X, &x)); 4119566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(F, &f)); 4123ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 413c4762a1bSJed Brown } 414f6dfbefdSBarry Smith 415c4762a1bSJed Brown /* 416c4762a1bSJed Brown FormJacobian - Evaluates Jacobian matrix. 417c4762a1bSJed Brown 418c4762a1bSJed Brown Input Parameters: 419c4762a1bSJed Brown . snes - the SNES context 420c4762a1bSJed Brown . x - input vector 421c4762a1bSJed Brown . ptr - optional user-defined context, as set by SNESSetJacobian() 422c4762a1bSJed Brown 423c4762a1bSJed Brown Output Parameters: 424c4762a1bSJed Brown . A - Jacobian matrix 425c4762a1bSJed Brown . B - optionally different preconditioning matrix 426*0b4b7b1cSBarry Smith 427c4762a1bSJed Brown */ 428d71ae5a4SJacob Faibussowitsch PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr) 429d71ae5a4SJacob Faibussowitsch { 430da81f932SPierre Jolivet AppCtx *user = (AppCtx *)ptr; /* user-defined application context */ 431c4762a1bSJed Brown PetscInt i, j, row, mx, my, col[5]; 432c4762a1bSJed Brown PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc; 433c4762a1bSJed Brown const PetscScalar *x; 434c4762a1bSJed Brown PetscReal hx, hy, hxdhy, hydhx; 435c4762a1bSJed Brown 4363ba16761SJacob Faibussowitsch PetscFunctionBeginUser; 437c4762a1bSJed Brown mx = user->mx; 438c4762a1bSJed Brown my = user->my; 439c4762a1bSJed Brown lambda = user->param; 440c4762a1bSJed Brown hx = 1.0 / (PetscReal)(mx - 1); 441c4762a1bSJed Brown hy = 1.0 / (PetscReal)(my - 1); 442c4762a1bSJed Brown sc = hx * hy; 443c4762a1bSJed Brown hxdhy = hx / hy; 444c4762a1bSJed Brown hydhx = hy / hx; 445c4762a1bSJed Brown 446c4762a1bSJed Brown /* 447c4762a1bSJed Brown Get pointer to vector data 448c4762a1bSJed Brown */ 4499566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(X, &x)); 450c4762a1bSJed Brown 451c4762a1bSJed Brown /* 452c4762a1bSJed Brown Compute entries of the Jacobian 453c4762a1bSJed Brown */ 454c4762a1bSJed Brown for (j = 0; j < my; j++) { 455c4762a1bSJed Brown for (i = 0; i < mx; i++) { 456c4762a1bSJed Brown row = i + j * mx; 457c4762a1bSJed Brown if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) { 4589566063dSJacob Faibussowitsch PetscCall(MatSetValues(jac, 1, &row, 1, &row, &one, INSERT_VALUES)); 459c4762a1bSJed Brown continue; 460c4762a1bSJed Brown } 4619371c9d4SSatish Balay v[0] = -hxdhy; 4629371c9d4SSatish Balay col[0] = row - mx; 4639371c9d4SSatish Balay v[1] = -hydhx; 4649371c9d4SSatish Balay col[1] = row - 1; 4659371c9d4SSatish Balay v[2] = two * (hydhx + hxdhy) - sc * lambda * PetscExpScalar(x[row]); 4669371c9d4SSatish Balay col[2] = row; 4679371c9d4SSatish Balay v[3] = -hydhx; 4689371c9d4SSatish Balay col[3] = row + 1; 4699371c9d4SSatish Balay v[4] = -hxdhy; 4709371c9d4SSatish Balay col[4] = row + mx; 4719566063dSJacob Faibussowitsch PetscCall(MatSetValues(jac, 1, &row, 5, col, v, INSERT_VALUES)); 472c4762a1bSJed Brown } 473c4762a1bSJed Brown } 474c4762a1bSJed Brown 475c4762a1bSJed Brown /* 476c4762a1bSJed Brown Restore vector 477c4762a1bSJed Brown */ 4789566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(X, &x)); 479c4762a1bSJed Brown 480c4762a1bSJed Brown /* 481c4762a1bSJed Brown Assemble matrix 482c4762a1bSJed Brown */ 4839566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 4849566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 485c4762a1bSJed Brown 486c4762a1bSJed Brown if (jac != J) { 4879566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 4889566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 489c4762a1bSJed Brown } 4903ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 491c4762a1bSJed Brown } 492c4762a1bSJed Brown 493d71ae5a4SJacob Faibussowitsch PetscErrorCode ConvergenceTest(KSP ksp, PetscInt it, PetscReal nrm, KSPConvergedReason *reason, void *ctx) 494d71ae5a4SJacob Faibussowitsch { 4953ba16761SJacob Faibussowitsch PetscFunctionBeginUser; 496c4762a1bSJed Brown *reason = KSP_CONVERGED_ITERATING; 497c4762a1bSJed Brown if (it > 1) { 498c4762a1bSJed Brown *reason = KSP_CONVERGED_ITS; 4999566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL, "User provided convergence test returning after 2 iterations\n")); 500c4762a1bSJed Brown } 5013ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 502c4762a1bSJed Brown } 503c4762a1bSJed Brown 504d71ae5a4SJacob Faibussowitsch PetscErrorCode ConvergenceDestroy(void *ctx) 505d71ae5a4SJacob Faibussowitsch { 5063ba16761SJacob Faibussowitsch PetscFunctionBeginUser; 5079566063dSJacob Faibussowitsch PetscCall(PetscInfo(NULL, "User provided convergence destroy called\n")); 5089566063dSJacob Faibussowitsch PetscCall(PetscFree(ctx)); 5093ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 510c4762a1bSJed Brown } 511c4762a1bSJed Brown 512d71ae5a4SJacob Faibussowitsch PetscErrorCode postcheck(SNES snes, Vec x, Vec y, Vec w, PetscBool *changed_y, PetscBool *changed_w, void *ctx) 513d71ae5a4SJacob Faibussowitsch { 514c4762a1bSJed Brown PetscReal norm; 515c4762a1bSJed Brown Vec tmp; 516c4762a1bSJed Brown 5173ba16761SJacob Faibussowitsch PetscFunctionBeginUser; 5189566063dSJacob Faibussowitsch PetscCall(VecDuplicate(x, &tmp)); 5199566063dSJacob Faibussowitsch PetscCall(VecWAXPY(tmp, -1.0, x, w)); 5209566063dSJacob Faibussowitsch PetscCall(VecNorm(tmp, NORM_2, &norm)); 5219566063dSJacob Faibussowitsch PetscCall(VecDestroy(&tmp)); 5229566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of search step %g\n", (double)norm)); 5233ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 524c4762a1bSJed Brown } 525c4762a1bSJed Brown 526c4762a1bSJed Brown /*TEST 527c4762a1bSJed Brown 528c4762a1bSJed Brown build: 529c4762a1bSJed Brown requires: !single 530c4762a1bSJed Brown 531c4762a1bSJed Brown test: 532c4762a1bSJed Brown args: -ksp_gmres_cgs_refinement_type refine_always 533c4762a1bSJed Brown 534c4762a1bSJed Brown test: 535c4762a1bSJed Brown suffix: 2 536c4762a1bSJed Brown args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always 537c4762a1bSJed Brown 538c4762a1bSJed Brown test: 539c4762a1bSJed Brown suffix: 2a 540c4762a1bSJed Brown filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault" 541c4762a1bSJed Brown args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info 542dfd57a17SPierre Jolivet requires: defined(PETSC_USE_INFO) 543c4762a1bSJed Brown 544c4762a1bSJed Brown test: 545c4762a1bSJed Brown suffix: 2b 546c4762a1bSJed Brown filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test" 547c4762a1bSJed Brown args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info 548dfd57a17SPierre Jolivet requires: defined(PETSC_USE_INFO) 549c4762a1bSJed Brown 550c4762a1bSJed Brown test: 55124fb275aSStefano Zampini suffix: 2c 55224fb275aSStefano Zampini 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 55324fb275aSStefano Zampini 55424fb275aSStefano Zampini test: 555c4762a1bSJed Brown suffix: 3 556c4762a1bSJed Brown args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always 557c4762a1bSJed Brown 558c4762a1bSJed Brown test: 559c4762a1bSJed Brown suffix: 4 560c4762a1bSJed Brown args: -pc -par 6.807 -snes_monitor -snes_converged_reason 561c4762a1bSJed Brown 56278e7fe0eSHong Zhang test: 56378e7fe0eSHong Zhang suffix: 5 56478e7fe0eSHong Zhang args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always -prune_jacobian 56578e7fe0eSHong Zhang output_file: output/ex1_3.out 566ccb5f961SBarry Smith 567ccb5f961SBarry Smith test: 568ccb5f961SBarry Smith suffix: 6 569ccb5f961SBarry Smith args: -snes_monitor draw:image:testfile -viewer_view 570ccb5f961SBarry Smith 5719be84c52SStefano Zampini test: 5729be84c52SStefano Zampini suffix: python 5739be84c52SStefano Zampini requires: petsc4py 5749be84c52SStefano Zampini args: -python -snes_type python -snes_python_type ex1.py:MySNES -snes_view -null_appctx {{0 1}separate output} 5759be84c52SStefano Zampini localrunfiles: ex1.py 5769be84c52SStefano Zampini 577c4762a1bSJed Brown TEST*/ 578