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 162! View solver converged reason; we could instead use the option -snes_converged_reason 163 call SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr) 164 165 call SNESGetIterationNumber(snes,its,ierr); 166 if (rank .eq. 0) then 167 write(6,100) its 168 endif 169 100 format('Number of SNES iterations = ',i5) 170 171! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 172! Free work space. All PETSc objects should be destroyed when they 173! are no longer needed. 174! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 175 176 call VecDestroy(x,ierr) 177 call VecDestroy(r,ierr) 178 call MatDestroy(J,ierr) 179 call SNESDestroy(snes,ierr) 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 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 306 subroutine MyLineSearch(linesearch, lctx, ierr) 307 use petscsnes 308 implicit none 309 310 SNESLineSearch linesearch 311 SNES snes 312 integer lctx 313 Vec x, f,g, y, w 314 PetscReal ynorm,gnorm,xnorm 315 PetscBool flag 316 PetscErrorCode ierr 317 318 PetscScalar mone 319 320 mone = -1.0 321 call SNESLineSearchGetSNES(linesearch, snes, ierr) 322 call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr) 323 call VecNorm(y,NORM_2,ynorm,ierr) 324 call VecAXPY(x,mone,y,ierr) 325 call SNESComputeFunction(snes,x,f,ierr) 326 call VecNorm(f,NORM_2,gnorm,ierr) 327 call VecNorm(x,NORM_2,xnorm,ierr) 328 call VecNorm(y,NORM_2,ynorm,ierr) 329 call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, & 330 & ierr) 331 flag = PETSC_FALSE 332 return 333 end 334 335!/*TEST 336! 337! test: 338! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 339! requires: !single 340! 341!TEST*/ 342