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