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