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