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! flag - flag indicating matrix structure 252! 253 subroutine FormJacobian(snes,X,jac,B,dummy,ierr) 254 use petscsnes 255 implicit none 256 257 SNES snes 258 Vec X 259 Mat jac,B 260 PetscScalar A(4) 261 PetscErrorCode ierr 262 PetscInt idx(2),i2 263 integer dummy(*) 264 265! Declarations for use with local arrays 266 267 PetscScalar lx_v(2) 268 PetscOffset lx_i 269 270! Get pointer to vector data 271 272 i2 = 2 273 call VecGetArrayRead(x,lx_v,lx_i,ierr) 274 275! Compute Jacobian entries and insert into matrix. 276! - Since this is such a small problem, we set all entries for 277! the matrix at once. 278! - Note that MatSetValues() uses 0-based row and column numbers 279! in Fortran as well as in C (as set here in the array idx). 280 281 idx(1) = 0 282 idx(2) = 1 283 A(1) = 2.0*lx_a(1) + lx_a(2) 284 A(2) = lx_a(1) 285 A(3) = lx_a(2) 286 A(4) = lx_a(1) + 2.0*lx_a(2) 287 call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr) 288 289! Restore vector 290 291 call VecRestoreArrayRead(x,lx_v,lx_i,ierr) 292 293! Assemble matrix 294 295 call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr) 296 call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr) 297 if (B .ne. jac) then 298 call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr) 299 call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr) 300 endif 301 302 return 303 end 304 305 subroutine MyLineSearch(linesearch, lctx, ierr) 306 use petscsnes 307 implicit none 308 309 SNESLineSearch linesearch 310 SNES snes 311 integer lctx 312 Vec x, f,g, y, w 313 PetscReal ynorm,gnorm,xnorm 314 PetscBool flag 315 PetscErrorCode ierr 316 317 PetscScalar mone 318 319 mone = -1.0 320 call SNESLineSearchGetSNES(linesearch, snes, ierr) 321 call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr) 322 call VecNorm(y,NORM_2,ynorm,ierr) 323 call VecAXPY(x,mone,y,ierr) 324 call SNESComputeFunction(snes,x,f,ierr) 325 call VecNorm(f,NORM_2,gnorm,ierr) 326 call VecNorm(x,NORM_2,xnorm,ierr) 327 call VecNorm(y,NORM_2,ynorm,ierr) 328 call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, & 329 & ierr) 330 flag = PETSC_FALSE 331 return 332 end 333 334!/*TEST 335! 336! test: 337! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 338! requires: !single 339! 340!TEST*/ 341