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! dummy test of logging a reduction 55#if defined(PETSC_USE_LOG) 56 ierr = PetscAReduce() 57#endif 58 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)) 59 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)) 60 PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example') 61 62 i2 = 2 63 i20 = 20 64! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - 65! Create nonlinear solver context 66! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - 67 68 PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,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_DEFAULT_REAL,PETSC_DEFAULT_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(SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,0,ierr)) 122 endif 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! View solver converged reason; we could instead use the option -snes_converged_reason 138 PetscCallA(SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr)) 139 140 PetscCallA(SNESGetIterationNumber(snes,its,ierr)) 141 if (rank .eq. 0) then 142 write(6,100) its 143 endif 144 100 format('Number of SNES iterations = ',i5) 145 146! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147! Free work space. All PETSc objects should be destroyed when they 148! are no longer needed. 149! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 151 PetscCallA(VecDestroy(x,ierr)) 152 PetscCallA(VecDestroy(r,ierr)) 153 PetscCallA(MatDestroy(J,ierr)) 154 PetscCallA(SNESDestroy(snes,ierr)) 155#if defined(PETSC_USE_LOG) 156 PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr)) 157 PetscCallA(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr)) 158 PetscCallA(PetscLogView(viewer,ierr)) 159 PetscCallA(PetscViewerDestroy(viewer,ierr)) 160#endif 161 PetscCallA(PetscFinalize(ierr)) 162 end 163! 164! ------------------------------------------------------------------------ 165! 166! FormFunction - Evaluates nonlinear function, F(x). 167! 168! Input Parameters: 169! snes - the SNES context 170! x - input vector 171! dummy - optional user-defined context (not used here) 172! 173! Output Parameter: 174! f - function vector 175! 176 subroutine FormFunction(snes,x,f,dummy,ierr) 177 use petscsnes 178 implicit none 179 180 SNES snes 181 Vec x,f 182 PetscErrorCode ierr 183 integer dummy(*) 184 185! Declarations for use with local arrays 186 PetscScalar,pointer :: lx_v(:),lf_v(:) 187 188! Get pointers to vector data. 189! - VecGetArrayF90() returns a pointer to the data array. 190! - You MUST call VecRestoreArrayF90() when you no longer need access to 191! the array. 192 193 PetscCall(VecGetArrayReadF90(x,lx_v,ierr)) 194 PetscCall(VecGetArrayF90(f,lf_v,ierr)) 195 196! Compute function 197 198 lf_v(1) = lx_v(1)*lx_v(1) + lx_v(1)*lx_v(2) - 3.0 199 lf_v(2) = lx_v(1)*lx_v(2) + lx_v(2)*lx_v(2) - 6.0 200 201! Restore vectors 202 203 PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr)) 204 PetscCall(VecRestoreArrayF90(f,lf_v,ierr)) 205 206 return 207 end 208 209! --------------------------------------------------------------------- 210! 211! FormJacobian - Evaluates Jacobian matrix. 212! 213! Input Parameters: 214! snes - the SNES context 215! x - input vector 216! dummy - optional user-defined context (not used here) 217! 218! Output Parameters: 219! A - Jacobian matrix 220! B - optionally different preconditioning matrix 221! 222 subroutine FormJacobian(snes,X,jac,B,dummy,ierr) 223 use petscsnes 224 implicit none 225 226 SNES snes 227 Vec X 228 Mat jac,B 229 PetscScalar A(4) 230 PetscErrorCode ierr 231 PetscInt idx(2),i2 232 integer dummy(*) 233 234! Declarations for use with local arrays 235 236 PetscScalar,pointer :: lx_v(:) 237 238! Get pointer to vector data 239 240 i2 = 2 241 PetscCall(VecGetArrayReadF90(x,lx_v,ierr)) 242 243! Compute Jacobian entries and insert into matrix. 244! - Since this is such a small problem, we set all entries for 245! the matrix at once. 246! - Note that MatSetValues() uses 0-based row and column numbers 247! in Fortran as well as in C (as set here in the array idx). 248 249 idx(1) = 0 250 idx(2) = 1 251 A(1) = 2.0*lx_v(1) + lx_v(2) 252 A(2) = lx_v(1) 253 A(3) = lx_v(2) 254 A(4) = lx_v(1) + 2.0*lx_v(2) 255 PetscCall(MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)) 256 257! Restore vector 258 259 PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr)) 260 261! Assemble matrix 262 263 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)) 264 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)) 265 if (B .ne. jac) then 266 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)) 267 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)) 268 endif 269 270 return 271 end 272 273 subroutine MyLineSearch(linesearch, lctx, ierr) 274 use petscsnes 275 implicit none 276 277 SNESLineSearch linesearch 278 SNES snes 279 integer lctx 280 Vec x, f,g, y, w 281 PetscReal ynorm,gnorm,xnorm 282 PetscErrorCode ierr 283 284 PetscScalar mone 285 286 mone = -1.0 287 PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr)) 288 PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)) 289 PetscCall(VecNorm(y,NORM_2,ynorm,ierr)) 290 PetscCall(VecAXPY(x,mone,y,ierr)) 291 PetscCall(SNESComputeFunction(snes,x,f,ierr)) 292 PetscCall(VecNorm(f,NORM_2,gnorm,ierr)) 293 PetscCall(VecNorm(x,NORM_2,xnorm,ierr)) 294 PetscCall(VecNorm(y,NORM_2,ynorm,ierr)) 295 PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,ierr)) 296 return 297 end 298 299!/*TEST 300! 301! test: 302! args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short 303! requires: !single 304! 305!TEST*/ 306