PETSc - University of Saskatchewan
Contents
1. Include petscts h to use the PETSc timestepping routines Note that this file automatically includes petsc h and other lower level PETSc include files include petscts h Create an application context to contain data needed by the application provided call back routines FormJacobian and FormFunction typedef struct AppCtx 86 User defined routines extern PetscErrorCode RHSFunction TS PetscReal Vec Vec void undef __FUNCT__ define __FUNCT__ main int main int argc char argv TS ts timestepping context Vec y solution residual vectors PetscErrorCode ierr PetscReal end_time dt PetscReal ftime PetscInt num_time_steps AppCtx user user defined work context PetscScalar xinitial_condition PetscInitialize amp argc amp argv PETSC_NULL help Set up the timestep can be an option from command line dt 0 5 end_time 1 0 ierr Petsc0ptionsGetReal PETSC_NULL dt amp dt PETSC_NULL CHKERRQ ierr num_time_steps round end_time dt Create vector to hold the solution ierr VecCreateSeq PETSC_COMM_SELF 1 amp y CHKERRQ ierr Create timestepper context ierr TSCreate PETSC_COMM_WORLD amp ts CHKERRQ ierr ierr TSSetProblemType ts TS_NONLINEAR CHKERRQ ierr 87 Set initial condition ierr VecGetArray y amp initial_condition CHKER2. Monitoring routines for SNES can be cancelled at runtime via the snes_monitor_cancel option 83 Solving nonlinear equations The solution vector and function vector from an SNES context can be accessed via SNESGetSolution SNES snes Vec x SNESGetFunction SNES snes Vec r void ctx int func SNES Vec Vec void Not only is this useful to see the solution and function values yourself e g during the iteration or after convergence but this would be how such information would be passed to a custom convergence test that required them 84 Solving IVPs PETSc can be used to solve initial value problems IVPs in ordinary differential equations ODEs and differential algebraic equations DAEs Consider the following example R dt y y 0 1 First we set up a makefile to build the executable Makefile to build ode_example c include PETSC_DIR conf base ode_example ode_example o chkopts CLINKER o ode_example ode_example o PETSC_TS_LIB RM ode_example o 85 Solving ODEs DAEs The following code now solves the above IVP using Euler s method static char help Solves the simple ODE dy dt y y 0 1 Concepts solving ordinary differential equations Processors 1 This code demonstrates how one may solve an ordinary differential equation using the built in ODE solvers
3. Linear functions solved with an implicit solver are invoked via TSSetMatrices TS ts Mat A PetscErrorCode frhs TS PetscReal Mat Mat MatStructure void Mat B PetscErrorCode flhs TS PetscReal Mat Mat MatStructure void MatStructure flag void ctx where the solver is assumed to take the form The functions frhs and flhs are used to form A and B at each time step if they are time dependent If not the user should pass in PETSC_NULL 102 Solving ODEs DAEs The variable ctx allows users to pass in an application context that is passed to the frhs or flhs function whenever they are called as the final argument The user must provide the matrix A If B I the user should pass in PETSC_NULL If the right hand side is provided only as a linear function the user must construct a MATSHELL matrix shell to provide a user defined matrix type We note that this is the same interface as that for SNESSet Jacobian 103 Solving ODEs DAEs The right hand side for nonlinear problems or linear problems solved using explicit methods is set up via TSSetRHSFunction TS ts Vec R PetscErrorCode f TS double Vec Vec void void fP where the vector R is an optional location to store the result The arguments to the function are the time step context the current time the input for the function the output for the function and the optional user provided context variable P
4. Suppose we have an initial guess x for approximating a desired solution x to F x 0 The classical version of Newton s method is formally defined by the following iteration ee e Jg x F x 1 le eee where x is the nth approximation to the solution of F x O and Jp x is the Jacobian matrix evaluated at x 47 Solving nonlinear equations In practice this is implemented via solution of the following system of linear equations at each iteration Ip x a F x followed by the update XOH x n 4 gr where d is called the Newton direction 48 Solving nonlinear equations In order for Newton s method to be practical we need e a termination criterion that ensures a sufficiently accurate approximate solution e an efficient linear system solver to compute the Newton direction d e a globalization strategy to aid the convergence of the iteration to a solution of F x 0 for any x 49 Solving nonlinear equations In practice a choice must be made as to how accurately to solve for d e When x is far from the solution d may not be useful so it may not make sense to solve for it accurately If you do that is called oversolving e When x is close to the solution d may be useful so it may make sense to solve for it accurately If you do not that is called undersolving A balance between these undesirable extremes is achieved through the use
5. ksp_type lt type gt pc_type lt type gt ksp_monitor ksp_rtol lt rtol gt These options will override those specified above as long as KSPSetFrom0ptions is called _after_ any other customization routines ierr KSPSetFrom0ptions ksp CHKERRQ ierr if monzeroguess PetscScalar p 5 ierr VecSet x p CHKERRQ ierr ierr KSPSetInitialGuessNonzero ksp PETSC_TRUE CHKERRQ ierr x Solve the linear system Solve linear system ierr KSPSolve ksp b x CHKERRQ ierr 25 View solver info we could instead use the option ksp_view to print this info to the screen at the conclusion of KSPSolve ierr KSPView ksp PETSC_VIEWER_STDOUT_WORLD CHKERRQ ierr Check the error ierr VecAXPY x neg_one u CHKERRQ ierr ierr VecNorm x NORM_2 amp norm CHKERRQ ierr ierr KSPGetIterationNumber ksp amp its CHKERRQ ierr if norm gt tol ierr PetscPrintf PETSC_COMM_WORLD Norm of error G Iterations D n norm its CHKERRQ ierr Free work space All PETSc objects should be destroyed when they are no longer needed ierr VecDestroy amp x CHKERRQ ierr ierr VecDestroy amp u CHKERRQ ierr ierr VecDestroy amp b CHKERRQ ierr ierr MatDestroy amp A CHKERRQ ierr ierr KSPDestroy amp ksp CHKERRQ ierr Always call PetscFinalize before exiting a p
6. 104 Solving ODEs DAEs For nonlinear problems the user must also provide the approximate Jacobian matrix of F t u and a function to compute it at each Newton iteration This is done via TSSetRHSJacobian TS ts Mat A Mat P PetscErrorCode fjac TS double Vec Mat Mat MatStructure void void fP where the arguments for the function fjac are the time step context the current time the location where the Jacobian is to be computed the Jacobian matrix an alternative approximate Jacobian matrix used as a preconditioner and the optional user provided context passed in as fP The user must provide the Jacobian as a matrix thus if a matrix free approach is used the user must create a MATSHELL matrix We again note the similarity to SNESSetJacobian 105 Solving ODEs DAEs To solve a DAE instead of TSSetRHSFunction and TSSetRHSJacobian one uses TSSetIFunction TS ts Vec R PetscErrorCode f TS PetscReal Vec Vec Vec void void funP where the vector R is an optional location to store the residual The arguments to are the time step context current time input state u input time derivative u_t and the optional user provided context funP 106 Solving ODEs DAEs Unless one is using a matrix free method without preconditioning the user must also provide an approximate Jacobian matrix of G u F t u w au because the integrator internally sets w au whe
7. int indices PetscScalar values INSERT_VALUES Vector components can be set to the same scalar value via VecSet Vec x PetscScalar value Note the use of the PET Sc variable type PetscScalar The PetscScalar type is defined to be double in C C or double precision in Fortran for versions of PETSc that have not been compiled for use with complex numbers The PetscScalar data type enables code to be re used when the PETSc libraries have been compiled for use with complex numbers Values can be added to existing components of a vector by using the ADD_VALUES insert mode 31 Solving Ax b in serial Assigning individual values to components of a vector is more complicated two step process but with the advantage that the end result is more efficient underlying parallel code The first call is to VecSetValues Vec x int n int indices PetscScalar values INSERT VALUES where n gives the number of components to be set in this call The integer array indices contains the global component indices and values is the array of values to be set Any process can set any components of the vector PETSc insures that they are automatically stored in the correct location 32 Solving Ax b in serial However once all of the values have been inserted with VecSetValues one must call VecAssemblyBegin Vec x followed by VecAssemblyEnd Vec x so components can be transferred to the appropriate process as
8. 0 xx 1i 3 0 ierr VecRestoreArray x amp xx CHKERRQ ierr 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 ierr SNESSolve snes PETSC_NULL x CHKERRQ ierr ierr SNESGetIterationNumber snes kits CHKERRQ ierr if flg 63 Vec f ierr VecView x PETSC_VIEWER_STDOUT_WORLD CHKERRQ ierr SNESGetFunction snes amp f 0 0 CHKERRQ ierr VecView r PETSC_VIEWER_STDOUT_WORLD CHKERRQ ierr ierr ierr ierr PetscPrintf PETSC_COMM_WORLD number of SNES iterations D n its CHKERRQ Free work space All PETSc objects should be destroyed when they are no longer needed ierr VecDestroy amp x CHKERRQ ierr ierr VecDestroy amp r CHKERRQ ierr ierr MatDestroy amp J CHKERRQ ierr ierr SNESDestroy amp snes CHKERRQ ierr if size gt 1 ierr VecDestroy kuser xloc CHKERRQ ierr ierr VecDestroy user rloc CHKERRQ ierr ierr VecScatterDestroy amp user scatter CHKERRQ ierr ierr PetscFinalize return 0 undef __FUNCT__ define __FUNCT__ FormFunctioni FormFunction1 Evaluates nonlinear function F x Input Parameters snes the SNES context x input
9. h viewers petscpc h preconditioners Note The corresponding uniprocessor example is exl c include lt petscksp h gt undef __FUNCT__ define __FUNCT__ main 39 int main int argc char args Vec x b u approx solution RHS exact solution Mat A linear system matrix KSP ksp linear solver context PC pc preconditioner context PetscReal norm tol 1 e 11 norm of solution error PetscErrorCode ierr PetscInt i n 10 col 3 its rstart rend nlocal PetscScalar neg_one 1 0 one 1 0 value 3 PetscInitialize amp argc amp args char 0 help ierr PetscOptionsGetInt PETSC_NULL n amp n PETSC_NULL CHKERRQ ierr Compute the matrix and right hand side vector that define the linear system Ax b Create vectors Note that we form 1 vector from scratch and then duplicate as needed For this simple case let PETSc decide how many elements of the vector are stored on each processor The second argument to VecSetSizes below causes PETSc to decide ierr VecCreate PETSC_COMM_WORLD amp x CHKERRQ ierr ierr VecSetSizes x PETSC_DECIDE n CHKERRQ ierr ierr VecSetFrom0ptions x CHKERRQ ierr ierr VecDuplicate x amp b CHKERRQ ierr ierr VecDuplicate x amp u CHKERRQ ierr Identify the starting and ending mesh points on each processor for the interior part of the mesh We let PETSc decide above ierr VecGetOwnership
10. linear and nonlinear via SNESSetFromOptions snes This function provides an interface to the PETSc options database so that choices can be made at runtime for things such as the e particular nonlinear solver e various parameters and customized routines e g specialized line search variants e convergence tolerance e monitoring routines e linear solver options in the KSP and PC modules 71 Solving nonlinear equations The problem is then solved via a call to SNESSolve SNES snes Vec b Vec x where x will contain the solution vector x should be initialized to the initial guess for the nonlinear solver before calling SNESSolve For example to employ an initial guess of zero not always a good idea x should be explicitly set to zero by calling VecSet Finally after all the nonlinear system solving is done the SNES context should be destroyed via SNESDestroy SNES snes 12 Solving nonlinear equations The user specifies the system of nonlinear equations F x 0 to be solved via SNESSetFunction SNES snes Vec f PetscErrorCode FormFunction SNES snes Vec x Vec f void ctx void ctx The argument ctx is an optional user defined context that can store any private application specific data required by the function evaluation routine Use PETSC_NULL if such information is not needed Some kind of routine to approximate Jp x is typically specified with SNESSetJacobian SNES snes Mat
11. undef __FUNCT__ define __FUNCT__ FormJacobian2 PetscErrorCode FormJacobian2 SNES snes Vec x Mat jac Mat B MatStructure flag void PetscScalar xx A 4 PetscErrorCode ierr PetscInt idx 2 0 1 Get pointer to vector data 68 ierr VecGetArray x amp xx CHKERRQ ierr Compute Jacobian entries and insert into matrix Since this is such a small problem we set all entries for the matrix at once A O 3 0 PetscCosScalar 3 0 xx 0 1 0 A 1 0 0 A 2 0 0 A 3 1 0 ierr MatSetValues B 2 idx 2 idx A INSERT_VALUES CHKERRQ ierr flag SAME_NONZERO_PATTERN Restore vector ierr VecRestoreArray x amp xx CHKERRQ ierr Assemble matrix ierr MatAssemblyBegin B MAT_FINAL_ASSEMBLY CHKERRQ ierr ierr MatAssemblyEnd B MAT_FINAL_ASSEMBLY CHKERRQ ierr if jac B ierr MatAssemblyBegin jac MAT_FINAL_ASSEMBLY CHKERRQ ierr ierr MatAssemblyEnd jac MAT_FINAL_ASSEMBLY CHKERRQ ierr return 0 69 Solving nonlinear equations Notes e An SNES solver is created via a call of the form SNESCreate MPI Comm comm SNES snes e The routines to evaluate F x and Jp x are set e The nonlinear solution method is chosen via SNESSetType SNES snes SNESType method or from command line via snes_type lt method gt 70 Solving nonlinear equations An application code can take complete control of the solvers both
12. when trying to solve Problem 1 with np gt 1 it simply returns an incorrect answer 59 Solving nonlinear equations static char help Newton s method for a two variable system sequential n n T Concepts SNES basic example Tx 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 lt petscsnes h gt typedef struct Vec xloc rloc local solution residual vectors VecScatter scatter AppCtx User defined routines extern PetscErrorCode FormJacobiani SNES Vec Mat Mat MatStructure void extern PetscErrorCode FormFunction1 SNES Vec Vec void extern PetscErrorCode FormJacobian2 SNES Vec Mat Mat MatStructure void extern PetscErrorCode FormFunction2 SNES Vec Vec void undef __FUNCT__ define __FUNCT__ main 60 int main int argc char argv SNES snes nonlinear solver context KSP ksp linear solver context PC pc preconditioner context Vec X t solution residual vectors Mat J Jacobian matrix PetscErrorCode ierr PetscInt its PetscMPIInt size rank PetscScalar pfive 5 xx PetscBool flg AppCtx user user defined work context IS isglobal islo
13. x 1 _ The term ine search refers to a procedure for finding the step length A The accepted strategy is to first try A 1 If this fails to satisfy the criterion that makes x acceptable then we backtrack in a systematic fashion along d 15 Solving nonlinear equations Experience has shown the importance of taking the full Newton step whenever possible Failure to do so leads to sub optimal convergence rates near the solution Although no criterion will always be optimal it seems to be good sense to require E x FY lt F If we compute Jp x along with F x then we will have 4 pieces of information function values and derivatives at each of the two end points A 0 and A 1 of the line search This means we can fit a cubic interpolating polynomial to these points and find its minimizer analytically This value of A is taken as the value to be used when determining x 76 Solving nonlinear equations Alternative line search routines can be set via SNESSetLineSearch SNES snes PetscErrorCode 1s SNES Vec Vec Vec Vec double double double void lsctx PETSc also provides other line search methods SNESSearchQuadraticLine SNESLineSearchNo and SNESLineSearchNoNorms which can be set via snes_ls cubic quadratic basic basicnonorms where cubic is the default method just discussed quadratic is the corresponding backtracking line search method based on
14. 1 0 value 1 2 0 value 2 1 0 for i rstart i lt rend i col 0 i 1 col 1 i col 2 it1 ierr MatSetValues A 1 amp i1 3 col value INSERT_VALUES CHKERRQ ierr Assemble the matrix ierr MatAssemblyBegin A MAT_FINAL_ASSEMBLY CHKERRQ ierr ierr MatAssemblyEnd A MAT_FINAL_ASSEMBLY CHKERRQ ierr Set exact solution then compute right hand side vector ierr VecSet u one CHKERRQ ierr ierr MatMult A u b CHKERRQ ierr x Create the linear solver and set various options Create linear solver context ierr KSPCreate PETSC_COMM_WORLD amp ksp CHKERRQ ierr Set operators Here the matrix that defines the linear system also serves as the preconditioning matrix ierr KSPSetOperators ksp A A DIFFERENT_NONZERO_PATTERN CHKERRQ ierr Set linear solver defaults for this problem optional By extracting the KSP and PC contexts from the KSP context we can then directly call any KSP and PC routines to set various options The following four statements are optional all of these parameters could alternatively be specified at runtime via 42 KSPSetFromOptions ierr KSPGetPC ksp amp pc CHKERRQ ierr ierr PCSetType pc PCJACOBI CHKERRQ ierr ierr KSPSetTolerances ksp 1 e 7 PETSC_DEFAULT PETSC_DEFAULT PETSC_DEFAULT CHKERRQ ierr Set r
15. A Mat B PetscErrorCode FormJacobian SNES snes Vec x Mat A Mat B MatStructure flag void ctx void ctx where x is the current iterate A is the Jacobian matrix B is the preconditioner usually A flag gives information about the preconditioner structure identical to the options for the flag of KSPSetOperators 13 Solving nonlinear equations ctx Is an optional user defined Jacobian context for application specific data All SNES solvers are data structure neutral so all PETSc matrix formats including matrix free methods can be used It Is common practice to assemble the Jacobian into the preconditioner B This has no effect in the common case that A and B are identical But it allows us to check A by passing the snes_mf_operator flag PETSc then constructs a finite difference approximation to Jp in A and B remains as the preconditioner Even if B is incorrect the iteration will likely converge using A thus this procedure can be used to verify the correctness of B 74 Solving nonlinear equations The globalization strategies built in to PETSc are based on line search and trust region methods The method SNESLS which can also be invoked via the flag snes_type ls provides a line search for Newton s method The default strategy is called cubic backtracking The idea of a line search is natural given a descent direction d take a step in that direction that yields an acceptable
16. GetInt PETSC_NULL n amp n PETSC_NULL CHKERRQ ierr ierr PetscOptionsGetBool PETSC_NULL nonzero_guess knonzeroguess PETSC_NULL CHKERRQ ierr Compute the matrix and right hand side vector that define the linear system Ax b Create vectors Note that we form 1 vector from scratch and then duplicate as needed ierr VecCreate PETSC_COMM_WORLD amp x CHKERRQ ierr ierr PetscObjectSetName PetscObject x Solution CHKERRQ ierr ierr VecSetSizes x PETSC_DECIDE n CHKERRQ ierr ierr VecSetFrom0ptions x CHKERRQ ierr ierr VecDuplicate x amp b CHKERRQ ierr ierr VecDuplicate x amp u CHKERRQ ierr Create matrix When using MatCreate the matrix format can be specified at runtime Performance tuning note For problems of substantial size preallocation of matrix memory is crucial for attaining good performance See the matrix chapter of the users manual for details ierr MatCreate PETSC_COMM_WORLD amp A CHKERRQ ierr 23 ierr MatSetSizes A PETSC_DECIDE PETSC_DECIDE n n CHKERRQ ierr ierr MatSetFrom0ptions A CHKERRQ ierr ierr MatSetUp A CHKERRQ ierr Assemble matrix value 0 1 0 value 1 2 0 value 2 1 0 for i 1 i lt n 1 i col 0O i 1 col 1 i col 2 it1 ierr MatSetValues A 1 amp i1 3 col value INSERT_VALUES CHKERRQ ierr i n 1 colf O n 2 col i n 1 ierr MatSet
17. Nicolson method and the solvers from the SUNDIALS package CVODE and IDA SUNDIALS stands for the SUite of Nonlinear DIfferential A Lgebraic equation Solvers It was developed at the Lawrence Livermore National Laboratory where it continues to be supported developed and maintained CVODE Is a powerful general purpose IVP solver that includes methods for non stiff and stiff systems In the case of stiff systems the implicit solvers can use direct full or banded linear system solvers for the Newton iteration as well as iterative preconditioned Krylov subspace methods At present CVODE contains three Krylov subspace methods GMRES Bi CGStab and TFQMR 91 Solving ODEs DAEs In general GMRES is considered to be the most robust solver and hence the one to try first An advantage of CVODE as in PETSc is that it is easy to try the other Krylov space solvers if convergence issues are encountered with GMRES Historically there used to be a parallel version of CVODE called PVODE but it has been incorporated Into CVODE The family of solvers known as SUNDIALS actually consists of CVODE for IVPs in ODEs KINSOL for nonlinear algebraic equations IDA for IVPs in DAEs and the extensions CVODES and IDAS that can perform sensitivity analyses both forward and adjoint 92 Solving ODEs DAEs The TS library in PETSc provides a framework for solving ODEs and DAEs in parallel The main motivation is to solve such p
18. PETSc Portable Extensible Toolkit for Scientific Computing Raymond J Spiteri Lecture Notes for CMPT 898 Numerical Software University of Saskatchewan Objectives Introduction to PETSc Installation and use Linear system solvers Nonlinear system solvers ODE DAE solvers Introduction to PETSc PETSc petsi stands for the Portable Extensible Toolkit for Scientific Computing e suite of libraries data structures objects methods functions subroutines produced and maintained by Argonne National Laboratory for the solution of scientific applications modelled by partial differential equations PDEs e intended for use in large scale HPC applications e includes a large suite of parallel linear and nonlinear algebraic equation solvers as well as ODE solvers e sets up many of the mechanisms required for parallel application codes including parallel matrix and vector assembly e integrates BLAS Basic Linear Algebra Subprograms and LAPACK Linear Algebra PACKage as the back end for basic vector and matrix operations Introduction to PETSc e can be used in Windows Mac and Linux e written in C e can be used within C C Fortran Python and MATLAB applications e interfaces with many other software packages making it easy to try new methods and solvers e used in many applications including heart simulation carbon sequestration ground water flow medical imaging and quantum computing Se
19. RQ ierr initial_condition 0 1 0 ierr VecRestoreArray y amp initial_condition CHKERRQ ierr ierr TSSetSolution ts y CHKERRQ ierr Provide the call back for the nonlinear function we are evaluating dy dt f t y Thus whenever the timestepping routines need the function they will call this routine Note the final argument is the application context used by the call back functions ierr TSSetRHSFunction ts PETSC_NULL RHSFunction amp user CHKERRQ ierr This indicates that we are using Euler s method ierr TSSetType ts TSEULER CHKERRQ ierr Set the initial time and the initial timestep given above ierr TSSetInitialTimeStep ts 0 0 dt CHKERRQ ierr Set a maximum number of timesteps and final simulation time ierr TSSetDuration ts num_time_steps end_time Set any additional options from the options database This includes all options for the nonlinear and linear solvers used internally the the timestepping routines ierr TSSetFromOptions ts CHKERRQ ierr TSSetUp ts CHKERRQ ierr ierr 88 Perform the solve This is where the timestepping takes place ierr TSSolve ts y amp ftime CHKERRQ ierr View information about the time stepping method and the solution at the end time TSView ts PETSC_VIEWER_STDOUT_SELF VecView y PETSC_VIEWER_STDOUT_SELF printf nThis is ftime f n ftime Free the da
20. Range x amp rstart amp rend CHKERRQ ierr VecGetLocalSize x amp nlocal CHKERRQ ierr ierr 40 Create matrix When using MatCreate the matrix format can be specified at runtime Performance tuning note For problems of substantial size preallocation of matrix memory is crucial for attaining good performance See the matrix chapter of the users manual for details We pass in nlocal as the local size of the matrix to force it to have the same parallel layout as the vector created above ierr MatCreate PETSC_COMM_WORLD amp A CHKERRQ ierr ierr MatSetSizes A nlocal nlocal n n CHKERRQ ierr ierr MatSetFrom0ptions A CHKERRQ ierr ierr MatSetUp A CHKERRQ ierr if if Assemble matrix The linear system is distributed across the processors by chunks of contiguous rows which correspond to contiguous sections of the mesh on which the problem is discretized For matrix assembly each processor contributes entries for the part that it owns locally rstart rstart 1 i 0 col 0 0 col 1 1 value 0 2 0 value 1 1 0 ierr MatSetValues A 1 amp i 2 col value INSERT_VALUES CHKERRQ ierr rend n rend n 1 i n 1 col 0 n 2 col 1 n 1 value 0 1 0 value 1i 2 0 ierr MatSetValues A 1 amp i1 2 col value INSERT_VALUES CHKERRQ ierr 41 Set entries corresponding to the mesh interior value O
21. TS ts TSType type where the supported types are TSEULER TSRK Runge Kutta TSBEULER TSCN Crank Nicolson TSTHETA TSGL generalized linear TSPSEUDO and TSSUNDIALS if the Sundials package is installed or the command line option ts_type euler rk beuler cn theta gl pseudo sundials Normally it is best to use the TSSetFrom0ptions command and then set the TS type from the options database rather than by using this routine 99 Solving ODEs DAEs The initial time and time step are set via TSSetInitialTimeStep TS ts double time double dt The time step can be changed via TSSetTimeStep TS ts double dt The current time step can be determined via TSGetTimeStep TS ts double dt The maximum number of steps or time to use whichever occurs first is set via TSSetDuration TS ts int maxsteps double maxtime A specific number of steps can be executed via TSSolve TS ts Vec U PetscReal ftime A single step can be taken via TSStep TS ts 100 Solving ODEs DAEs The TSSolve call implicitly sets up the timestep context but it can be done explicitly with TSSetUp TS ts The context can be viewed via TSView TS ts PetscViewer viewer and destroyed via TSDestroy TS ts Initial conditions are set via TSSetSolution TS ts Vec initialsolution 101 Solving ODEs DAEs The right hand side function for ODEs is set according to whether it is linear and an implicit solver is being used
22. Values A 1 amp i 2 col value INSERT_VALUES CHKERRQ ierr i 0 col 0 0 col 1 1 value 0 2 0 value 1i 1 0 ierr MatSetValues A 1 amp i 2 col value INSERT_VALUES CHKERRQ ierr ierr MatAssemblyBegin A MAT_FINAL_ASSEMBLY CHKERRQ ierr ierr MatAssemblyEnd A MAT_FINAL_ASSEMBLY CHKERRQ ierr Set exact solution then compute right hand side vector ierr VecSet u one CHKERRQ ierr ierr MatMult A u b CHKERRQ ierr x Create the linear solver and set various options Create linear solver context ierr KSPCreate PETSC_COMM_WORLD amp ksp CHKERRQ ierr Set operators Here the matrix that defines the linear system also serves as the preconditioning matrix ierr KSPSetOperators ksp A A DIFFERENT_NONZERO_PATTERN CHKERRQ ierr 24 Set linear solver defaults for this problem optional By extracting the KSP and PC contexts from the KSP context we can then directly call any KSP and PC routines to set various options The following four statements are optional all of these parameters could alternatively be specified at runtime via KSPSetFrom0ptions ierr KSPGetPC ksp amp pc CHKERRQ ierr PCSetType pc PCJACOBI CHKERRQ ierr ierr KSPSetTolerances ksp 1 e 5 PETSC_DEFAULT PETSC_DEFAULT PETSC_DEFAULT CHKERRQ ierr ierr Set runtime options e g
23. _VALUES SCATTER_REVERSE CHKERRQ ierr ierr VecScatterEnd scatter floc f INSERT_VALUES SCATTER_REVERSE CHKERRQ ierr else Get pointers to vector data 65 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 ierr VecGetArrayRead x amp xx CHKERRQ ierr VecGetArray f amp ff CHKERRQ ierr ierr Compute function ff 0 xx 0 xx 0 xx 0 xx 1 3 ff 1 xx O xx 1 xx 1 xx 1i 6 0 0 gt Restore vectors ierr VecRestoreArrayRead x amp xx CHKERRQ ierr ierr VecRestoreArray f amp ff CHKERRQ ierr return 0 undef __FUNCT__ define __FUNCT__ FormJacobiani FormJacobianl Evaluates Jacobian matrix Input Parameters snes the SNES context x input vector dummy optional user defined context not used here Output Parameters jac Jacobian matrix B optionally different preconditioning matrix flag flag indicating matrix structure PetscErrorCode FormJacobiani SNES snes Vec x Mat jac Mat B MatStructure flag void PetscScalar xx A 4 PetscErrorCode ierr 66 PetscInt idx 2 0 1 Get pointer to vector data ierr VecGetArray x amp xx CHKERRQ ierr Compute Jacobian entries and inse
24. a quadratic interpolant basic is really only a template for a classical Newton iteration and basicnonorms is a classical Newton iteration T7 Solving nonlinear equations It is not recommended to use either basic or basicnonorms options in general lt may make sense to use the basicnonorms option when you know you want to perform a fixed number of full Newton iterations This might be the case when trying to achieve an efficient approximation to a solution starting from a highly accurate initial guess Examples include predictor corrector methods and fully implicit Runge Kutta methods for initial value problems and error estimators and deferred corrections for boundary value problems The line search routines involve several parameters with defaults that are often reasonable in practice they can be overridden via snes_ls_alpha lt alpha gt snes_ls_maxstep lt max gt and snes_ls_steptol lt tol gt 78 Solving nonlinear equations Line search methods attempt to compute an acceptable step length in a given search direction The underlying assumptions were that the direction would be d and A 1 would be the first trial step In trust region methods we relax the assumption that the step be in the direction d Indeed there is more information available from Jp x and trust region methods try to use it Trust region methods use a pre specified step length For example it may not be so unrea
25. braries as well as MPI provides summary and diagnostic information if certain runtime options are chosen e g log_summary ierr PetscFinalize return 0 44 Solving Ax b in parallel It is significant to realize that serial and parallel code can look the same the behaviour changes based on the definition of the communicator In practice there are some things to keep in mind e he creation routines are collective over all processes in the communicator thus all processes in the communicator must call the creation routine e Any sequence of collective routines must be called in the same order on each process e It is possible to have finer control over matrix vector distribution in the parallel version 45 Solving nonlinear equations It is also possible to specify and solve systems of nonlinear algebraic equations in PETSc Systems of nonlinear algebraic equations take the form where F R is the nonlinear residual function x R is the vector of unknowns and O is a vector of zeros Notes e There may a unique solution multiple solutions or no solution e Analysis of existence and uniqueness of solutions is generally difficult but especially for large m e Finding a suitable initial guess can also be difficult 46 Solving nonlinear equations The only robust methods for solving systems of nonlinear algebraic equations in practice are based on Newton s method
26. c petsc current docs manual pdf PETSc Library Hierarchy Level of Abstraction Application Codes a TS SNES Time Stepping Nonlinear Equations Solvers PC KSP Preconditioners Krylov Subspace Methods Matrices Vectors Index Sets BLAS MPI Figure 1 Hierarchy of PETSc libraries PETSc libraries Parallel Numerical Components of PETSc Nonlinear Solvers i Time Steppers Newton based Methods IMEX Pseudo_Ti seudo Time Other General Lingar Runge Kutt Line Search Trust Region Stepping Krylov Subspace Methods GMRES CG CGS Bi CG Stab TFQMR Richardson Chebychev Other Preconditioners Additive Block LU Schwarz Tacobi Jacobi ILU ICC sequential only Other Matrices Compressed Block Compressed Symmetric Sparse Row Sparse Row Block Compressed Row Dense Other AIJ BAID SBAIJ Index Sets Vectors Indices Block Indices Stride Other Figure 2 PETSc libraries PET Sc components PETSc provides components such as e IS Index sets for vectors matrices e Vec Vector operations parallel scatter and gather e Mat Parallel sparse matrix operations including four parallel matrix data structures e DM Data management between vectors matrices and
27. cal PetscInitialize amp argc amp argv char 0 help ierr MPI_Comm_size PETSC_COMM_WORLD amp size CHKERRQ ierr ierr MPI_Comm_rank PETSC_COMM_WORLD amp rank CHKERRQ ierr ierr SNESCreate PETSC_COMM_WORLD amp snes CHKERRQ ierr Create vectors for solution and nonlinear function ierr VecCreate PETSC_COMM_WORLD amp x CHKERRQ ierr ierr VecSetSizes x PETSC_DECIDE 2 CHKERRQ ierr ierr VecSetFrom0ptions x CHKERRQ ierr ierr VecDuplicate x amp r CHKERRQ ierr if size gt 1 ierr VecCreateSeq PETSC_COMM_SELF 2 amp user xloc CHKERRQ ierr ierr VecDuplicate user xloc user rloc CHKERRQ ierr 61 Create the scatter between the global x and local xloc ierr ISCreateStride MPI_COMM_SELF 2 0 1 amp islocal CHKERRQ ierr ierr ISCreateStride MPI_COMM_SELF 2 0 1 amp isglobal CHKERRQ ierr ierr VecScatterCreate x isglobal user xloc islocal user scatter CHKERRQ ierr ierr ISDestroy amp isglobal CHKERRQ ierr ierr ISDestroy amp islocal CHKERRQ ierr Create Jacobian matrix data structure ierr MatCreate PETSC_COMM_WORLD amp J CHKERRQ ierr ierr MatSetSizes J PETSC_DECIDE PETSC_DECIDE 2 2 CHKERRQ ierr ierr MatSetFrom0ptions J CHKERRQ ierr ierr MatSetUp J CHKERRQ ierr ierr PetscOptionsHasName PETSC_NULL hard amp flg CHKERRQ ierr if flg Set function evaluation routine a
28. d with the command MatSetValues Mat A int m int im int n int in PetscScalar values INSERT_VALUES After all elements have been inserted into the matrix it must be processed with the pair of commands MatAssemblyBegin Mat A MAT_FINAL_ASSEMBLY MatAssemblyEnd Mat A MAT_FINAL_ASSEMBLY This is necessary because values may be buffered so this pair of calls ensures all relevant buffers are flushed Again other code can be placed between these two calls to perform computation while data are in transit And again calls to MatSetValues with the INSERT_VALUES and ADD_VALUES options cannot be mixed without calling the assembly routines Intermediate assembly calls should use the cheaper routine MAT_FLUSH_ASSEMBLY MAT_FINAL_ASSEMBLY is required only before a matrix is used 36 Solving Ax b in serial After having created the system Ax b we can now solve it using the KSP package This can be achieved by the following sequence of commands KSPCreate MPI Comm comm KSP ksp KSPSetOperators KSP ksp Mat A Mat PrecA MatStructure flag KSPSetFromOptions KSP ksp KSPSolve KSP ksp Vec b Vec x KSPDestroy KSP ksp First the KSP context is created then the system operators the coefficient matrix A and optionally a preconditioning matrix are set Various options for customized solution are then assigned the linear system is solved and finally the KSP context is destroyed The command KSPS
29. distributed meshes e PC Sequential and parallel preconditioners e KSP Parallel Krylov subspace iterative methods e SNES Nonlinear algebraic equation solvers e TS Time steppers for ordinary differential equations and differential algebraic equations Parallelism in PETSc PETSc uses the Message Passing Interface MPI for its parallelization The calls made to MPI in PETSc are at a high level and mostly hidden from the user typically users should not have to call MPI directly themselves although they are free to do so Nonetheless familiarity with the basic concepts of message passing and distributed memory computing is essential when using PET Sc In these notes we assume users have this basic level of familiarity 10 Installation PETSc has comprehensive instructions for installation on Linux and Windows with a variety of configurations PETSc can download and install its own dependencies making the installation process easy for new users One way to install PETSc is to download it from the main website http www mcs anl gov petsc download index html Assuming PETSc is at version xxxx we unpack the downloaded tar file via the following command tar xzvf petsc xxxx tar gz Another way is to download it using via hg clone https bitbucket org petsc petsc 3 3 petsc 3 3 hg clone https bitbucket org petsc buildsystem 3 3 petsc 3 3 config BuildSystem 11 Installation Next assuming a PETSc installa
30. e http www mcs anl gov petsc for details Software accessible through PE TSc PETSc can interface with other software applications including e ADIC ADIFOR for computing sparse Jacobian matrices by automatic differentiation e The Unsymmetric MultiFrontal Sparse LU factorization PACKage UMFPACK for the direct solution of unsymmetric sparse linear systems e The MuUltifrontal Massively Parallel sparse direct Solver MUMPS for large sparse linear systems e SuperLU a library for the direct solution of large Sparse nonsymmetric linear systems e Prometheus a scalable unstructured finite element solver e Trilinos ML a multigrid preconditioner package Software currently using PETSc PETSc is used by many software packages including e TAO Toolkit for Advanced Optimization e SLEPc Scalable Library for Eigenvalue Problems e deal ll an adaptive finite element solver e Chaste Cancer Heart and Soft Tissue Environment Features of PETSc Features of PETSc include intended to be easy to learn and use for beginners extensive examples and tutorials e control over numerical methods for advanced users e parallel vectors and matrices e parallel linear system solvers based on Krylov subspace methods including Conjugate Gradient e parallel preconditioners e parallel nonlinear solvers based on Newton s method e parallel ODE solvers e extensive documentation complete user s manual http www mcs anl gov pets
31. etFrom0ptions enables customization of the linear solution method at runtime via the options database e g choice of iterative method preconditioner convergence tolerance etc 37 Error Checking All PETSc routines return an integer ierr indicating whether an error has occurred during the call The PETSc macro CHKERRQ ierr checks ierr and calls the PETSc error handler upon error detection It is strongly advised that CHKERRQ ierr be used in all calls to PETSc to enable a complete error trace The debug version of PETSc does a great deal of checking for memory corruption e g writing outside of array bounds etc The macros CHKMEMQ can be called anywhere in the code to check the current memory status for corruption By strategically placing these macros in your code you can usually pinpoint any problematic segment of code 38 Solving Ax b in parallel We now go through the same example of solving Ax b but this time in parallel The program we are considering is ex23 c Program usage mpiexec ex23 help all PETSc options static char help Solves a tridiagonal linear system n n T Concepts KSP basic parallel example Processors n T Include petscksp h so that we can use KSP 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
32. har argv char file char xhelp where argc argv are the usual command line arguments or in Fortran call PetscInitialize character file integer ierr PetscInitialize calls MPI_Init if MPI has not been previously initialized If MPI needs to be initialized directly or is initialized by some other library MPI_Init can be called before PetscInitialize By default PetscInitialize sets the PETSc world communicator PETSC_COMM_WORLD to MPI_COMM_WORLD 19 Writing PETSc Programs e All PETSc routines return an integer indicating whether an error has occurred during the call The error code is nonzero if an error has been detected otherwise it is zero It is a good idea to check this error code after every call to a PETSc routine e All PETSc programs should call PetscFinalize before finishing The C and Fortran syntax respectively are PetscFinalize call PetscFinalize ierr If MPI was initialized from PetscInitialize PetscFinalize calls MPI_Finalize otherwise MPI_Finalize must be called directly 20 Solving Ax b in serial We now examine ex1 c in detail This example solves the linear system Ax b of dimension 10 where A is a tridiagonal matrix with elements 1 2 1 and b has been constructed so that the exact solution x is a vector of ones i e x 1 1 1 1 1 1 1 1 1 1 SSS SSS Set SoS Cie 6 6 wie COCCCORN FO SSO SO Ose Nw Hoo oo oo Ne So O 5 S
33. nd vector ierr SNESSetFunction snes r FormFunctioni user CHKERRQ ierr Set Jacobian matrix data structure and Jacobian evaluation routine ierr SNESSetJacobian snes J J FormJacobian1 PETSC_NULL CHKERRQ ierr else if size 1 SETERRQ PETSC_COMM_SELF 1 This case is a uniprocessor example only ierr SNESSetFunction snes r FormFunction2 PETSC_NULL CHKERRQ ierr ierr SNESSetJacobian snes J J FormJacobian2 PETSC_NULL CHKERRQ ierr 62 Set linear solver defaults for this problem By extracting the KSP KSP and PC contexts from the SNES context we can then directly call any KSP KSP and PC routines to set various options ierr SNESGetKSP snes amp ksp CHKERRQ ierr ierr KSPGetPC ksp amp pc CHKERRQ ierr ierr PCSetType pc PCNONE CHKERRQ ierr ierr KSPSetTolerances ksp 1 e 4 PETSC_DEFAULT PETSC_DEFAULT 20 CHKERRQ ierr Set SNES KSP KSP PC runtime options e g snes_view snes_monitor ksp_type lt ksp gt pc_type lt pc gt These options will override those specified above as long as SNESSetFrom0ptions is called _after_ any other customization routines ierr SNESSetFrom0ptions snes CHKERRQ ierr x Evaluate initial guess then solve nonlinear system if flg 4 ierr VecSet x pfive CHKERRQ ierr else ierr VecGetArray x amp xx CHKERRQ ierr xx 0 2
34. needed Communication and calculation can overlap in between i e your code can perform other actions between these two calls while any data are in transition Note that PETSc does not allow simultaneous insertion and addition because the operations could be performed in either order thus resulting in non deterministic behaviour Insertions and additions must be separated by calls to VecAssemblyBegin and VecAssemblyEnd 33 Solving Ax b in serial One can examine a vector with the command VecView Vec x PetscViewer v To print the vector to the screen one can use the viewer PETSC_VIEWER_STDOUT_WORLD which ensures that parallel vectors are printed correctly to stdout To display the vector in an X window one can use the default X windows viewer PETSC_VIEWER_DRAW_WORLD or one can create a viewer with the routine PetscViewerDrawOpenXx 34 Solving Ax b in serial e New parallel or sequential matrices A with M global rows and N global columns can be created via MatCreate MPI_Comm comm Mat A MatSetSizes Mat A int m int n int M int N where the matrix format can be specified at runtime Alternatively the number of local rows and columns can be specified using m and n Generally the matrix type is specified with e g MatSetType Mat A MATATJ This stores the matrix entries in the compressed Sparse row storage format 35 Solving Ax b in serial Values can then be inserte
35. of a dynamic variable called the forcing term 7 50 Solving nonlinear equations Compute the Newton direction d that satisfies the inexact Newton condition Le x lt nE x where Lp x F x Ip x d is called the local linear model We assume that L x closely approximates F x t The inexact Newton condition requires that Lp x be smaller than F x by a factor of n Solving nonlinear equations e n controls the level of accuracy of d In the case of direct methods 7 0 the inexact Newton condition holds with equality In the case of iterative methods nh determines the number of inner iterations to be performed e n influences the rate of convergence and performance of Newton s method If L x approximates F x poorly choosing n too small imposes too much accuracy on d thus leading to oversolving Choosing 7 too large may fail to approximate a sufficiently accurate d thus leading to undersolving Oversolving or undersolving may result in little or no reduction in the residual 52 Solving nonlinear equations A practical Newton algorithm looks something like this Input initial iterate x residual function F x absolute tolerance T and relative tolerance T Output the approximate solution x 1 2 3 4 5 6 T 8 9 X x9 while termination criterion is not met do Cho
36. oo Sw eo oo Ss CO Ce iH 6 O O 045 OrPNFODOCCSO PFPNFODOCOCCCSO NFODOCOOCCCSO on TN a as js B 21 Solving Ax b in serial Program usage mpiexec ex1 help all PETSc options static char help Solves a tridiagonal linear system with KSP n n T Concepts KSP solving a system of linear equations Processors 1 Tx Include petscksp h so that we can use KSP 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 Note The corresponding parallel example is ex23 c include lt petscksp h gt undef __FUNCT__ define FUNCT__ main int main int argc char args Vec x b u approx solution RHS exact solution Mat A linear system matrix KSP ksp linear solver context PC pc preconditioner context PetscReal norm tol 1 e 14 norm of solution error PetscErrorCode ierr 22 PetscInt i n 10 col 3 its PetscMPIInt size PetscScalar neg_one 1 0 one 1 0 value 3 PetscBool nonzeroguess PETSC_FALSE PetscInitialize amp argc amp args char 0 help ierr MPI_Comm_size PETSC_COMM_WORLD amp size CHKERRQ ierr if size 1 SETERRQ PETSC_COMM_WORLD 1 This is a uniprocessor example only ierr PetscOptions
37. ose a forcing term 7 Find d so that F x Jp x d lt 7 F x If d cannot be found terminate with failure Find a step length A x x Ad end while return x We discuss the strategies for choosing 7 and A shortly 53 Solving nonlinear equations The suite of routines for solving systems of nonlinear algebraic equations in PETSc is called SNES In the simplest usage of SNES the user merely provides a C C or Fortran routine to evaluate F x By default the corresponding Jacobian matrix can be approximated with finite differences For improved efficiency and robustness the user can provide a routine to compute the Jacobian analytically or call an automatic differentiation package The user has available a well trusted forcing term strategy due to Eisenstat and Walker line search or trust region methods as globalization strategies all the KSP methods and preconditioners for solving linear systems and a modified Newton method strategy for updating the Jacobian For most options including for the termination criterion users are free to define their own functions 54 Solving nonlinear equations We now go through the solution procedure for the following two nonlinear problems each consisting of two variables and two equations Problem 1 r xy 3 y zty 6 with exact solution els yy S and initial guess 1 1 Seo yy I 9 55 Solving nonlinear equations Problem 2 with exac
38. ote that this example is for a single process the value passed to the number of processes flag np is 1 14 Installation For the corresponding parallel example we can try make ex23 mpiexec np 2 ex23 We can pass any value to np up to the number of processes allowed on the computer being used In this directory of tutorials and examples each example has C or Fortran code and an HTML file that describes the code in detail For complete installation instructions including instructions for installation under Windows see the installation page available from PETSc http www mcs anl gov petsc documentation installation html 15 Running PETSc All PETSc programs support the h or help option as well as the v or version option The following is a list of some of the other more useful options a complete list is available via the h option e log_summary summarize program performance e fp_trap stop on floating point exception e malloc_dump list free memory at end of run e malloc_debug enable memory tracing e info print detailed output as program runs e options_file lt filename gt read options from file It is also possible start all processes in a debugger or start the debugger only when encountering an error 16 Writing PETSc Programs It is easy to run the examples that come with PETSc but to create and run applications of your own additional work is required A practical starting stra
39. quation differentiating the boundary conditions is usually OK 96 Solving ODEs DAEs Different spatial discretizations such as those obtained from the finite element method would yield systems of the form Mu f t u where M is invertible Often M has properties such as symmetry or positive definiteness that make it more amenable to inversion It may be convenient to treat such a system as a differential algebraic equation F t u 0 If M M t or M M y there may be no choice In general if cr is non singular for all t u the system u can be written as an ODE otherwise it is a DAE 97 Solving ODEs DAEs The solution of steady state problems 0 F u can be solved in some cases by pseudo transient continuation by solving F u Generally this is done with the backward Euler method because we would like to take large steps and reach a suitably large final time such that u converges to its steady state It is important to note that transient solution values i e at for intermediate times that have yet to converge are essentially meaningless 98 Solving ODEs DAEs The TS library is used in a manner similar to other libraries in PETSc A TS object is created with the command int TSCreate MPI Comm comm TSProblemType problemtype TS ts The TSProblemType is set to TS_LINEAR if u A t u or TS_NONLINEAR if F t u The solution method is set via TSSetType
40. re the positive shift a and vector w depend on the integration method step size and past states The function to evaluate G u Fu aFy is set via TSSetIJacobian TS ts Mat A Mat B PetscErrorCode fjac TS PetscReal Vec Vec PetscReal Mat Mat MatStructure void void jacP where fjac takes the time step context current time input state u input derivative u_t shift a matrix A preconditioning matrix B flag describing structure of preconditioning matrix see KSPSetOperators and the optional user provided context jacP 107 Summary PETSc scientific library used to solve applications modelled by partial differential equations Large suite of sequential and parallel nonlinear and linear algebraic equation solvers Used in many applications including heart simulation and carbon sequestration One of the go to libraries for solving large sparse linear systems MPI used for parallelization High level use of MPI for users Easy for beginners fine grained control for experts 108
41. roblems as resulting from the method of lines applied to time dependent PDEs It is also possible to solve time independent stationary steady state PDEs through the use of the method of lines combined with pseudo timestepping For example suppose a temperature distribution u x t satisfies the heat equation w Ups OO t gt O subject to the boundary and initial data ula S Ta UO S TTo U SE T 93 Solving ODEs DAEs This PDE can be discretized by the method of lines on a uniform spatial mesh with m 1 subintervals using centred differences for the spatial derivative to yield the initial value DAE 0 UO Ss Ta m e a Ui 1 Ui Ax y 0 Um 1 Lp i 1 2 m subject to u 0 To 2 t 0 1 m 1 94 Solving ODEs DAEs In this case one could take advantage of the explicit knowledge we have of the Dirichlet boundary conditions to write the problem as an IVP for an ODE ug 2u Ta y l U4 2u T U Nt M A m i E Tp z 2Um Um 1 mo Az yj subject to u 0 Totes A Pereak 95 Solving ODEs DAEs Sometimes people differentiate the BCs to obtain o 0 Wi 2U Ui 5 M asl E 6 2 Um 1 0 subject to u 0 Toa 2 0 1 m 1 This is generally not recommended for hyperbolic problems because noise generally propagates in from the boundary polluting the solution In the case of a parabolic problem such as the heat e
42. rogram This routine finalizes the PETSc libraries as well as MPI provides summary and diagnostic information if certain runtime options are chosen e g log_summary ierr PetscFinalize return 0 26 Solving Ax b in serial Notes e The instruction include petscksp h includes the header file for the Krylov subspace iterative linear solvers A PETSc program must specify a header file that corresponds to the highest level PETSc objects needed within the program all of the necessary lower level header files are automatically included from the higher level call For example petscksp h includes petscmat h for matrices petscvec h for vectors and petscsys h for base PETSc routines The PETSc header files are located in the directory PETSC_DIR include 27 Solving Ax b in serial e Input control data can be set at run time using the options database In this example the command PetscOptionsGet Int PETSC_NULL n amp n amp flg checks whether the user has provided a command line option to set the size of the problem n If so the variable n is set accordingly otherwise n remains unchanged at its default value of 10 For example mpiexec np 1 ex1 n 100 would run the corresponding problem but with dimension 100 28 Solving Ax b in serial e Parallel or sequential vectors x of global dimension M can be created via VecCreate MPI_Comm comm Vec x VecSetSize
43. rt into matrix Since this is such a small problem we set all entries for the matrix at once ACO 2 0 xx 0 xx 1 A 1 xx 0 A 2 xx 1 A 3 xx 0 2 0 xx 1 ierr MatSetValues B 2 idx 2 idx A INSERT_VALUES CHKERRQ ierr flag SAME_NONZERO_PATTERN Restore vector ierr VecRestoreArray x amp xx CHKERRQ ierr Assemble matrix ierr MatAssemblyBegin B MAT_FINAL_ASSEMBLY CHKERRQ ierr ierr MatAssemblyEnd B MAT_FINAL_ASSEMBLY CHKERRQ ierr if jac B ierr MatAssemblyBegin jac MAT_FINAL_ASSEMBLY CHKERRQ ierr ierr MatAssemblyEnd jac MAT_FINAL_ASSEMBLY CHKERRQ ierr return 0 undef __FUNCT__ define __FUNCT__ FormFunction2 Pe tscErrorCode FormFunction2 SNES snes Vec x Vec f void dummy 67 PetscErrorCode ierr PetscScalar xx ff Get pointers 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 ierr VecGetArray x amp xx CHKERRQ ierr ierr VecGetArray f amp ff CHKERRQ ierr Compute function f 0 PetscSinScalar 3 0 xx 0 xx 0 ff 1 xx 1 Restore vectors ierr VecRestoreArray x amp xx CHKERRQ ierr VecRestoreArray f amp ff CHKERRQ ierr return 0 ierr
44. s Vec x int m int M where comm denotes the MPI communicator and m is the local size which may be PETSC_DECIDE As such the vectors created as described are empty The type of storage for a vector may be set by calling either VecSetType or VecSetFromOptions Vectors of the same type can be formed via VecDuplicate Vec old Vec new or VecDuplicateVecs Vec old int n Vec new Sequential or distributed vectors can also be created directly using VecCreateSeq and VecCreateMPI 29 Solving Ax b in serial These constructs are useful because they allow one to write library code that does not depend on the particular format of the vectors being used Instead PETSc automatically and correctly creates work vectors based on the existing vector As usual when a vector or array of vectors is no longer needed it should be destroyed via VecDestroy Vec x or VecDestroyVecs PetscInt n Vec vecs Vectors that use an array provided by the user rather than having PETSc internally allocate the array space can be created via VecCreateSeqWithArray PETSC COMM SELF int bs int n PetscScalar array Vec V and VecCreateMPIWithArray MPI Comm comm int bs int n int N PetscScalar array Vec vv where n must be given i e not PETSC_DECIDE and there must be enough space n sizeof PetscScalar 30 Solving Ax b in serial Vector components can be set inserted via VecSetValues Vec x int n
45. sonable to begin with a step length on the order of the previous one The general idea is define a region in which the local linear model approximates the actual model well and adjust its size as the iteration proceeds Several effective heuristics are available for this 79 Solving nonlinear equations The trust region method used by SNES is invoked via SNESTR or the snes_type tr flag The strategy itself is taken from the MINPACK project There are several parameters that control the size of the trust region during the solution process For example the user can control the initial trust region radius A in A Ao Fol e by setting Ao via the flag snes_tr_delta0 lt delta0 gt 80 Solving nonlinear equations Convergence of the SNES solvers can be tested in various ways including by a user defined test Each convergence test involves several parameters with default values that should work well in practice One way to check convergence is to test whether x0 xD lt gtol Another way is to test E x lt atol x rtol Often a combination of the two is recommended These parameters are set via SNESSetTolerances SNES snes double atol double rtol double stol int its int fcts 81 Solving nonlinear equations We note this also sets the maximum numbers of nonlinear iterations its and function evaluations fcts The command line arguments for these settings are snes_atol lt a
46. t solution and initial guess 56 Solving nonlinear equations The nonlinear solver directory also contains many examples and tutorials We consider exi_snes c a program that solves each of the two nonlinear equations described above To solve Problem 1 by default we execute cd PETSC_DIR src snes examples tutorials make exi_snes mpiexec np 1 exi_snes The output is simply number of SNES iterations 6 57 Solving nonlinear equations To make it more informative we can add negate the flag in line 147 to view the complete solution if flg Vec f ierr VecView x PETSC_VIEWER_STDOUT_WORLD CHKERRQ ierr ierr SNESGetFunction snes amp f 0 0 CHKERRQ ierr ierr VecView r PETSC_VIEWER_STDOUT_WORLD CHKERRQ ierr j The default output is then Vector Object 1 MPI processes type seq 1 2 Vector Object 1 MPI processes type seq 4 44089e 16 O number of SNES iterations 6 58 Solving nonlinear equations We solve Problem 2 by including the flag hard when the executable file is run make exi_snes mpiexec np 1 exi_snes hard This gives the output Vector Object 1 MPI processes type seq 2 58086e 13 3 87421e 16 Vector Object 1 MPI processes type seq 1 03234e 12 3 87421e 16 number of SNES iterations 6 Note that this example does not work with np gt 1 A peculiarity of the code however is that it only reports an error when trying to solve Problem 2
47. ta structures constructed above ierr VecDestroy amp y CHKERRQ ierr ierr TSDestroy amp ts CHKERRQ ierr ierr PetscFinalize CHKERRQ ierr return 0 undef __FUNCT__ define __FUNCT__ RHSFunction PetscErrorCode RHSFunction TS ts PetscReal t Vec X Vec F void ptr x PetscErrorCode ierr PetscScalar y f ierr VecGetArray X amp y CHKERRQ ierr ierr VecGetArray F amp f CHKERRQ ierr f 0 y 0 dy dt y ierr VecRestoreArray X amp y CHKERRQ ierr VecRestoreArray F amp f CHKERRQ ierr return 0 ierr 89 Command line options It is possible to specify non default methods used to solve a problem via command line options when the executable is launched e g the preconditioner and the KSP solver It is also possible to specify that an external package such as MUMPS be used as the linear system solver Additional PETSc information about the run can also be displayed For example after running the makefile to build ex23 c the executable can be run with the ilu preconditioner cg solver and additional output including viewing the matrix A mpiexec np 2 ex23 ksp_type cg pc_type ilu log_summary mat_view The ability to change the methods used to solve a linear system makes it convenient test a variety of numerical methods on a given problem 90 Solving ODEs DAEs Other ODE solvers are available including the backward Euler method the Crank
48. tegy is to find an example in the many examples provided that most closely matches the problem you wish to solve After ensuring this example can be built and run successfully copy the contents into a new directory then use the programs and makefile as a starting point from which to being solving your problem An important part of this process is the construction of the makefile A makefile is a file containing a set of instructions required to build an executable from source code An excellent introduction to makefiles as well as other basic computing skills is available at http software carpentry org 17 Writing PETSc Programs Consider again the first example ex1 The makefile to build this example looks like Makefile to build exi include PETSC_DIR conf base ex1 exl o chkopts CLINKER o ex1 ex1 o PETSC_KSP_LIB RM ext o Because this example is written in C the C linker CLINKER is used This is an example of an internal environment variable known to and used by PETSc Another environment variable PETSC_KSP_LIB gives the linker the information required to build against the PET Sc libraries Finally recall that the PETSC_DIR environment variable was set when PETSc was installed and stores the path to the main PETSc directory 18 Writing PETSc Programs The general form of a PETSc program includes the following structure e Programs in C begin with PetscInitialize int argc c
49. tion on a Linux machine we set two important environment variables e PETSC_DIR gives the full path to the PETSc home install directory e PETSC_ARCH specifies the architecture particular set of compiler options and machine type In a bash script on a linux machine these variables can be set via cd petsc xxxx export PETSC_DIR pwd export PETSC_ARCH linux gnu 12 Installation The next step is to configure PETSc so that it knows where libraries are installed and what additional libraries are required the installation For example config configure py download f blas lapack 1 download openmpi yes with cc gcc with fc gfortran This example downloads and installs the MPI library openmpi independent of any existing installation This is the recommended approach assuming sufficient disk space etc However PETSc can be made to work with existing MPI installations 13 Installation Now PETSc must be built and tested This is done by executing the following command make all test It is also easy to try running one of the many PETSc examples provided making sure the environment variables PETSC_DIR and PETSC_ARCH are set cd PETSC_DIR src ksp ksp examples tutorials make exl mpiexec np 1 ex1 If you run into trouble using mpiexec it may be that the path to mpiexec has not been set From the bash shell the following command may help export PATH PATH PETSC_DIR PETSC_ARCH bin N
50. tol gt snes_rtol lt rtol gt snes_stol lt stol gt snes_max_it lt its gt and snes_max_funcs lt fcts gt The settings for the tolerances can be obtained via SNESGetTolerances Convergence tests for trust region methods have a setting for the minimum allowable trust region radius This parameter can be set via snes_trtol lt trtol gt or using SNESSetTrustRegionTolerance SNES snes double trtol Customized convergence tests can be defined via SNESSetConvergenceTest SNES snes PetscErrorCode test SNES snes int it double xnorm double gnorm double f SNESConvergedReason reason void cctx void cctx PetscErrorCode destroy void cctx 82 Solving nonlinear equations By default the SNES solvers do not display information about the iterations Monitoring can be established via SNESMonitorSet SNES snes PetscErrorCode mon SNES int its double norm void mctx void mctx PetscErrorCode monitordestroy void where mon is a user defined monitoring routine its and mctx are the iteration number and an optional user defined context for private data for the monitor routine and norm is the function norm SNESMonitorSet is called once after every outer SNES iteration hence any application specific computations can be done after each iteration The option snes_monitor activates the SNES monitor routine SNESMonitorDefault snes_monitor_draw draws a simple graph of F x
51. untime options e g ksp_type lt type gt pc_type lt type gt ksp_monitor ksp_rtol lt rtol gt These options will override those specified above as long as KSPSetFrom0ptions is called _after_ any other customization routines ierr KSPSetFrom0ptions ksp CHKERRQ ierr x Solve the linear system Solve linear system ierr KSPSolve ksp b x CHKERRQ ierr View solver info we could instead use the option ksp_view to print this info to the screen at the conclusion of KSPSolve ierr KSPView ksp PETSC_VIEWER_STDOUT_WORLD CHKERRQ ierr x Check solution and clean up Check the error ierr VecAXPY x neg_one u CHKERRQ ierr ierr VecNorm x NORM_2 amp norm CHKERRQ ierr 43 ierr KSPGetIterationNumber ksp amp its CHKERRQ ierr if norm gt tol ierr PetscPrintf PETSC_COMM_WORLD Norm of error G Iterations D n norm its Free work space All PETSc objects should be destroyed when they are no longer needed ierr VecDestroy amp x CHKERRQ ierr ierr VecDestroy amp u CHKERRQ ierr ierr VecDestroy amp b CHKERRQ ierr ierr MatDestroy amp A CHKERRQ ierr ierr KSPDestroy amp ksp CHKERRQ ierr Always call PetscFinalize before exiting a program This routine finalizes the PETSc li
52. vector ctx optional user defined context Output Parameter f function vector PetscErrorCode FormFunction1 SNES snes Vec x Vec f void ctx 64 PetscErrorCode ierr const PetscScalar xx PetscScalar f Ff AppCtx xuser AppCtx ctx Vec xloc user gt xloc floc user gt rloc VecScatter scatter user gt scatter MPI_Comm comm PetscMPIInt size rank PetscInt rstart rend ierr PetscObjectGetComm PetscObject snes amp comm CHKERRQ ierr ierr MPI_Comm_size comm amp size CHKERRQ ierr ierr MPI_Comm_rank comm amp rank CHKERRQ ierr if size gt 1 This is a ridiculous case for testing intermidiate steps from sequential code development to parallel implementation 1 scatter x into a sequetial vector 2 each process evaluates all values of floc 3 scatter floc back to the parallel f ierr VecScatterBegin scatter x xloc INSERT_VALUES SCATTER_FORWARD CHKERRQ ierr ierr VecScatterEnd scatter x xloc INSERT_VALUES SCATTER_FORWARD CHKERRQ ierr ierr VecGetOwnershipRange f amp rstart amp rend CHKERRQ ierr ierr VecGetArrayRead xloc amp xx CHKERRQ ierr ierr VecGetArray floc amp ff CHKERRQ ierr f 0 xx 0 xx 0 xx 0 xx 1 3 0 f 1 xx 0 xx 1 xx 1 xx 1 6 0 ierr VecRestoreArray floc amp ff CHKERRQ ierr ierr VecRestoreArrayRead xloc amp xx CHKERRQ ierr ierr VecScatterBegin scatter floc f INSERT
Download Pdf Manuals
Related Search