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