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