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 74d21b9a37SPierre 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) 161*91f3e32bSBarry Smith 162*91f3e32bSBarry Smith! View solver converged reason; we could instead use the option -snes_converged_reason 163*91f3e32bSBarry Smith call SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr) 164*91f3e32bSBarry Smith 165c4762a1bSJed Brown call SNESGetIterationNumber(snes,its,ierr); 166c4762a1bSJed Brown if (rank .eq. 0) then 167c4762a1bSJed Brown write(6,100) its 168c4762a1bSJed Brown endif 169c4762a1bSJed Brown 100 format('Number of SNES iterations = ',i5) 170c4762a1bSJed Brown 171c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 172c4762a1bSJed Brown! Free work space. All PETSc objects should be destroyed when they 173c4762a1bSJed Brown! are no longer needed. 174c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 175c4762a1bSJed Brown 176c4762a1bSJed Brown call VecDestroy(x,ierr) 177c4762a1bSJed Brown call VecDestroy(r,ierr) 178c4762a1bSJed Brown call MatDestroy(J,ierr) 179c4762a1bSJed Brown call SNESDestroy(snes,ierr) 180c4762a1bSJed Brown call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr) 181c4762a1bSJed Brown call PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr) 182c4762a1bSJed Brown call PetscLogView(viewer,ierr) 183c4762a1bSJed Brown call PetscViewerDestroy(viewer,ierr) 184c4762a1bSJed Brown call PetscFinalize(ierr) 185c4762a1bSJed Brown end 186c4762a1bSJed Brown! 187c4762a1bSJed Brown! ------------------------------------------------------------------------ 188c4762a1bSJed Brown! 189c4762a1bSJed Brown! FormFunction - Evaluates nonlinear function, F(x). 190c4762a1bSJed Brown! 191c4762a1bSJed Brown! Input Parameters: 192c4762a1bSJed Brown! snes - the SNES context 193c4762a1bSJed Brown! x - input vector 194c4762a1bSJed Brown! dummy - optional user-defined context (not used here) 195c4762a1bSJed Brown! 196c4762a1bSJed Brown! Output Parameter: 197c4762a1bSJed Brown! f - function vector 198c4762a1bSJed Brown! 199c4762a1bSJed Brown subroutine FormFunction(snes,x,f,dummy,ierr) 200c4762a1bSJed Brown use petscsnes 201c4762a1bSJed Brown implicit none 202c4762a1bSJed Brown 203c4762a1bSJed Brown SNES snes 204c4762a1bSJed Brown Vec x,f 205c4762a1bSJed Brown PetscErrorCode ierr 206c4762a1bSJed Brown integer dummy(*) 207c4762a1bSJed Brown 208c4762a1bSJed Brown! Declarations for use with local arrays 209c4762a1bSJed Brown 210c4762a1bSJed Brown PetscScalar lx_v(2),lf_v(2) 211c4762a1bSJed Brown PetscOffset lx_i,lf_i 212c4762a1bSJed Brown 213c4762a1bSJed Brown! Get pointers to vector data. 214c4762a1bSJed Brown! - For default PETSc vectors, VecGetArray() returns a pointer to 215c4762a1bSJed Brown! the data array. Otherwise, the routine is implementation dependent. 216c4762a1bSJed Brown! - You MUST call VecRestoreArray() when you no longer need access to 217c4762a1bSJed Brown! the array. 218c4762a1bSJed Brown! - Note that the Fortran interface to VecGetArray() differs from the 219c4762a1bSJed Brown! C version. See the Fortran chapter of the users manual for details. 220c4762a1bSJed Brown 221c4762a1bSJed Brown call VecGetArrayRead(x,lx_v,lx_i,ierr) 222c4762a1bSJed Brown call VecGetArray(f,lf_v,lf_i,ierr) 223c4762a1bSJed Brown 224c4762a1bSJed Brown! Compute function 225c4762a1bSJed Brown 226c4762a1bSJed Brown lf_a(1) = lx_a(1)*lx_a(1) & 227c4762a1bSJed Brown & + lx_a(1)*lx_a(2) - 3.0 228c4762a1bSJed Brown lf_a(2) = lx_a(1)*lx_a(2) & 229c4762a1bSJed Brown & + lx_a(2)*lx_a(2) - 6.0 230c4762a1bSJed Brown 231c4762a1bSJed Brown! Restore vectors 232c4762a1bSJed Brown 233c4762a1bSJed Brown call VecRestoreArrayRead(x,lx_v,lx_i,ierr) 234c4762a1bSJed Brown call VecRestoreArray(f,lf_v,lf_i,ierr) 235c4762a1bSJed Brown 236c4762a1bSJed Brown return 237c4762a1bSJed Brown end 238c4762a1bSJed Brown 239c4762a1bSJed Brown! --------------------------------------------------------------------- 240c4762a1bSJed Brown! 241c4762a1bSJed Brown! FormJacobian - Evaluates Jacobian matrix. 242c4762a1bSJed Brown! 243c4762a1bSJed Brown! Input Parameters: 244c4762a1bSJed Brown! snes - the SNES context 245c4762a1bSJed Brown! x - input vector 246c4762a1bSJed Brown! dummy - optional user-defined context (not used here) 247c4762a1bSJed Brown! 248c4762a1bSJed Brown! Output Parameters: 249c4762a1bSJed Brown! A - Jacobian matrix 250c4762a1bSJed Brown! B - optionally different preconditioning matrix 251c4762a1bSJed Brown! flag - flag indicating matrix structure 252c4762a1bSJed Brown! 253c4762a1bSJed Brown subroutine FormJacobian(snes,X,jac,B,dummy,ierr) 254c4762a1bSJed Brown use petscsnes 255c4762a1bSJed Brown implicit none 256c4762a1bSJed Brown 257c4762a1bSJed Brown SNES snes 258c4762a1bSJed Brown Vec X 259c4762a1bSJed Brown Mat jac,B 260c4762a1bSJed Brown PetscScalar A(4) 261c4762a1bSJed Brown PetscErrorCode ierr 262c4762a1bSJed Brown PetscInt idx(2),i2 263c4762a1bSJed Brown integer dummy(*) 264c4762a1bSJed Brown 265c4762a1bSJed Brown! Declarations for use with local arrays 266c4762a1bSJed Brown 267c4762a1bSJed Brown PetscScalar lx_v(2) 268c4762a1bSJed Brown PetscOffset lx_i 269c4762a1bSJed Brown 270c4762a1bSJed Brown! Get pointer to vector data 271c4762a1bSJed Brown 272c4762a1bSJed Brown i2 = 2 273c4762a1bSJed Brown call VecGetArrayRead(x,lx_v,lx_i,ierr) 274c4762a1bSJed Brown 275c4762a1bSJed Brown! Compute Jacobian entries and insert into matrix. 276c4762a1bSJed Brown! - Since this is such a small problem, we set all entries for 277c4762a1bSJed Brown! the matrix at once. 278c4762a1bSJed Brown! - Note that MatSetValues() uses 0-based row and column numbers 279c4762a1bSJed Brown! in Fortran as well as in C (as set here in the array idx). 280c4762a1bSJed Brown 281c4762a1bSJed Brown idx(1) = 0 282c4762a1bSJed Brown idx(2) = 1 283c4762a1bSJed Brown A(1) = 2.0*lx_a(1) + lx_a(2) 284c4762a1bSJed Brown A(2) = lx_a(1) 285c4762a1bSJed Brown A(3) = lx_a(2) 286c4762a1bSJed Brown A(4) = lx_a(1) + 2.0*lx_a(2) 287c4762a1bSJed Brown call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr) 288c4762a1bSJed Brown 289c4762a1bSJed Brown! Restore vector 290c4762a1bSJed Brown 291c4762a1bSJed Brown call VecRestoreArrayRead(x,lx_v,lx_i,ierr) 292c4762a1bSJed Brown 293c4762a1bSJed Brown! Assemble matrix 294c4762a1bSJed Brown 295c4762a1bSJed Brown call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr) 296c4762a1bSJed Brown call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr) 297c4762a1bSJed Brown if (B .ne. jac) then 298c4762a1bSJed Brown call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr) 299c4762a1bSJed Brown call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr) 300c4762a1bSJed Brown endif 301c4762a1bSJed Brown 302c4762a1bSJed Brown return 303c4762a1bSJed Brown end 304c4762a1bSJed Brown 305c4762a1bSJed Brown 306c4762a1bSJed Brown subroutine MyLineSearch(linesearch, lctx, ierr) 307c4762a1bSJed Brown use petscsnes 308c4762a1bSJed Brown implicit none 309c4762a1bSJed Brown 310c4762a1bSJed Brown SNESLineSearch linesearch 311c4762a1bSJed Brown SNES snes 312c4762a1bSJed Brown integer lctx 313c4762a1bSJed Brown Vec x, f,g, y, w 314c4762a1bSJed Brown PetscReal ynorm,gnorm,xnorm 315c4762a1bSJed Brown PetscBool flag 316c4762a1bSJed Brown PetscErrorCode ierr 317c4762a1bSJed Brown 318c4762a1bSJed Brown PetscScalar mone 319c4762a1bSJed Brown 320c4762a1bSJed Brown mone = -1.0 321c4762a1bSJed Brown call SNESLineSearchGetSNES(linesearch, snes, ierr) 322c4762a1bSJed Brown call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr) 323c4762a1bSJed Brown call VecNorm(y,NORM_2,ynorm,ierr) 324c4762a1bSJed Brown call VecAXPY(x,mone,y,ierr) 325c4762a1bSJed Brown call SNESComputeFunction(snes,x,f,ierr) 326c4762a1bSJed Brown call VecNorm(f,NORM_2,gnorm,ierr) 327c4762a1bSJed Brown call VecNorm(x,NORM_2,xnorm,ierr) 328c4762a1bSJed Brown call VecNorm(y,NORM_2,ynorm,ierr) 329c4762a1bSJed Brown call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, & 330c4762a1bSJed Brown & ierr) 331c4762a1bSJed Brown flag = PETSC_FALSE 332c4762a1bSJed Brown return 333c4762a1bSJed Brown end 334c4762a1bSJed Brown 335c4762a1bSJed Brown!/*TEST 336c4762a1bSJed Brown! 337c4762a1bSJed Brown! test: 338c4762a1bSJed Brown! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 339c4762a1bSJed Brown! requires: !single 340c4762a1bSJed Brown! 341c4762a1bSJed Brown!TEST*/ 342