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