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