1c4762a1bSJed Brown! 2c4762a1bSJed Brown! 3c4762a1bSJed Brown! Description: Uses the Newton method to solve a two-variable system. 4c4762a1bSJed Brown! 5c4762a1bSJed Brown!!/*T 6c4762a1bSJed Brown! Concepts: SNES^basic uniprocessor example 7c4762a1bSJed Brown! Processors: 1 8c4762a1bSJed Brown!T*/ 9c4762a1bSJed Brown 10c4762a1bSJed Brown 11c4762a1bSJed Brown! 12c4762a1bSJed Brown! ----------------------------------------------------------------------- 13c4762a1bSJed Brown 14c4762a1bSJed Brown program main 15c4762a1bSJed Brown#include <petsc/finclude/petsc.h> 16c4762a1bSJed Brown use petsc 17c4762a1bSJed Brown implicit none 18c4762a1bSJed Brown 19c4762a1bSJed Brown 20c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 21c4762a1bSJed Brown! Variable declarations 22c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 23c4762a1bSJed Brown! 24c4762a1bSJed Brown! Variables: 25c4762a1bSJed Brown! snes - nonlinear solver 26c4762a1bSJed Brown! ksp - linear solver 27c4762a1bSJed Brown! pc - preconditioner context 28c4762a1bSJed Brown! ksp - Krylov subspace method context 29c4762a1bSJed Brown! x, r - solution, residual vectors 30c4762a1bSJed Brown! J - Jacobian matrix 31c4762a1bSJed Brown! its - iterations for convergence 32c4762a1bSJed Brown! 33c4762a1bSJed Brown SNES snes 34c4762a1bSJed Brown PC pc 35c4762a1bSJed Brown KSP ksp 36c4762a1bSJed Brown Vec x,r 37c4762a1bSJed Brown Mat J 38c4762a1bSJed Brown SNESLineSearch linesearch 39c4762a1bSJed Brown PetscErrorCode ierr 40c4762a1bSJed Brown PetscInt its,i2,i20 41c4762a1bSJed Brown PetscMPIInt size,rank 42c4762a1bSJed Brown PetscScalar pfive 43c4762a1bSJed Brown PetscReal tol 44c4762a1bSJed Brown PetscBool setls 45c4762a1bSJed Brown PetscViewer viewer 46c4762a1bSJed Brown double precision threshold,oldthreshold 47c4762a1bSJed Brown 48c4762a1bSJed Brown! Note: Any user-defined Fortran routines (such as FormJacobian) 49c4762a1bSJed Brown! MUST be declared as external. 50c4762a1bSJed Brown 51c4762a1bSJed Brown external FormFunction, FormJacobian, MyLineSearch 52c4762a1bSJed Brown 53c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 54c4762a1bSJed Brown! Macro definitions 55c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 56c4762a1bSJed Brown! 57c4762a1bSJed Brown! Macros to make clearer the process of setting values in vectors and 58c4762a1bSJed Brown! getting values from vectors. These vectors are used in the routines 59c4762a1bSJed Brown! FormFunction() and FormJacobian(). 60c4762a1bSJed Brown! - The element lx_a(ib) is element ib in the vector x 61c4762a1bSJed Brown! 62c4762a1bSJed Brown#define lx_a(ib) lx_v(lx_i + (ib)) 63c4762a1bSJed Brown#define lf_a(ib) lf_v(lf_i + (ib)) 64c4762a1bSJed Brown! 65c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 66c4762a1bSJed Brown! Beginning of program 67c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 68c4762a1bSJed Brown 69c4762a1bSJed Brown call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 70c4762a1bSJed Brown if (ierr .ne. 0) then 71c4762a1bSJed Brown print*,'Unable to initialize PETSc' 72c4762a1bSJed Brown stop 73c4762a1bSJed Brown endif 74*d21b9a37SPierre Jolivet call PetscLogNestedBegin(ierr);CHKERRA(ierr) 75c4762a1bSJed Brown threshold = 1.0 76c4762a1bSJed Brown call PetscLogSetThreshold(threshold,oldthreshold,ierr) 77c4762a1bSJed Brown! dummy test of logging a reduction 78c4762a1bSJed Brown ierr = PetscAReduce() 79c4762a1bSJed Brown call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) 80c4762a1bSJed Brown call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) 81c4762a1bSJed Brown if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif 82c4762a1bSJed Brown 83c4762a1bSJed Brown i2 = 2 84c4762a1bSJed Brown i20 = 20 85c4762a1bSJed Brown! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - 86c4762a1bSJed Brown! Create nonlinear solver context 87c4762a1bSJed Brown! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - 88c4762a1bSJed Brown 89c4762a1bSJed Brown call SNESCreate(PETSC_COMM_WORLD,snes,ierr) 90c4762a1bSJed Brown 91c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92c4762a1bSJed Brown! Create matrix and vector data structures; set corresponding routines 93c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 94c4762a1bSJed Brown 95c4762a1bSJed Brown! Create vectors for solution and nonlinear function 96c4762a1bSJed Brown 97c4762a1bSJed Brown call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr) 98c4762a1bSJed Brown call VecDuplicate(x,r,ierr) 99c4762a1bSJed Brown 100c4762a1bSJed Brown! Create Jacobian matrix data structure 101c4762a1bSJed Brown 102c4762a1bSJed Brown call MatCreate(PETSC_COMM_SELF,J,ierr) 103c4762a1bSJed Brown call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr) 104c4762a1bSJed Brown call MatSetFromOptions(J,ierr) 105c4762a1bSJed Brown call MatSetUp(J,ierr) 106c4762a1bSJed Brown 107c4762a1bSJed Brown! Set function evaluation routine and vector 108c4762a1bSJed Brown 109c4762a1bSJed Brown call SNESSetFunction(snes,r,FormFunction,0,ierr) 110c4762a1bSJed Brown 111c4762a1bSJed Brown! Set Jacobian matrix data structure and Jacobian evaluation routine 112c4762a1bSJed Brown 113c4762a1bSJed Brown call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr) 114c4762a1bSJed Brown 115c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 116c4762a1bSJed Brown! Customize nonlinear solver; set runtime options 117c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 118c4762a1bSJed Brown 119c4762a1bSJed Brown! Set linear solver defaults for this problem. By extracting the 120c4762a1bSJed Brown! KSP, KSP, and PC contexts from the SNES context, we can then 121c4762a1bSJed Brown! directly call any KSP, KSP, and PC routines to set various options. 122c4762a1bSJed Brown 123c4762a1bSJed Brown call SNESGetKSP(snes,ksp,ierr) 124c4762a1bSJed Brown call KSPGetPC(ksp,pc,ierr) 125c4762a1bSJed Brown call PCSetType(pc,PCNONE,ierr) 126c4762a1bSJed Brown tol = 1.e-4 127c4762a1bSJed Brown call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, & 128c4762a1bSJed Brown & PETSC_DEFAULT_REAL,i20,ierr) 129c4762a1bSJed Brown 130c4762a1bSJed Brown! Set SNES/KSP/KSP/PC runtime options, e.g., 131c4762a1bSJed Brown! -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 132c4762a1bSJed Brown! These options will override those specified above as long as 133c4762a1bSJed Brown! SNESSetFromOptions() is called _after_ any other customization 134c4762a1bSJed Brown! routines. 135c4762a1bSJed Brown 136c4762a1bSJed Brown 137c4762a1bSJed Brown call SNESSetFromOptions(snes,ierr) 138c4762a1bSJed Brown 139c4762a1bSJed Brown call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 140c4762a1bSJed Brown & '-setls',setls,ierr) 141c4762a1bSJed Brown 142c4762a1bSJed Brown if (setls) then 143c4762a1bSJed Brown call SNESGetLineSearch(snes, linesearch, ierr) 144c4762a1bSJed Brown call SNESLineSearchSetType(linesearch, 'shell', ierr) 145c4762a1bSJed Brown call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch, & 146c4762a1bSJed Brown & 0, ierr) 147c4762a1bSJed Brown endif 148c4762a1bSJed Brown 149c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150c4762a1bSJed Brown! Evaluate initial guess; then solve nonlinear system 151c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152c4762a1bSJed Brown 153c4762a1bSJed Brown! Note: The user should initialize the vector, x, with the initial guess 154c4762a1bSJed Brown! for the nonlinear solver prior to calling SNESSolve(). In particular, 155c4762a1bSJed Brown! to employ an initial guess of zero, the user should explicitly set 156c4762a1bSJed Brown! this vector to zero by calling VecSet(). 157c4762a1bSJed Brown 158c4762a1bSJed Brown pfive = 0.5 159c4762a1bSJed Brown call VecSet(x,pfive,ierr) 160c4762a1bSJed Brown call SNESSolve(snes,PETSC_NULL_VEC,x,ierr) 161c4762a1bSJed Brown call SNESGetIterationNumber(snes,its,ierr); 162c4762a1bSJed Brown if (rank .eq. 0) then 163c4762a1bSJed Brown write(6,100) its 164c4762a1bSJed Brown endif 165c4762a1bSJed Brown 100 format('Number of SNES iterations = ',i5) 166c4762a1bSJed Brown 167c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 168c4762a1bSJed Brown! Free work space. All PETSc objects should be destroyed when they 169c4762a1bSJed Brown! are no longer needed. 170c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 171c4762a1bSJed Brown 172c4762a1bSJed Brown call VecDestroy(x,ierr) 173c4762a1bSJed Brown call VecDestroy(r,ierr) 174c4762a1bSJed Brown call MatDestroy(J,ierr) 175c4762a1bSJed Brown call SNESDestroy(snes,ierr) 176c4762a1bSJed Brown call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr) 177c4762a1bSJed Brown call PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr) 178c4762a1bSJed Brown call PetscLogView(viewer,ierr) 179c4762a1bSJed Brown call PetscViewerDestroy(viewer,ierr) 180c4762a1bSJed Brown call PetscFinalize(ierr) 181c4762a1bSJed Brown end 182c4762a1bSJed Brown! 183c4762a1bSJed Brown! ------------------------------------------------------------------------ 184c4762a1bSJed Brown! 185c4762a1bSJed Brown! FormFunction - Evaluates nonlinear function, F(x). 186c4762a1bSJed Brown! 187c4762a1bSJed Brown! Input Parameters: 188c4762a1bSJed Brown! snes - the SNES context 189c4762a1bSJed Brown! x - input vector 190c4762a1bSJed Brown! dummy - optional user-defined context (not used here) 191c4762a1bSJed Brown! 192c4762a1bSJed Brown! Output Parameter: 193c4762a1bSJed Brown! f - function vector 194c4762a1bSJed Brown! 195c4762a1bSJed Brown subroutine FormFunction(snes,x,f,dummy,ierr) 196c4762a1bSJed Brown use petscsnes 197c4762a1bSJed Brown implicit none 198c4762a1bSJed Brown 199c4762a1bSJed Brown SNES snes 200c4762a1bSJed Brown Vec x,f 201c4762a1bSJed Brown PetscErrorCode ierr 202c4762a1bSJed Brown integer dummy(*) 203c4762a1bSJed Brown 204c4762a1bSJed Brown! Declarations for use with local arrays 205c4762a1bSJed Brown 206c4762a1bSJed Brown PetscScalar lx_v(2),lf_v(2) 207c4762a1bSJed Brown PetscOffset lx_i,lf_i 208c4762a1bSJed Brown 209c4762a1bSJed Brown! Get pointers to vector data. 210c4762a1bSJed Brown! - For default PETSc vectors, VecGetArray() returns a pointer to 211c4762a1bSJed Brown! the data array. Otherwise, the routine is implementation dependent. 212c4762a1bSJed Brown! - You MUST call VecRestoreArray() when you no longer need access to 213c4762a1bSJed Brown! the array. 214c4762a1bSJed Brown! - Note that the Fortran interface to VecGetArray() differs from the 215c4762a1bSJed Brown! C version. See the Fortran chapter of the users manual for details. 216c4762a1bSJed Brown 217c4762a1bSJed Brown call VecGetArrayRead(x,lx_v,lx_i,ierr) 218c4762a1bSJed Brown call VecGetArray(f,lf_v,lf_i,ierr) 219c4762a1bSJed Brown 220c4762a1bSJed Brown! Compute function 221c4762a1bSJed Brown 222c4762a1bSJed Brown lf_a(1) = lx_a(1)*lx_a(1) & 223c4762a1bSJed Brown & + lx_a(1)*lx_a(2) - 3.0 224c4762a1bSJed Brown lf_a(2) = lx_a(1)*lx_a(2) & 225c4762a1bSJed Brown & + lx_a(2)*lx_a(2) - 6.0 226c4762a1bSJed Brown 227c4762a1bSJed Brown! Restore vectors 228c4762a1bSJed Brown 229c4762a1bSJed Brown call VecRestoreArrayRead(x,lx_v,lx_i,ierr) 230c4762a1bSJed Brown call VecRestoreArray(f,lf_v,lf_i,ierr) 231c4762a1bSJed Brown 232c4762a1bSJed Brown return 233c4762a1bSJed Brown end 234c4762a1bSJed Brown 235c4762a1bSJed Brown! --------------------------------------------------------------------- 236c4762a1bSJed Brown! 237c4762a1bSJed Brown! FormJacobian - Evaluates Jacobian matrix. 238c4762a1bSJed Brown! 239c4762a1bSJed Brown! Input Parameters: 240c4762a1bSJed Brown! snes - the SNES context 241c4762a1bSJed Brown! x - input vector 242c4762a1bSJed Brown! dummy - optional user-defined context (not used here) 243c4762a1bSJed Brown! 244c4762a1bSJed Brown! Output Parameters: 245c4762a1bSJed Brown! A - Jacobian matrix 246c4762a1bSJed Brown! B - optionally different preconditioning matrix 247c4762a1bSJed Brown! flag - flag indicating matrix structure 248c4762a1bSJed Brown! 249c4762a1bSJed Brown subroutine FormJacobian(snes,X,jac,B,dummy,ierr) 250c4762a1bSJed Brown use petscsnes 251c4762a1bSJed Brown implicit none 252c4762a1bSJed Brown 253c4762a1bSJed Brown SNES snes 254c4762a1bSJed Brown Vec X 255c4762a1bSJed Brown Mat jac,B 256c4762a1bSJed Brown PetscScalar A(4) 257c4762a1bSJed Brown PetscErrorCode ierr 258c4762a1bSJed Brown PetscInt idx(2),i2 259c4762a1bSJed Brown integer dummy(*) 260c4762a1bSJed Brown 261c4762a1bSJed Brown! Declarations for use with local arrays 262c4762a1bSJed Brown 263c4762a1bSJed Brown PetscScalar lx_v(2) 264c4762a1bSJed Brown PetscOffset lx_i 265c4762a1bSJed Brown 266c4762a1bSJed Brown! Get pointer to vector data 267c4762a1bSJed Brown 268c4762a1bSJed Brown i2 = 2 269c4762a1bSJed Brown call VecGetArrayRead(x,lx_v,lx_i,ierr) 270c4762a1bSJed Brown 271c4762a1bSJed Brown! Compute Jacobian entries and insert into matrix. 272c4762a1bSJed Brown! - Since this is such a small problem, we set all entries for 273c4762a1bSJed Brown! the matrix at once. 274c4762a1bSJed Brown! - Note that MatSetValues() uses 0-based row and column numbers 275c4762a1bSJed Brown! in Fortran as well as in C (as set here in the array idx). 276c4762a1bSJed Brown 277c4762a1bSJed Brown idx(1) = 0 278c4762a1bSJed Brown idx(2) = 1 279c4762a1bSJed Brown A(1) = 2.0*lx_a(1) + lx_a(2) 280c4762a1bSJed Brown A(2) = lx_a(1) 281c4762a1bSJed Brown A(3) = lx_a(2) 282c4762a1bSJed Brown A(4) = lx_a(1) + 2.0*lx_a(2) 283c4762a1bSJed Brown call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr) 284c4762a1bSJed Brown 285c4762a1bSJed Brown! Restore vector 286c4762a1bSJed Brown 287c4762a1bSJed Brown call VecRestoreArrayRead(x,lx_v,lx_i,ierr) 288c4762a1bSJed Brown 289c4762a1bSJed Brown! Assemble matrix 290c4762a1bSJed Brown 291c4762a1bSJed Brown call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr) 292c4762a1bSJed Brown call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr) 293c4762a1bSJed Brown if (B .ne. jac) then 294c4762a1bSJed Brown call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr) 295c4762a1bSJed Brown call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr) 296c4762a1bSJed Brown endif 297c4762a1bSJed Brown 298c4762a1bSJed Brown return 299c4762a1bSJed Brown end 300c4762a1bSJed Brown 301c4762a1bSJed Brown 302c4762a1bSJed Brown subroutine MyLineSearch(linesearch, lctx, ierr) 303c4762a1bSJed Brown use petscsnes 304c4762a1bSJed Brown implicit none 305c4762a1bSJed Brown 306c4762a1bSJed Brown SNESLineSearch linesearch 307c4762a1bSJed Brown SNES snes 308c4762a1bSJed Brown integer lctx 309c4762a1bSJed Brown Vec x, f,g, y, w 310c4762a1bSJed Brown PetscReal ynorm,gnorm,xnorm 311c4762a1bSJed Brown PetscBool flag 312c4762a1bSJed Brown PetscErrorCode ierr 313c4762a1bSJed Brown 314c4762a1bSJed Brown PetscScalar mone 315c4762a1bSJed Brown 316c4762a1bSJed Brown mone = -1.0 317c4762a1bSJed Brown call SNESLineSearchGetSNES(linesearch, snes, ierr) 318c4762a1bSJed Brown call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr) 319c4762a1bSJed Brown call VecNorm(y,NORM_2,ynorm,ierr) 320c4762a1bSJed Brown call VecAXPY(x,mone,y,ierr) 321c4762a1bSJed Brown call SNESComputeFunction(snes,x,f,ierr) 322c4762a1bSJed Brown call VecNorm(f,NORM_2,gnorm,ierr) 323c4762a1bSJed Brown call VecNorm(x,NORM_2,xnorm,ierr) 324c4762a1bSJed Brown call VecNorm(y,NORM_2,ynorm,ierr) 325c4762a1bSJed Brown call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, & 326c4762a1bSJed Brown & ierr) 327c4762a1bSJed Brown flag = PETSC_FALSE 328c4762a1bSJed Brown return 329c4762a1bSJed Brown end 330c4762a1bSJed Brown 331c4762a1bSJed Brown!/*TEST 332c4762a1bSJed Brown! 333c4762a1bSJed Brown! test: 334c4762a1bSJed Brown! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 335c4762a1bSJed Brown! requires: !single 336c4762a1bSJed Brown! 337c4762a1bSJed Brown!TEST*/ 338