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