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