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(SNESLineSearchShellSetApply(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 end 203 204! --------------------------------------------------------------------- 205! 206! FormJacobian - Evaluates Jacobian matrix. 207! 208! Input Parameters: 209! snes - the SNES context 210! x - input vector 211! dummy - optional user-defined context (not used here) 212! 213! Output Parameters: 214! A - Jacobian matrix 215! B - optionally different preconditioning matrix 216! 217 subroutine FormJacobian(snes,X,jac,B,dummy,ierr) 218 use petscsnes 219 implicit none 220 221 SNES snes 222 Vec X 223 Mat jac,B 224 PetscScalar A(4) 225 PetscErrorCode ierr 226 PetscInt idx(2),i2 227 integer dummy(*) 228 229! Declarations for use with local arrays 230 231 PetscScalar,pointer :: lx_v(:) 232 233! Get pointer to vector data 234 235 i2 = 2 236 PetscCall(VecGetArrayReadF90(x,lx_v,ierr)) 237 238! Compute Jacobian entries and insert into matrix. 239! - Since this is such a small problem, we set all entries for 240! the matrix at once. 241! - Note that MatSetValues() uses 0-based row and column numbers 242! in Fortran as well as in C (as set here in the array idx). 243 244 idx(1) = 0 245 idx(2) = 1 246 A(1) = 2.0*lx_v(1) + lx_v(2) 247 A(2) = lx_v(1) 248 A(3) = lx_v(2) 249 A(4) = lx_v(1) + 2.0*lx_v(2) 250 PetscCall(MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)) 251 252! Restore vector 253 254 PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr)) 255 256! Assemble matrix 257 258 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)) 259 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)) 260 if (B .ne. jac) then 261 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 262 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 263 endif 264 265 end 266 267 subroutine MyLineSearch(linesearch, lctx, ierr) 268 use petscsnes 269 implicit none 270 271 SNESLineSearch linesearch 272 SNES snes 273 integer lctx 274 Vec x, f,g, y, w 275 PetscReal ynorm,gnorm,xnorm 276 PetscErrorCode ierr 277 278 PetscScalar mone 279 280 mone = -1.0 281 PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr)) 282 PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)) 283 PetscCall(VecNorm(y,NORM_2,ynorm,ierr)) 284 PetscCall(VecAXPY(x,mone,y,ierr)) 285 PetscCall(SNESComputeFunction(snes,x,f,ierr)) 286 PetscCall(VecNorm(f,NORM_2,gnorm,ierr)) 287 PetscCall(VecNorm(x,NORM_2,xnorm,ierr)) 288 PetscCall(VecNorm(y,NORM_2,ynorm,ierr)) 289 PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,ierr)) 290 end 291 292!/*TEST 293! 294! test: 295! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 296! requires: !single 297! 298!TEST*/ 299