1c4762a1bSJed Brown! 2c4762a1bSJed Brown! Description: This example solves a nonlinear system on 1 processor with SNES. 3c4762a1bSJed Brown! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular 4c4762a1bSJed Brown! domain. The command line options include: 5c4762a1bSJed Brown! -par <parameter>, where <parameter> indicates the nonlinearity of the problem 6c4762a1bSJed Brown! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 7c4762a1bSJed Brown! -mx <xg>, where <xg> = number of grid points in the x-direction 8c4762a1bSJed Brown! -my <yg>, where <yg> = number of grid points in the y-direction 9c4762a1bSJed Brown! 10c4762a1bSJed Brown!!/*T 11c4762a1bSJed Brown! Concepts: SNES^sequential Bratu example 12c4762a1bSJed Brown! Processors: 1 13c4762a1bSJed Brown!T*/ 14c4762a1bSJed Brown 15c4762a1bSJed Brown 16c4762a1bSJed Brown! 17c4762a1bSJed Brown! -------------------------------------------------------------------------- 18c4762a1bSJed Brown! 19c4762a1bSJed Brown! Solid Fuel Ignition (SFI) problem. This problem is modeled by 20c4762a1bSJed Brown! the partial differential equation 21c4762a1bSJed Brown! 22c4762a1bSJed Brown! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 23c4762a1bSJed Brown! 24c4762a1bSJed Brown! with boundary conditions 25c4762a1bSJed Brown! 26c4762a1bSJed Brown! u = 0 for x = 0, x = 1, y = 0, y = 1. 27c4762a1bSJed Brown! 28c4762a1bSJed Brown! A finite difference approximation with the usual 5-point stencil 29c4762a1bSJed Brown! is used to discretize the boundary value problem to obtain a nonlinear 30c4762a1bSJed Brown! system of equations. 31c4762a1bSJed Brown! 32c4762a1bSJed Brown! The parallel version of this code is snes/tutorials/ex5f.F 33c4762a1bSJed Brown! 34c4762a1bSJed Brown! -------------------------------------------------------------------------- 35c4762a1bSJed Brown subroutine postcheck(snes,x,y,w,changed_y,changed_w,ctx,ierr) 36c4762a1bSJed Brown#include <petsc/finclude/petscsnes.h> 37c4762a1bSJed Brown use petscsnes 38c4762a1bSJed Brown implicit none 39c4762a1bSJed Brown SNES snes 40c4762a1bSJed Brown PetscReal norm 41c4762a1bSJed Brown Vec tmp,x,y,w 42c4762a1bSJed Brown PetscBool changed_w,changed_y 43c4762a1bSJed Brown PetscErrorCode ierr 44c4762a1bSJed Brown PetscInt ctx 45c4762a1bSJed Brown PetscScalar mone 46c4762a1bSJed Brown 47c4762a1bSJed Brown call VecDuplicate(x,tmp,ierr) 48c4762a1bSJed Brown mone = -1.0 49c4762a1bSJed Brown call VecWAXPY(tmp,mone,x,w,ierr) 50c4762a1bSJed Brown call VecNorm(tmp,NORM_2,norm,ierr) 51c4762a1bSJed Brown call VecDestroy(tmp,ierr) 52c4762a1bSJed Brown print*, 'Norm of search step ',norm 53c4762a1bSJed Brown changed_y = PETSC_FALSE 54c4762a1bSJed Brown changed_w = PETSC_FALSE 55c4762a1bSJed Brown return 56c4762a1bSJed Brown end 57c4762a1bSJed Brown 58c4762a1bSJed Brown program main 59c4762a1bSJed Brown#include <petsc/finclude/petscdraw.h> 60c4762a1bSJed Brown use petscsnes 61c4762a1bSJed Brown implicit none 62*17a42bb7SSatish Balay interface SNESSetJacobian 63*17a42bb7SSatish Balay subroutine SNESSetJacobian1(a,b,c,d,e,z) 64*17a42bb7SSatish Balay use petscsnes 65*17a42bb7SSatish Balay SNES a 66*17a42bb7SSatish Balay Mat b 67*17a42bb7SSatish Balay Mat c 68*17a42bb7SSatish Balay external d 69*17a42bb7SSatish Balay MatFDColoring e 70*17a42bb7SSatish Balay PetscErrorCode z 71*17a42bb7SSatish Balay end subroutine 72*17a42bb7SSatish Balay subroutine SNESSetJacobian2(a,b,c,d,e,z) 73*17a42bb7SSatish Balay use petscsnes 74*17a42bb7SSatish Balay SNES a 75*17a42bb7SSatish Balay Mat b 76*17a42bb7SSatish Balay Mat c 77*17a42bb7SSatish Balay external d 78*17a42bb7SSatish Balay integer e 79*17a42bb7SSatish Balay PetscErrorCode z 80*17a42bb7SSatish Balay end subroutine 81*17a42bb7SSatish Balay end interface 82c4762a1bSJed Brown! 83c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 84c4762a1bSJed Brown! Variable declarations 85c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86c4762a1bSJed Brown! 87c4762a1bSJed Brown! Variables: 88c4762a1bSJed Brown! snes - nonlinear solver 89c4762a1bSJed Brown! x, r - solution, residual vectors 90c4762a1bSJed Brown! J - Jacobian matrix 91c4762a1bSJed Brown! its - iterations for convergence 92c4762a1bSJed Brown! matrix_free - flag - 1 indicates matrix-free version 93c4762a1bSJed Brown! lambda - nonlinearity parameter 94c4762a1bSJed Brown! draw - drawing context 95c4762a1bSJed Brown! 96c4762a1bSJed Brown SNES snes 97c4762a1bSJed Brown MatColoring mc 98c4762a1bSJed Brown Vec x,r 99c4762a1bSJed Brown PetscDraw draw 100c4762a1bSJed Brown Mat J 101c4762a1bSJed Brown PetscBool matrix_free,flg,fd_coloring 102c4762a1bSJed Brown PetscErrorCode ierr 103c4762a1bSJed Brown PetscInt its,N, mx,my,i5 104c4762a1bSJed Brown PetscMPIInt size,rank 105c4762a1bSJed Brown PetscReal lambda_max,lambda_min,lambda 106c4762a1bSJed Brown MatFDColoring fdcoloring 107c4762a1bSJed Brown ISColoring iscoloring 108c4762a1bSJed Brown PetscBool pc 109c4762a1bSJed Brown external postcheck 110c4762a1bSJed Brown 111c4762a1bSJed Brown PetscScalar lx_v(0:1) 112c4762a1bSJed Brown PetscOffset lx_i 113c4762a1bSJed Brown 114c4762a1bSJed Brown! Store parameters in common block 115c4762a1bSJed Brown 116c4762a1bSJed Brown common /params/ lambda,mx,my,fd_coloring 117c4762a1bSJed Brown 118c4762a1bSJed Brown! Note: Any user-defined Fortran routines (such as FormJacobian) 119c4762a1bSJed Brown! MUST be declared as external. 120c4762a1bSJed Brown 121c4762a1bSJed Brown external FormFunction,FormInitialGuess,FormJacobian 122c4762a1bSJed Brown 123c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 124c4762a1bSJed Brown! Initialize program 125c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 126c4762a1bSJed Brown 127c4762a1bSJed Brown call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 128c4762a1bSJed Brown if (ierr .ne. 0) then 129c4762a1bSJed Brown print*,'Unable to initialize PETSc' 130c4762a1bSJed Brown stop 131c4762a1bSJed Brown endif 132c4762a1bSJed Brown call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) 133c4762a1bSJed Brown call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) 134c4762a1bSJed Brown 135c4762a1bSJed Brown if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif 136c4762a1bSJed Brown 137c4762a1bSJed Brown! Initialize problem parameters 138c4762a1bSJed Brown i5 = 5 139c4762a1bSJed Brown lambda_max = 6.81 140c4762a1bSJed Brown lambda_min = 0.0 141c4762a1bSJed Brown lambda = 6.0 142c4762a1bSJed Brown mx = 4 143c4762a1bSJed Brown my = 4 144c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr) 145c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr) 146c4762a1bSJed Brown call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr) 147c4762a1bSJed Brown if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda out of range '); endif 148c4762a1bSJed Brown N = mx*my 149c4762a1bSJed Brown pc = PETSC_FALSE 150c4762a1bSJed Brown call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-pc',pc,PETSC_NULL_BOOL,ierr); 151c4762a1bSJed Brown 152c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 153c4762a1bSJed Brown! Create nonlinear solver context 154c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155c4762a1bSJed Brown 156c4762a1bSJed Brown call SNESCreate(PETSC_COMM_WORLD,snes,ierr) 157c4762a1bSJed Brown 158c4762a1bSJed Brown if (pc .eqv. PETSC_TRUE) then 159c4762a1bSJed Brown call SNESSetType(snes,SNESNEWTONTR,ierr) 160c4762a1bSJed Brown call SNESNewtonTRSetPostCheck(snes, postcheck,snes,ierr) 161c4762a1bSJed Brown endif 162c4762a1bSJed Brown 163c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 164c4762a1bSJed Brown! Create vector data structures; set function evaluation routine 165c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 166c4762a1bSJed Brown 167c4762a1bSJed Brown call VecCreate(PETSC_COMM_WORLD,x,ierr) 168c4762a1bSJed Brown call VecSetSizes(x,PETSC_DECIDE,N,ierr) 169c4762a1bSJed Brown call VecSetFromOptions(x,ierr) 170c4762a1bSJed Brown call VecDuplicate(x,r,ierr) 171c4762a1bSJed Brown 172c4762a1bSJed Brown! Set function evaluation routine and vector. Whenever the nonlinear 173c4762a1bSJed Brown! solver needs to evaluate the nonlinear function, it will call this 174c4762a1bSJed Brown! routine. 175c4762a1bSJed Brown! - Note that the final routine argument is the user-defined 176c4762a1bSJed Brown! context that provides application-specific data for the 177c4762a1bSJed Brown! function evaluation routine. 178c4762a1bSJed Brown 179c4762a1bSJed Brown call SNESSetFunction(snes,r,FormFunction,fdcoloring,ierr) 180c4762a1bSJed Brown 181c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 182c4762a1bSJed Brown! Create matrix data structure; set Jacobian evaluation routine 183c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184c4762a1bSJed Brown 185c4762a1bSJed Brown! Create matrix. Here we only approximately preallocate storage space 186c4762a1bSJed Brown! for the Jacobian. See the users manual for a discussion of better 187c4762a1bSJed Brown! techniques for preallocating matrix memory. 188c4762a1bSJed Brown 189c4762a1bSJed Brown call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_mf',matrix_free,ierr) 190c4762a1bSJed Brown if (.not. matrix_free) then 191c4762a1bSJed Brown call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,i5,PETSC_NULL_INTEGER,J,ierr) 192c4762a1bSJed Brown endif 193c4762a1bSJed Brown 194c4762a1bSJed Brown! 195c4762a1bSJed Brown! This option will cause the Jacobian to be computed via finite differences 196c4762a1bSJed Brown! efficiently using a coloring of the columns of the matrix. 197c4762a1bSJed Brown! 198c4762a1bSJed Brown fd_coloring = .false. 199c4762a1bSJed Brown call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-snes_fd_coloring',fd_coloring,ierr) 200c4762a1bSJed Brown if (fd_coloring) then 201c4762a1bSJed Brown 202c4762a1bSJed Brown! 203c4762a1bSJed Brown! This initializes the nonzero structure of the Jacobian. This is artificial 204c4762a1bSJed Brown! because clearly if we had a routine to compute the Jacobian we won't need 205c4762a1bSJed Brown! to use finite differences. 206c4762a1bSJed Brown! 207c4762a1bSJed Brown call FormJacobian(snes,x,J,J,0,ierr) 208c4762a1bSJed Brown! 209c4762a1bSJed Brown! Color the matrix, i.e. determine groups of columns that share no common 210c4762a1bSJed Brown! rows. These columns in the Jacobian can all be computed simulataneously. 211c4762a1bSJed Brown! 212c4762a1bSJed Brown call MatColoringCreate(J,mc,ierr) 213c4762a1bSJed Brown call MatColoringSetType(mc,MATCOLORINGNATURAL,ierr) 214c4762a1bSJed Brown call MatColoringSetFromOptions(mc,ierr) 215c4762a1bSJed Brown call MatColoringApply(mc,iscoloring,ierr) 216c4762a1bSJed Brown call MatColoringDestroy(mc,ierr) 217c4762a1bSJed Brown! 218c4762a1bSJed Brown! Create the data structure that SNESComputeJacobianDefaultColor() uses 219c4762a1bSJed Brown! to compute the actual Jacobians via finite differences. 220c4762a1bSJed Brown! 221c4762a1bSJed Brown call MatFDColoringCreate(J,iscoloring,fdcoloring,ierr) 222c4762a1bSJed Brown call MatFDColoringSetFunction(fdcoloring,FormFunction,fdcoloring,ierr) 223c4762a1bSJed Brown call MatFDColoringSetFromOptions(fdcoloring,ierr) 224c4762a1bSJed Brown call MatFDColoringSetUp(J,iscoloring,fdcoloring,ierr) 225c4762a1bSJed Brown! 226c4762a1bSJed Brown! Tell SNES to use the routine SNESComputeJacobianDefaultColor() 227c4762a1bSJed Brown! to compute Jacobians. 228c4762a1bSJed Brown! 229c4762a1bSJed Brown call SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring,ierr) 230c4762a1bSJed Brown call ISColoringDestroy(iscoloring,ierr) 231c4762a1bSJed Brown 232c4762a1bSJed Brown else if (.not. matrix_free) then 233c4762a1bSJed Brown 234c4762a1bSJed Brown! Set Jacobian matrix data structure and default Jacobian evaluation 235c4762a1bSJed Brown! routine. Whenever the nonlinear solver needs to compute the 236c4762a1bSJed Brown! Jacobian matrix, it will call this routine. 237c4762a1bSJed Brown! - Note that the final routine argument is the user-defined 238c4762a1bSJed Brown! context that provides application-specific data for the 239c4762a1bSJed Brown! Jacobian evaluation routine. 240c4762a1bSJed Brown! - The user can override with: 241c4762a1bSJed Brown! -snes_fd : default finite differencing approximation of Jacobian 242c4762a1bSJed Brown! -snes_mf : matrix-free Newton-Krylov method with no preconditioning 243c4762a1bSJed Brown! (unless user explicitly sets preconditioner) 244c4762a1bSJed Brown! -snes_mf_operator : form preconditioning matrix as set by the user, 245c4762a1bSJed Brown! but use matrix-free approx for Jacobian-vector 246c4762a1bSJed Brown! products within Newton-Krylov method 247c4762a1bSJed Brown! 248c4762a1bSJed Brown call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr) 249c4762a1bSJed Brown endif 250c4762a1bSJed Brown 251c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 252c4762a1bSJed Brown! Customize nonlinear solver; set runtime options 253c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 254c4762a1bSJed Brown 255c4762a1bSJed Brown! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 256c4762a1bSJed Brown 257c4762a1bSJed Brown call SNESSetFromOptions(snes,ierr) 258c4762a1bSJed Brown 259c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 260c4762a1bSJed Brown! Evaluate initial guess; then solve nonlinear system. 261c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 262c4762a1bSJed Brown 263c4762a1bSJed Brown! Note: The user should initialize the vector, x, with the initial guess 264c4762a1bSJed Brown! for the nonlinear solver prior to calling SNESSolve(). In particular, 265c4762a1bSJed Brown! to employ an initial guess of zero, the user should explicitly set 266c4762a1bSJed Brown! this vector to zero by calling VecSet(). 267c4762a1bSJed Brown 268c4762a1bSJed Brown call FormInitialGuess(x,ierr) 269c4762a1bSJed Brown call SNESSolve(snes,PETSC_NULL_VEC,x,ierr) 270c4762a1bSJed Brown call SNESGetIterationNumber(snes,its,ierr); 271c4762a1bSJed Brown if (rank .eq. 0) then 272c4762a1bSJed Brown write(6,100) its 273c4762a1bSJed Brown endif 274c4762a1bSJed Brown 100 format('Number of SNES iterations = ',i1) 275c4762a1bSJed Brown 276c4762a1bSJed Brown! PetscDraw contour plot of solution 277c4762a1bSJed Brown 278c4762a1bSJed Brown call PetscDrawCreate(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'Solution',300,0,300,300,draw,ierr) 279c4762a1bSJed Brown call PetscDrawSetFromOptions(draw,ierr) 280c4762a1bSJed Brown 281c4762a1bSJed Brown call VecGetArrayRead(x,lx_v,lx_i,ierr) 282c4762a1bSJed Brown call PetscDrawTensorContour(draw,mx,my,PETSC_NULL_REAL,PETSC_NULL_REAL,lx_v(lx_i+1),ierr) 283c4762a1bSJed Brown call VecRestoreArrayRead(x,lx_v,lx_i,ierr) 284c4762a1bSJed Brown 285c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 286c4762a1bSJed Brown! Free work space. All PETSc objects should be destroyed when they 287c4762a1bSJed Brown! are no longer needed. 288c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 289c4762a1bSJed Brown 290c4762a1bSJed Brown if (.not. matrix_free) call MatDestroy(J,ierr) 291c4762a1bSJed Brown if (fd_coloring) call MatFDColoringDestroy(fdcoloring,ierr) 292c4762a1bSJed Brown 293c4762a1bSJed Brown call VecDestroy(x,ierr) 294c4762a1bSJed Brown call VecDestroy(r,ierr) 295c4762a1bSJed Brown call SNESDestroy(snes,ierr) 296c4762a1bSJed Brown call PetscDrawDestroy(draw,ierr) 297c4762a1bSJed Brown call PetscFinalize(ierr) 298c4762a1bSJed Brown end 299c4762a1bSJed Brown 300c4762a1bSJed Brown! --------------------------------------------------------------------- 301c4762a1bSJed Brown! 302c4762a1bSJed Brown! FormInitialGuess - Forms initial approximation. 303c4762a1bSJed Brown! 304c4762a1bSJed Brown! Input Parameter: 305c4762a1bSJed Brown! X - vector 306c4762a1bSJed Brown! 307c4762a1bSJed Brown! Output Parameters: 308c4762a1bSJed Brown! X - vector 309c4762a1bSJed Brown! ierr - error code 310c4762a1bSJed Brown! 311c4762a1bSJed Brown! Notes: 312c4762a1bSJed Brown! This routine serves as a wrapper for the lower-level routine 313c4762a1bSJed Brown! "ApplicationInitialGuess", where the actual computations are 314c4762a1bSJed Brown! done using the standard Fortran style of treating the local 315c4762a1bSJed Brown! vector data as a multidimensional array over the local mesh. 316c4762a1bSJed Brown! This routine merely accesses the local vector data via 317c4762a1bSJed Brown! VecGetArray() and VecRestoreArray(). 318c4762a1bSJed Brown! 319c4762a1bSJed Brown subroutine FormInitialGuess(X,ierr) 320c4762a1bSJed Brown use petscsnes 321c4762a1bSJed Brown implicit none 322c4762a1bSJed Brown 323c4762a1bSJed Brown! Input/output variables: 324c4762a1bSJed Brown Vec X 325c4762a1bSJed Brown PetscErrorCode ierr 326c4762a1bSJed Brown 327c4762a1bSJed Brown! Declarations for use with local arrays: 328c4762a1bSJed Brown PetscScalar lx_v(0:1) 329c4762a1bSJed Brown PetscOffset lx_i 330c4762a1bSJed Brown 331c4762a1bSJed Brown ierr = 0 332c4762a1bSJed Brown 333c4762a1bSJed Brown! Get a pointer to vector data. 334c4762a1bSJed Brown! - For default PETSc vectors, VecGetArray() returns a pointer to 335c4762a1bSJed Brown! the data array. Otherwise, the routine is implementation dependent. 336c4762a1bSJed Brown! - You MUST call VecRestoreArray() when you no longer need access to 337c4762a1bSJed Brown! the array. 338c4762a1bSJed Brown! - Note that the Fortran interface to VecGetArray() differs from the 339c4762a1bSJed Brown! C version. See the users manual for details. 340c4762a1bSJed Brown 341c4762a1bSJed Brown call VecGetArray(X,lx_v,lx_i,ierr) 342c4762a1bSJed Brown 343c4762a1bSJed Brown! Compute initial guess 344c4762a1bSJed Brown 345c4762a1bSJed Brown call ApplicationInitialGuess(lx_v(lx_i),ierr) 346c4762a1bSJed Brown 347c4762a1bSJed Brown! Restore vector 348c4762a1bSJed Brown 349c4762a1bSJed Brown call VecRestoreArray(X,lx_v,lx_i,ierr) 350c4762a1bSJed Brown 351c4762a1bSJed Brown return 352c4762a1bSJed Brown end 353c4762a1bSJed Brown 354c4762a1bSJed Brown! --------------------------------------------------------------------- 355c4762a1bSJed Brown! 356c4762a1bSJed Brown! ApplicationInitialGuess - Computes initial approximation, called by 357c4762a1bSJed Brown! the higher level routine FormInitialGuess(). 358c4762a1bSJed Brown! 359c4762a1bSJed Brown! Input Parameter: 360c4762a1bSJed Brown! x - local vector data 361c4762a1bSJed Brown! 362c4762a1bSJed Brown! Output Parameters: 363c4762a1bSJed Brown! f - local vector data, f(x) 364c4762a1bSJed Brown! ierr - error code 365c4762a1bSJed Brown! 366c4762a1bSJed Brown! Notes: 367c4762a1bSJed Brown! This routine uses standard Fortran-style computations over a 2-dim array. 368c4762a1bSJed Brown! 369c4762a1bSJed Brown subroutine ApplicationInitialGuess(x,ierr) 370c4762a1bSJed Brown use petscksp 371c4762a1bSJed Brown implicit none 372c4762a1bSJed Brown 373c4762a1bSJed Brown! Common blocks: 374c4762a1bSJed Brown PetscReal lambda 375c4762a1bSJed Brown PetscInt mx,my 376c4762a1bSJed Brown PetscBool fd_coloring 377c4762a1bSJed Brown common /params/ lambda,mx,my,fd_coloring 378c4762a1bSJed Brown 379c4762a1bSJed Brown! Input/output variables: 380c4762a1bSJed Brown PetscScalar x(mx,my) 381c4762a1bSJed Brown PetscErrorCode ierr 382c4762a1bSJed Brown 383c4762a1bSJed Brown! Local variables: 384c4762a1bSJed Brown PetscInt i,j 385c4762a1bSJed Brown PetscReal temp1,temp,hx,hy,one 386c4762a1bSJed Brown 387c4762a1bSJed Brown! Set parameters 388c4762a1bSJed Brown 389c4762a1bSJed Brown ierr = 0 390c4762a1bSJed Brown one = 1.0 391c4762a1bSJed Brown hx = one/(mx-1) 392c4762a1bSJed Brown hy = one/(my-1) 393c4762a1bSJed Brown temp1 = lambda/(lambda + one) 394c4762a1bSJed Brown 395c4762a1bSJed Brown do 20 j=1,my 396c4762a1bSJed Brown temp = min(j-1,my-j)*hy 397c4762a1bSJed Brown do 10 i=1,mx 398c4762a1bSJed Brown if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 399c4762a1bSJed Brown x(i,j) = 0.0 400c4762a1bSJed Brown else 401c4762a1bSJed Brown x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,temp)) 402c4762a1bSJed Brown endif 403c4762a1bSJed Brown 10 continue 404c4762a1bSJed Brown 20 continue 405c4762a1bSJed Brown 406c4762a1bSJed Brown return 407c4762a1bSJed Brown end 408c4762a1bSJed Brown 409c4762a1bSJed Brown! --------------------------------------------------------------------- 410c4762a1bSJed Brown! 411c4762a1bSJed Brown! FormFunction - Evaluates nonlinear function, F(x). 412c4762a1bSJed Brown! 413c4762a1bSJed Brown! Input Parameters: 414c4762a1bSJed Brown! snes - the SNES context 415c4762a1bSJed Brown! X - input vector 416c4762a1bSJed Brown! dummy - optional user-defined context, as set by SNESSetFunction() 417c4762a1bSJed Brown! (not used here) 418c4762a1bSJed Brown! 419c4762a1bSJed Brown! Output Parameter: 420c4762a1bSJed Brown! F - vector with newly computed function 421c4762a1bSJed Brown! 422c4762a1bSJed Brown! Notes: 423c4762a1bSJed Brown! This routine serves as a wrapper for the lower-level routine 424c4762a1bSJed Brown! "ApplicationFunction", where the actual computations are 425c4762a1bSJed Brown! done using the standard Fortran style of treating the local 426c4762a1bSJed Brown! vector data as a multidimensional array over the local mesh. 427c4762a1bSJed Brown! This routine merely accesses the local vector data via 428c4762a1bSJed Brown! VecGetArray() and VecRestoreArray(). 429c4762a1bSJed Brown! 430c4762a1bSJed Brown subroutine FormFunction(snes,X,F,fdcoloring,ierr) 431c4762a1bSJed Brown use petscsnes 432c4762a1bSJed Brown implicit none 433c4762a1bSJed Brown 434c4762a1bSJed Brown! Input/output variables: 435c4762a1bSJed Brown SNES snes 436c4762a1bSJed Brown Vec X,F 437c4762a1bSJed Brown PetscErrorCode ierr 438c4762a1bSJed Brown MatFDColoring fdcoloring 439c4762a1bSJed Brown 440c4762a1bSJed Brown! Common blocks: 441c4762a1bSJed Brown PetscReal lambda 442c4762a1bSJed Brown PetscInt mx,my 443c4762a1bSJed Brown PetscBool fd_coloring 444c4762a1bSJed Brown common /params/ lambda,mx,my,fd_coloring 445c4762a1bSJed Brown 446c4762a1bSJed Brown! Declarations for use with local arrays: 447c4762a1bSJed Brown PetscScalar lx_v(0:1),lf_v(0:1) 448c4762a1bSJed Brown PetscOffset lx_i,lf_i 449c4762a1bSJed Brown 450c4762a1bSJed Brown PetscInt, pointer :: indices(:) 451c4762a1bSJed Brown 452c4762a1bSJed Brown! Get pointers to vector data. 453c4762a1bSJed Brown! - For default PETSc vectors, VecGetArray() returns a pointer to 454c4762a1bSJed Brown! the data array. Otherwise, the routine is implementation dependent. 455c4762a1bSJed Brown! - You MUST call VecRestoreArray() when you no longer need access to 456c4762a1bSJed Brown! the array. 457c4762a1bSJed Brown! - Note that the Fortran interface to VecGetArray() differs from the 458c4762a1bSJed Brown! C version. See the Fortran chapter of the users manual for details. 459c4762a1bSJed Brown 460c4762a1bSJed Brown call VecGetArrayRead(X,lx_v,lx_i,ierr) 461c4762a1bSJed Brown call VecGetArray(F,lf_v,lf_i,ierr) 462c4762a1bSJed Brown 463c4762a1bSJed Brown! Compute function 464c4762a1bSJed Brown 465c4762a1bSJed Brown call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr) 466c4762a1bSJed Brown 467c4762a1bSJed Brown! Restore vectors 468c4762a1bSJed Brown 469c4762a1bSJed Brown call VecRestoreArrayRead(X,lx_v,lx_i,ierr) 470c4762a1bSJed Brown call VecRestoreArray(F,lf_v,lf_i,ierr) 471c4762a1bSJed Brown 472c4762a1bSJed Brown call PetscLogFlops(11.0d0*mx*my,ierr) 473c4762a1bSJed Brown! 474c4762a1bSJed Brown! fdcoloring is in the common block and used here ONLY to test the 475c4762a1bSJed Brown! calls to MatFDColoringGetPerturbedColumnsF90() and MatFDColoringRestorePerturbedColumnsF90() 476c4762a1bSJed Brown! 477c4762a1bSJed Brown if (fd_coloring) then 478c4762a1bSJed Brown call MatFDColoringGetPerturbedColumnsF90(fdcoloring,indices,ierr) 479c4762a1bSJed Brown print*,'Indices from GetPerturbedColumnsF90' 480c4762a1bSJed Brown write(*,1000) indices 481c4762a1bSJed Brown 1000 format(50i4) 482c4762a1bSJed Brown call MatFDColoringRestorePerturbedColumnsF90(fdcoloring,indices,ierr) 483c4762a1bSJed Brown endif 484c4762a1bSJed Brown return 485c4762a1bSJed Brown end 486c4762a1bSJed Brown 487c4762a1bSJed Brown! --------------------------------------------------------------------- 488c4762a1bSJed Brown! 489c4762a1bSJed Brown! ApplicationFunction - Computes nonlinear function, called by 490c4762a1bSJed Brown! the higher level routine FormFunction(). 491c4762a1bSJed Brown! 492c4762a1bSJed Brown! Input Parameter: 493c4762a1bSJed Brown! x - local vector data 494c4762a1bSJed Brown! 495c4762a1bSJed Brown! Output Parameters: 496c4762a1bSJed Brown! f - local vector data, f(x) 497c4762a1bSJed Brown! ierr - error code 498c4762a1bSJed Brown! 499c4762a1bSJed Brown! Notes: 500c4762a1bSJed Brown! This routine uses standard Fortran-style computations over a 2-dim array. 501c4762a1bSJed Brown! 502c4762a1bSJed Brown subroutine ApplicationFunction(x,f,ierr) 503c4762a1bSJed Brown use petscsnes 504c4762a1bSJed Brown implicit none 505c4762a1bSJed Brown 506c4762a1bSJed Brown! Common blocks: 507c4762a1bSJed Brown PetscReal lambda 508c4762a1bSJed Brown PetscInt mx,my 509c4762a1bSJed Brown PetscBool fd_coloring 510c4762a1bSJed Brown common /params/ lambda,mx,my,fd_coloring 511c4762a1bSJed Brown 512c4762a1bSJed Brown! Input/output variables: 513c4762a1bSJed Brown PetscScalar x(mx,my),f(mx,my) 514c4762a1bSJed Brown PetscErrorCode ierr 515c4762a1bSJed Brown 516c4762a1bSJed Brown! Local variables: 517c4762a1bSJed Brown PetscScalar two,one,hx,hy 518c4762a1bSJed Brown PetscScalar hxdhy,hydhx,sc 519c4762a1bSJed Brown PetscScalar u,uxx,uyy 520c4762a1bSJed Brown PetscInt i,j 521c4762a1bSJed Brown 522c4762a1bSJed Brown ierr = 0 523c4762a1bSJed Brown one = 1.0 524c4762a1bSJed Brown two = 2.0 525c4762a1bSJed Brown hx = one/(mx-1) 526c4762a1bSJed Brown hy = one/(my-1) 527c4762a1bSJed Brown sc = hx*hy*lambda 528c4762a1bSJed Brown hxdhy = hx/hy 529c4762a1bSJed Brown hydhx = hy/hx 530c4762a1bSJed Brown 531c4762a1bSJed Brown! Compute function 532c4762a1bSJed Brown 533c4762a1bSJed Brown do 20 j=1,my 534c4762a1bSJed Brown do 10 i=1,mx 535c4762a1bSJed Brown if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 536c4762a1bSJed Brown f(i,j) = x(i,j) 537c4762a1bSJed Brown else 538c4762a1bSJed Brown u = x(i,j) 539c4762a1bSJed Brown uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j)) 540c4762a1bSJed Brown uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1)) 541c4762a1bSJed Brown f(i,j) = uxx + uyy - sc*exp(u) 542c4762a1bSJed Brown endif 543c4762a1bSJed Brown 10 continue 544c4762a1bSJed Brown 20 continue 545c4762a1bSJed Brown 546c4762a1bSJed Brown return 547c4762a1bSJed Brown end 548c4762a1bSJed Brown 549c4762a1bSJed Brown! --------------------------------------------------------------------- 550c4762a1bSJed Brown! 551c4762a1bSJed Brown! FormJacobian - Evaluates Jacobian matrix. 552c4762a1bSJed Brown! 553c4762a1bSJed Brown! Input Parameters: 554c4762a1bSJed Brown! snes - the SNES context 555c4762a1bSJed Brown! x - input vector 556c4762a1bSJed Brown! dummy - optional user-defined context, as set by SNESSetJacobian() 557c4762a1bSJed Brown! (not used here) 558c4762a1bSJed Brown! 559c4762a1bSJed Brown! Output Parameters: 560c4762a1bSJed Brown! jac - Jacobian matrix 561c4762a1bSJed Brown! jac_prec - optionally different preconditioning matrix (not used here) 562c4762a1bSJed Brown! flag - flag indicating matrix structure 563c4762a1bSJed Brown! 564c4762a1bSJed Brown! Notes: 565c4762a1bSJed Brown! This routine serves as a wrapper for the lower-level routine 566c4762a1bSJed Brown! "ApplicationJacobian", where the actual computations are 567c4762a1bSJed Brown! done using the standard Fortran style of treating the local 568c4762a1bSJed Brown! vector data as a multidimensional array over the local mesh. 569c4762a1bSJed Brown! This routine merely accesses the local vector data via 570c4762a1bSJed Brown! VecGetArray() and VecRestoreArray(). 571c4762a1bSJed Brown! 572c4762a1bSJed Brown subroutine FormJacobian(snes,X,jac,jac_prec,dummy,ierr) 573c4762a1bSJed Brown use petscsnes 574c4762a1bSJed Brown implicit none 575c4762a1bSJed Brown 576c4762a1bSJed Brown! Input/output variables: 577c4762a1bSJed Brown SNES snes 578c4762a1bSJed Brown Vec X 579c4762a1bSJed Brown Mat jac,jac_prec 580c4762a1bSJed Brown PetscErrorCode ierr 581c4762a1bSJed Brown integer dummy 582c4762a1bSJed Brown 583c4762a1bSJed Brown! Common blocks: 584c4762a1bSJed Brown PetscReal lambda 585c4762a1bSJed Brown PetscInt mx,my 586c4762a1bSJed Brown PetscBool fd_coloring 587c4762a1bSJed Brown common /params/ lambda,mx,my,fd_coloring 588c4762a1bSJed Brown 589c4762a1bSJed Brown! Declarations for use with local array: 590c4762a1bSJed Brown PetscScalar lx_v(0:1) 591c4762a1bSJed Brown PetscOffset lx_i 592c4762a1bSJed Brown 593c4762a1bSJed Brown! Get a pointer to vector data 594c4762a1bSJed Brown 595c4762a1bSJed Brown call VecGetArrayRead(X,lx_v,lx_i,ierr) 596c4762a1bSJed Brown 597c4762a1bSJed Brown! Compute Jacobian entries 598c4762a1bSJed Brown 599c4762a1bSJed Brown call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr) 600c4762a1bSJed Brown 601c4762a1bSJed Brown! Restore vector 602c4762a1bSJed Brown 603c4762a1bSJed Brown call VecRestoreArrayRead(X,lx_v,lx_i,ierr) 604c4762a1bSJed Brown 605c4762a1bSJed Brown! Assemble matrix 606c4762a1bSJed Brown 607c4762a1bSJed Brown call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr) 608c4762a1bSJed Brown call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr) 609c4762a1bSJed Brown 610c4762a1bSJed Brown 611c4762a1bSJed Brown return 612c4762a1bSJed Brown end 613c4762a1bSJed Brown 614c4762a1bSJed Brown! --------------------------------------------------------------------- 615c4762a1bSJed Brown! 616c4762a1bSJed Brown! ApplicationJacobian - Computes Jacobian matrix, called by 617c4762a1bSJed Brown! the higher level routine FormJacobian(). 618c4762a1bSJed Brown! 619c4762a1bSJed Brown! Input Parameters: 620c4762a1bSJed Brown! x - local vector data 621c4762a1bSJed Brown! 622c4762a1bSJed Brown! Output Parameters: 623c4762a1bSJed Brown! jac - Jacobian matrix 624c4762a1bSJed Brown! jac_prec - optionally different preconditioning matrix (not used here) 625c4762a1bSJed Brown! ierr - error code 626c4762a1bSJed Brown! 627c4762a1bSJed Brown! Notes: 628c4762a1bSJed Brown! This routine uses standard Fortran-style computations over a 2-dim array. 629c4762a1bSJed Brown! 630c4762a1bSJed Brown subroutine ApplicationJacobian(x,jac,jac_prec,ierr) 631c4762a1bSJed Brown use petscsnes 632c4762a1bSJed Brown implicit none 633c4762a1bSJed Brown 634c4762a1bSJed Brown! Common blocks: 635c4762a1bSJed Brown PetscReal lambda 636c4762a1bSJed Brown PetscInt mx,my 637c4762a1bSJed Brown PetscBool fd_coloring 638c4762a1bSJed Brown common /params/ lambda,mx,my,fd_coloring 639c4762a1bSJed Brown 640c4762a1bSJed Brown! Input/output variables: 641c4762a1bSJed Brown PetscScalar x(mx,my) 642c4762a1bSJed Brown Mat jac,jac_prec 643c4762a1bSJed Brown PetscErrorCode ierr 644c4762a1bSJed Brown 645c4762a1bSJed Brown! Local variables: 646c4762a1bSJed Brown PetscInt i,j,row(1),col(5),i1,i5 647c4762a1bSJed Brown PetscScalar two,one, hx,hy 648c4762a1bSJed Brown PetscScalar hxdhy,hydhx,sc,v(5) 649c4762a1bSJed Brown 650c4762a1bSJed Brown! Set parameters 651c4762a1bSJed Brown 652c4762a1bSJed Brown i1 = 1 653c4762a1bSJed Brown i5 = 5 654c4762a1bSJed Brown one = 1.0 655c4762a1bSJed Brown two = 2.0 656c4762a1bSJed Brown hx = one/(mx-1) 657c4762a1bSJed Brown hy = one/(my-1) 658c4762a1bSJed Brown sc = hx*hy 659c4762a1bSJed Brown hxdhy = hx/hy 660c4762a1bSJed Brown hydhx = hy/hx 661c4762a1bSJed Brown 662c4762a1bSJed Brown! Compute entries of the Jacobian matrix 663c4762a1bSJed Brown! - Here, we set all entries for a particular row at once. 664c4762a1bSJed Brown! - Note that MatSetValues() uses 0-based row and column numbers 665c4762a1bSJed Brown! in Fortran as well as in C. 666c4762a1bSJed Brown 667c4762a1bSJed Brown do 20 j=1,my 668c4762a1bSJed Brown row(1) = (j-1)*mx - 1 669c4762a1bSJed Brown do 10 i=1,mx 670c4762a1bSJed Brown row(1) = row(1) + 1 671c4762a1bSJed Brown! boundary points 672c4762a1bSJed Brown if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then 673c4762a1bSJed Brown call MatSetValues(jac_prec,i1,row,i1,row,one,INSERT_VALUES,ierr) 674c4762a1bSJed Brown! interior grid points 675c4762a1bSJed Brown else 676c4762a1bSJed Brown v(1) = -hxdhy 677c4762a1bSJed Brown v(2) = -hydhx 678c4762a1bSJed Brown v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j)) 679c4762a1bSJed Brown v(4) = -hydhx 680c4762a1bSJed Brown v(5) = -hxdhy 681c4762a1bSJed Brown col(1) = row(1) - mx 682c4762a1bSJed Brown col(2) = row(1) - 1 683c4762a1bSJed Brown col(3) = row(1) 684c4762a1bSJed Brown col(4) = row(1) + 1 685c4762a1bSJed Brown col(5) = row(1) + mx 686c4762a1bSJed Brown call MatSetValues(jac_prec,i1,row,i5,col,v,INSERT_VALUES,ierr) 687c4762a1bSJed Brown endif 688c4762a1bSJed Brown 10 continue 689c4762a1bSJed Brown 20 continue 690c4762a1bSJed Brown 691c4762a1bSJed Brown return 692c4762a1bSJed Brown end 693c4762a1bSJed Brown 694c4762a1bSJed Brown! 695c4762a1bSJed Brown!/*TEST 696c4762a1bSJed Brown! 697c4762a1bSJed Brown! build: 698c4762a1bSJed Brown! requires: !single 699c4762a1bSJed Brown! 700c4762a1bSJed Brown! test: 701c4762a1bSJed Brown! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always 702c4762a1bSJed Brown! 703c4762a1bSJed Brown! test: 704c4762a1bSJed Brown! suffix: 2 705c4762a1bSJed Brown! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always 706c4762a1bSJed Brown! 707c4762a1bSJed Brown! test: 708c4762a1bSJed Brown! suffix: 3 709c4762a1bSJed Brown! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always 710c4762a1bSJed Brown! filter: sort -b 711c4762a1bSJed Brown! filter_output: sort -b 712c4762a1bSJed Brown! 713c4762a1bSJed Brown! test: 714c4762a1bSJed Brown! suffix: 4 715c4762a1bSJed Brown! args: -pc -par 6.807 -nox 716c4762a1bSJed Brown! 717c4762a1bSJed Brown!TEST*/ 718