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