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