Actual source code: ex73f90t.F90

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1: !
  2: !  Description: Solves a nonlinear system in parallel with SNES.
  3: !  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
  4: !  domain, using distributed arrays (DMDAs) to partition the parallel grid.
  5: !  The command line options include:
  6: !    -par <parameter>, where <parameter> indicates the nonlinearity of the problem
  7: !       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
  8: !
  9: !  This system (A) is augmented with constraints:
 10: !
 11: !    A -B   *  phi  =  rho
 12: !   -C  I      lam  = 0
 13: !
 14: !  where I is the identity, A is the "normal" Poisson equation, B is the "distributor" of the
 15: !  total flux (the first block equation is the flux surface averaging equation).  The second
 16: !  equation  lambda = C * x enforces the surface flux auxiliary equation.  B and C have all
 17: !  positive entries, areas in C and fraction of area in B.
 18: !
 19: !!/*T
 20: !  Concepts: SNES^parallel Bratu example
 21: !  Concepts: MatNest
 22: !  Processors: n
 23: !T*/


 26: !
 27: !  --------------------------------------------------------------------------
 28: !
 29: !  Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 30: !  the partial differential equation
 31: !
 32: !          -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 33: !
 34: !  with boundary conditions
 35: !
 36: !           u = 0  for  x = 0, x = 1, y = 0, y = 1.
 37: !
 38: !  A finite difference approximation with the usual 5-point stencil
 39: !  is used to discretize the boundary value problem to obtain a nonlinear
 40: !  system of equations.
 41: !
 42: !  --------------------------------------------------------------------------
 43: !  The following define must be used before including any PETSc include files
 44: !  into a module or interface. This is because they can't handle declarations
 45: !  in them
 46: !
 47:       module petsc_kkt_solver
 48:  #include <petsc/finclude/petscdm.h>
 49:  #include <petsc/finclude/petscmat.h>
 50:       use petscdm
 51:       use petscmat
 52:       type petsc_kkt_solver_type
 53:         DM::da
 54: !     temp A block stuff
 55:         PetscInt mx,my
 56:         PetscMPIInt rank
 57:         PetscReal lambda
 58: !     Mats
 59:         Mat::Amat,AmatLin,Bmat,CMat,Dmat
 60:         IS::isPhi,isLambda
 61:       end type petsc_kkt_solver_type

 63:       end module petsc_kkt_solver

 65:       module petsc_kkt_solver_interfaces
 66:         use petsc_kkt_solver

 68:       Interface SNESSetApplicationContext
 69:         Subroutine SNESSetApplicationContext(snesIn,ctx,ierr)
 70:  #include <petsc/finclude/petscsnes.h>
 71:         use petscsnes
 72:         use petsc_kkt_solver
 73:           SNES::    snesIn
 74:           type(petsc_kkt_solver_type) ctx
 75:           PetscErrorCode ierr
 76:         End Subroutine
 77:       End Interface SNESSetApplicationContext

 79:       Interface SNESGetApplicationContext
 80:         Subroutine SNESGetApplicationContext(snesIn,ctx,ierr)
 81:  #include <petsc/finclude/petscsnes.h>
 82:         use petscsnes
 83:         use petsc_kkt_solver
 84:           SNES::     snesIn
 85:           type(petsc_kkt_solver_type), pointer :: ctx
 86:           PetscErrorCode ierr
 87:         End Subroutine
 88:       End Interface SNESGetApplicationContext
 89:       end module petsc_kkt_solver_interfaces

 91:       program main
 92:  #include <petsc/finclude/petscdm.h>
 93:  #include <petsc/finclude/petscsnes.h>
 94:       use petscdm
 95:       use petscdmda
 96:       use petscsnes
 97:       use petsc_kkt_solver
 98:       use petsc_kkt_solver_interfaces
 99:       implicit none
100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: !                   Variable declarations
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: !
104: !  Variables:
105: !     mysnes      - nonlinear solver
106: !     x, r        - solution, residual vectors
107: !     J           - Jacobian matrix
108: !     its         - iterations for convergence
109: !     Nx, Ny      - number of preocessors in x- and y- directions
110: !
111:       SNES::       mysnes
112:       Vec::        x,r,x2,x1,x1loc,x2loc
113:       Mat::       Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4)
114: !      Mat::       tmat
115:       DM::       daphi,dalam
116:       IS::        isglobal(2)
117:       PetscErrorCode   ierr
118:       PetscInt         its,N1,N2,i,j,irow,row(1)
119:       PetscInt         col(1),low,high,lamlow,lamhigh
120:       PetscBool        flg
121:       PetscInt         ione,nfour,itwo,nloc,nloclam
122:       PetscReal lambda_max,lambda_min
123:       type(petsc_kkt_solver_type)  solver
124:       PetscScalar      bval(1),cval(1),one

126: !  Note: Any user-defined Fortran routines (such as FormJacobian)
127: !  MUST be declared as external.
128:       external FormInitialGuess,FormJacobian,FormFunction

130: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: !  Initialize program
132: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
134:      if (ierr .ne. 0) then
135:          print*,'PetscInitialize failed'
136:          stop
137:       endif
138:       call MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr);CHKERRA(ierr)

140: !  Initialize problem parameters
141:       lambda_max  = 6.81_PETSC_REAL_KIND
142:       lambda_min  = 0.0
143:       solver%lambda = 6.0
144:       ione = 1
145:       nfour = 4
146:       itwo = 2
147:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr);CHKERRA(ierr)
148:       if (solver%lambda .ge. lambda_max .or. solver%lambda .lt. lambda_min) then; SETERRA(PETSC_COMM_SELF,1,'Lambda provided with -par is out of range'); endif

150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: !  Create vector data structures; set function evaluation routine
152: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

154: !     just get size
155:       call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione, &
156:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,daphi,ierr);CHKERRA(ierr)
157:       call DMSetFromOptions(daphi,ierr);CHKERRA(ierr)
158:       call DMSetUp(daphi,ierr);CHKERRA(ierr)
159:       call DMDAGetInfo(daphi,PETSC_NULL_INTEGER,solver%mx,solver%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,                        &
160:      &                 PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
161:       N1 = solver%my*solver%mx
162:       N2 = solver%my
163:       flg = .false.
164:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr);CHKERRA(ierr)
165:       if (flg) then
166:          N2 = 0
167:       endif

169:       call DMDestroy(daphi,ierr);CHKERRA(ierr)

171: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: !  Create matrix data structure; set Jacobian evaluation routine
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174:       call DMShellCreate(PETSC_COMM_WORLD,daphi,ierr);CHKERRA(ierr)
175:       call DMSetOptionsPrefix(daphi,'phi_',ierr);CHKERRA(ierr)
176:       call DMSetFromOptions(daphi,ierr);CHKERRA(ierr)

178:       call VecCreate(PETSC_COMM_WORLD,x1,ierr);CHKERRA(ierr)
179:       call VecSetSizes(x1,PETSC_DECIDE,N1,ierr);CHKERRA(ierr)
180:       call VecSetFromOptions(x1,ierr);CHKERRA(ierr)

182:       call VecGetOwnershipRange(x1,low,high,ierr);CHKERRA(ierr)
183:       nloc = high - low

185:       call MatCreate(PETSC_COMM_WORLD,Amat,ierr);CHKERRA(ierr)
186:       call MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr)
187:       call MatSetUp(Amat,ierr);CHKERRA(ierr)

189:       call MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr);CHKERRA(ierr)
190:       call MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr)
191:       call MatSetUp(solver%AmatLin,ierr);CHKERRA(ierr)

193:       call FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr);CHKERRA(ierr)
194:       call MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
195:       call MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

197:       call DMShellSetGlobalVector(daphi,x1,ierr);CHKERRA(ierr)
198:       call DMShellSetMatrix(daphi,Amat,ierr);CHKERRA(ierr)

200:       call VecCreate(PETSC_COMM_SELF,x1loc,ierr);CHKERRA(ierr)
201:       call VecSetSizes(x1loc,nloc,nloc,ierr);CHKERRA(ierr)
202:       call VecSetFromOptions(x1loc,ierr);CHKERRA(ierr)
203:       call DMShellSetLocalVector(daphi,x1loc,ierr);CHKERRA(ierr)

205: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206: !  Create B, C, & D matrices
207: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208:       call MatCreate(PETSC_COMM_WORLD,Cmat,ierr);CHKERRA(ierr)
209:       call MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr);CHKERRA(ierr)
210:       call MatSetUp(Cmat,ierr);CHKERRA(ierr)
211: !      create data for C and B
212:       call MatCreate(PETSC_COMM_WORLD,Bmat,ierr);CHKERRA(ierr)
213:       call MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr);CHKERRA(ierr)
214:       call MatSetUp(Bmat,ierr);CHKERRA(ierr)
215: !     create data for D
216:       call MatCreate(PETSC_COMM_WORLD,Dmat,ierr);CHKERRA(ierr)
217:       call MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr);CHKERRA(ierr)
218:       call MatSetUp(Dmat,ierr);CHKERRA(ierr)

220:       call VecCreate(PETSC_COMM_WORLD,x2,ierr);CHKERRA(ierr)
221:       call VecSetSizes(x2,PETSC_DECIDE,N2,ierr);CHKERRA(ierr)
222:       call VecSetFromOptions(x2,ierr);CHKERRA(ierr)

224:       call VecGetOwnershipRange(x2,lamlow,lamhigh,ierr);CHKERRA(ierr)
225:       nloclam = lamhigh-lamlow

227: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228: !  Set fake B and C
229: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230:       one    = 1.0
231:       if( N2 .gt. 0 ) then
232:          bval(1) = -one/(solver%mx-2)
233: !     cval = -one/(solver%my*solver%mx)
234:          cval(1) = -one
235:          do 20 irow=low,high-1
236:             j = irow/solver%mx   ! row in domain
237:             i = mod(irow,solver%mx)
238:             row(1) = irow
239:             col(1) = j
240:             if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1 ) then
241:                !     no op
242:             else
243:                call MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr);CHKERRA(ierr)
244:             endif
245:             row(1) = j
246:             call MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr)
247:  20   continue
248:       endif
249:       call MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
250:       call MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
251:       call MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
252:       call MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

254: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: !  Set D (indentity)
256: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257:       do 30 j=lamlow,lamhigh-1
258:          row(1) = j
259:          cval(1) = one
260:          call MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr)
261:  30   continue
262:       call MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
263:       call MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266: !  DM for lambda (dalam) : temp driver for A block, setup A block solver data
267: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268:       call DMShellCreate(PETSC_COMM_WORLD,dalam,ierr);CHKERRA(ierr)
269:       call DMShellSetGlobalVector(dalam,x2,ierr);CHKERRA(ierr)
270:       call DMShellSetMatrix(dalam,Dmat,ierr);CHKERRA(ierr)

272:       call VecCreate(PETSC_COMM_SELF,x2loc,ierr);CHKERRA(ierr)
273:       call VecSetSizes(x2loc,nloclam,nloclam,ierr);CHKERRA(ierr)
274:       call VecSetFromOptions(x2loc,ierr);CHKERRA(ierr)
275:       call DMShellSetLocalVector(dalam,x2loc,ierr);CHKERRA(ierr)

277:       call DMSetOptionsPrefix(dalam,'lambda_',ierr);CHKERRA(ierr)
278:       call DMSetFromOptions(dalam,ierr);CHKERRA(ierr)
279: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280: !  Create field split DA
281: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
282:       call DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr);CHKERRA(ierr)
283:       call DMCompositeAddDM(solver%da,daphi,ierr);CHKERRA(ierr)
284:       call DMCompositeAddDM(solver%da,dalam,ierr);CHKERRA(ierr)
285:       call DMSetFromOptions(solver%da,ierr);CHKERRA(ierr)
286:       call DMSetUp(solver%da,ierr);CHKERRA(ierr)
287:       call DMCompositeGetGlobalISs(solver%da,isglobal,ierr);CHKERRA(ierr)
288:       solver%isPhi = isglobal(1)
289:       solver%isLambda = isglobal(2)

291: !     cache matrices
292:       solver%Amat = Amat
293:       solver%Bmat = Bmat
294:       solver%Cmat = Cmat
295:       solver%Dmat = Dmat

297:       matArray(1) = Amat
298:       matArray(2) = Bmat
299:       matArray(3) = Cmat
300:       matArray(4) = Dmat

302:       call MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr);CHKERRA(ierr)
303:       call MatSetFromOptions(KKTmat,ierr);CHKERRA(ierr)

305: !  Extract global and local vectors from DMDA; then duplicate for remaining
306: !     vectors that are the same types
307:       x = tVec(0)
308:       call MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
309:       call VecDuplicate(x,r,ierr);CHKERRA(ierr)

311:       call SNESCreate(PETSC_COMM_WORLD,mysnes,ierr);CHKERRA(ierr)

313:       call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)

315:       call SNESSetApplicationContext(mysnes,solver,ierr);CHKERRA(ierr)

317:       call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr)

319: !  Set function evaluation routine and vector
320:       call SNESSetFunction(mysnes,r,FormFunction,solver,ierr);CHKERRA(ierr)

322:       call SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr);CHKERRA(ierr)

324: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
325: !  Customize nonlinear solver; set runtime options
326: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
327: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
328:       call SNESSetFromOptions(mysnes,ierr);CHKERRA(ierr)

330: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
331: !  Evaluate initial guess; then solve nonlinear system.
332: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333: !  Note: The user should initialize the vector, x, with the initial guess
334: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
335: !  to employ an initial guess of zero, the user should explicitly set
336: !  this vector to zero by calling VecSet().

338:       call FormInitialGuess(mysnes,x,ierr);CHKERRA(ierr)
339:       call SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr);CHKERRA(ierr)
340:       call SNESGetIterationNumber(mysnes,its,ierr);CHKERRA(ierr)
341:       if (solver%rank .eq. 0) then
342:          write(6,100) its
343:       endif
344:   100 format('Number of SNES iterations = ',i5)

346: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
347: !  Free work space.  All PETSc objects should be destroyed when they
348: !  are no longer needed.
349: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
350:       call MatDestroy(KKTmat,ierr);CHKERRA(ierr)
351:       call MatDestroy(Amat,ierr);CHKERRA(ierr)
352:       call MatDestroy(Dmat,ierr);CHKERRA(ierr)
353:       call MatDestroy(Bmat,ierr);CHKERRA(ierr)
354:       call MatDestroy(Cmat,ierr);CHKERRA(ierr)
355:       call MatDestroy(solver%AmatLin,ierr);CHKERRA(ierr)
356:       call ISDestroy(solver%isPhi,ierr);CHKERRA(ierr)
357:       call ISDestroy(solver%isLambda,ierr);CHKERRA(ierr)
358:       call VecDestroy(x,ierr);CHKERRA(ierr)
359:       call VecDestroy(x2,ierr);CHKERRA(ierr)
360:       call VecDestroy(x1,ierr);CHKERRA(ierr)
361:       call VecDestroy(x1loc,ierr);CHKERRA(ierr)
362:       call VecDestroy(x2loc,ierr);CHKERRA(ierr)
363:       call VecDestroy(r,ierr);CHKERRA(ierr)
364:       call SNESDestroy(mysnes,ierr);CHKERRA(ierr)
365:       call DMDestroy(solver%da,ierr);CHKERRA(ierr)
366:       call DMDestroy(daphi,ierr);CHKERRA(ierr)
367:       call DMDestroy(dalam,ierr);CHKERRA(ierr)

369:       call PetscFinalize(ierr)
370:       end

372: ! ---------------------------------------------------------------------
373: !
374: !  FormInitialGuess - Forms initial approximation.
375: !
376: !  Input Parameters:
377: !  X - vector
378: !
379: !  Output Parameter:
380: !  X - vector
381: !
382: !  Notes:
383: !  This routine serves as a wrapper for the lower-level routine
384: !  "InitialGuessLocal", where the actual computations are
385: !  done using the standard Fortran style of treating the local
386: !  vector data as a multidimensional array over the local mesh.
387: !  This routine merely handles ghost point scatters and accesses
388: !  the local vector data via VecGetArrayF90() and VecRestoreArrayF90().
389: !
390:       subroutine FormInitialGuess(mysnes,Xnest,ierr)
391:  #include <petsc/finclude/petscsnes.h>
392:       use petscsnes
393:       use petsc_kkt_solver
394:       use petsc_kkt_solver_interfaces
395:       implicit none
396: !  Input/output variables:
397:       SNES::     mysnes
398:       Vec::      Xnest
399:       PetscErrorCode ierr

401: !  Declarations for use with local arrays:
402:       type(petsc_kkt_solver_type), pointer:: solver
403:       Vec::      Xsub(2)
404:       PetscInt::  izero,ione,itwo

406:       izero = 0
407:       ione = 1
408:       itwo = 2
409:       0
410:       call SNESGetApplicationContext(mysnes,solver,ierr);CHKERRQ(ierr)
411:       call DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

413:       call InitialGuessLocal(solver,Xsub(1),ierr);CHKERRQ(ierr)
414:       call VecAssemblyBegin(Xsub(1),ierr);CHKERRQ(ierr)
415:       call VecAssemblyEnd(Xsub(1),ierr);CHKERRQ(ierr)

417: !     zero out lambda
418:       call VecZeroEntries(Xsub(2),ierr);CHKERRQ(ierr)
419:       call DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

421:       return
422:       end subroutine FormInitialGuess

424: ! ---------------------------------------------------------------------
425: !
426: !  InitialGuessLocal - Computes initial approximation, called by
427: !  the higher level routine FormInitialGuess().
428: !
429: !  Input Parameter:
430: !  X1 - local vector data
431: !
432: !  Output Parameters:
433: !  x - local vector data
434: !  ierr - error code
435: !
436: !  Notes:
437: !  This routine uses standard Fortran-style computations over a 2-dim array.
438: !
439:       subroutine InitialGuessLocal(solver,X1,ierr)
440:  #include <petsc/finclude/petscsys.h>
441:       use petscsys
442:       use petsc_kkt_solver
443:       implicit none
444: !  Input/output variables:
445:       type (petsc_kkt_solver_type)         solver
446:       Vec::      X1
447:       PetscErrorCode ierr

449: !  Local variables:
450:       PetscInt      row,i,j,ione,low,high
451:       PetscReal   temp1,temp,hx,hy,v
452:       PetscReal   one

454: !  Set parameters
455:       ione = 1
456:       0
457:       one    = 1.0
458:       hx     = one/(solver%mx-1)
459:       hy     = one/(solver%my-1)
460:       temp1  = solver%lambda/(solver%lambda + one) + one

462:       call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)

464:       do 20 row=low,high-1
465:          j = row/solver%mx
466:          i = mod(row,solver%mx)
467:          temp = min(j,solver%my-j+1)*hy
468:          if (i .eq. 0 .or. j .eq. 0  .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1 ) then
469:             v = 0.0
470:          else
471:             v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp))
472:          endif
473:          call VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
474:  20   continue

476:       return
477:       end subroutine InitialGuessLocal


480: ! ---------------------------------------------------------------------
481: !
482: !  FormJacobian - Evaluates Jacobian matrix.
483: !
484: !  Input Parameters:
485: !  dummy     - the SNES context
486: !  x         - input vector
487: !  solver    - solver data
488: !
489: !  Output Parameters:
490: !  jac      - Jacobian matrix
491: !  jac_prec - optionally different preconditioning matrix (not used here)
492: !  flag     - flag indicating matrix structure
493: !
494:       subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr)
495:  #include <petsc/finclude/petscsnes.h>
496:       use petscsnes
497:       use petsc_kkt_solver
498:       implicit none
499: !  Input/output variables:
500:       SNES::     dummy
501:       Vec::      X
502:      Mat::     jac,jac_prec
503:       type(petsc_kkt_solver_type)  solver
504:       PetscErrorCode ierr

506: !  Declarations for use with local arrays:
507:       Vec::      Xsub(1)
508:      Mat::     Amat
509:       PetscInt       izero,ione

511:       izero = 0
512:       ione = 1

514:       call DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

516: !     Compute entries for the locally owned part of the Jacobian preconditioner.
517:       call MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr);CHKERRQ(ierr)

519:       call FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr);CHKERRQ(ierr)
520:       call MatDestroy(Amat,ierr);CHKERRQ(ierr) ! discard our reference
521:       call DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)

523:       ! the rest of the matrix is not touched
524:       call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
525:       call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
526:       if (jac .ne. jac_prec) then
527:          call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
528:          call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
529:       end if

531: !     Tell the matrix we will never add a new nonzero location to the
532: !     matrix. If we do it will generate an error.
533:       call MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr);CHKERRQ(ierr)

535:       return
536:       end subroutine FormJacobian
537: 
538: ! ---------------------------------------------------------------------
539: !
540: !  FormJacobianLocal - Computes Jacobian preconditioner matrix,
541: !  called by the higher level routine FormJacobian().
542: !
543: !  Input Parameters:
544: !  x        - local vector data
545: !
546: !  Output Parameters:
547: !  jac - Jacobian preconditioner matrix
548: !  ierr     - error code
549: !
550: !  Notes:
551: !  This routine uses standard Fortran-style computations over a 2-dim array.
552: !
553:       subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr)
554:  #include <petsc/finclude/petscmat.h>
555:       use petscmat
556:       use petsc_kkt_solver
557:       implicit none
558: !  Input/output variables:
559:       type (petsc_kkt_solver_type) solver
560:       Vec::      X1
561:      Mat::     jac
562:       logical        add_nl_term
563:       PetscErrorCode ierr

565: !  Local variables:
566:       PetscInt    irow,row(1),col(5),i,j
567:       PetscInt    ione,ifive,low,high,ii
568:       PetscScalar two,one,hx,hy,hy2inv
569:       PetscScalar hx2inv,sc,v(5)
570:       PetscScalar,pointer :: lx_v(:)

572: !  Set parameters
573:       ione   = 1
574:       ifive  = 5
575:       one    = 1.0
576:       two    = 2.0
577:       hx     = one/(solver%mx-1)
578:       hy     = one/(solver%my-1)
579:       sc     = solver%lambda
580:       hx2inv = one/(hx*hx)
581:       hy2inv = one/(hy*hy)

583:       call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)
584:       call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)

586:       ii = 0
587:       do 20 irow=low,high-1
588:          j = irow/solver%mx
589:          i = mod(irow,solver%mx)
590:          ii = ii + 1            ! one based local index
591: !     boundary points
592:          if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then
593:             col(1) = irow
594:             row(1) = irow
595:             v(1)   = one
596:             call MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
597: !     interior grid points
598:          else
599:             v(1) = -hy2inv
600:             if(j-1==0) v(1) = 0.0
601:             v(2) = -hx2inv
602:             if(i-1==0) v(2) = 0.0
603:             v(3) = two*(hx2inv + hy2inv)
604:             if(add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii))
605:             v(4) = -hx2inv
606:             if(i+1==solver%mx-1) v(4) = 0.0
607:             v(5) = -hy2inv
608:             if(j+1==solver%my-1) v(5) = 0.0
609:             col(1) = irow - solver%mx
610:             col(2) = irow - 1
611:             col(3) = irow
612:             col(4) = irow + 1
613:             col(5) = irow + solver%mx
614:             row(1) = irow
615:             call MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
616:          endif
617:  20   continue

619:       call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)

621:       return
622:       end subroutine FormJacobianLocal


625: ! ---------------------------------------------------------------------
626: !
627: !  FormFunction - Evaluates nonlinear function, F(x).
628: !
629: !  Input Parameters:
630: !  snes - the SNES context
631: !  X - input vector
632: !  dummy - optional user-defined context, as set by SNESSetFunction()
633: !          (not used here)
634: !
635: !  Output Parameter:
636: !  F - function vector
637: !
638:       subroutine FormFunction(snesIn,X,F,solver,ierr)
639:  #include <petsc/finclude/petscsnes.h>
640:       use petscsnes
641:       use petsc_kkt_solver
642:       implicit none
643: !  Input/output variables:
644:       SNES::     snesIn
645:      Vec::      X,F
646:       PetscErrorCode ierr
647:       type (petsc_kkt_solver_type) solver

649: !  Declarations for use with local arrays:
650:      Vec::              Xsub(2),Fsub(2)
651:       PetscInt               izero,ione,itwo

653: !  Scatter ghost points to local vector, using the 2-step process
654: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
655: !  By placing code between these two statements, computations can
656: !  be done while messages are in transition.
657: 
658:       izero = 0
659:       ione = 1
660:       itwo = 2
661:       call DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
662:       call DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)

664:       call FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr );CHKERRQ(ierr)
665:       call MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)

667: !     do rest of operator (linear)
668:       call MatMult(    solver%Cmat, Xsub(1),      Fsub(2), ierr);CHKERRQ(ierr)
669:       call MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr)
670:       call MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr);CHKERRQ(ierr)

672:       call DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr)
673:       call DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
674:       return
675:       end subroutine formfunction


678: ! ---------------------------------------------------------------------
679: !
680: !  FormFunctionNLTerm - Computes nonlinear function, called by
681: !  the higher level routine FormFunction().
682: !
683: !  Input Parameter:
684: !  x - local vector data
685: !
686: !  Output Parameters:
687: !  f - local vector data, f(x)
688: !  ierr - error code
689: !
690: !  Notes:
691: !  This routine uses standard Fortran-style computations over a 2-dim array.
692: !
693:       subroutine FormFunctionNLTerm(X1,F1,solver,ierr)
694:  #include <petsc/finclude/petscvec.h>
695:       use petscvec
696:       use petsc_kkt_solver
697:       implicit none
698: !  Input/output variables:
699:       type (petsc_kkt_solver_type) solver
700:      Vec::      X1,F1
701:       PetscErrorCode ierr
702: !  Local variables:
703:       PetscScalar one,sc
704:       PetscScalar u,v(1)
705:       PetscInt  i,j,low,high,ii,ione,irow,row(1)
706:       PetscScalar,pointer :: lx_v(:)

708:       one    = 1.0
709:       sc     = solver%lambda
710:       ione   = 1

712:       call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)
713:       call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr)

715: !     Compute function over the locally owned part of the grid
716:       ii = 0
717:       do 20 irow=low,high-1
718:          j = irow/solver%mx
719:          i = mod(irow,solver%mx)
720:          ii = ii + 1            ! one based local index
721:          row(1) = irow
722:          if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1 ) then
723:             v(1) = 0.0
724:          else
725:             u = lx_v(ii)
726:             v(1) = -sc*exp(u)
727:          endif
728:          call VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr)
729:  20   continue

731:       call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr)

733:       call VecAssemblyBegin(F1,ierr);CHKERRQ(ierr)
734:       call VecAssemblyEnd(F1,ierr);CHKERRQ(ierr)

736:       0
737:       return
738:       end subroutine FormFunctionNLTerm

740: !/*TEST
741: !
742: !   build:
743: !      requires: !single !libpgf90 !complex
744: !
745: !   test:
746: !      nsize: 4
747: !      args: -par 5.0 -da_grid_x 10 -da_grid_y 10 -snes_monitor_short -snes_linesearch_type basic -snes_converged_reason -ksp_type fgmres -ksp_norm_type unpreconditioned -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type upper -ksp_monitor_short -fieldsplit_lambda_ksp_type preonly -fieldsplit_lambda_pc_type jacobi -fieldsplit_phi_pc_type gamg -fieldsplit_phi_pc_gamg_agg_nsmooths 1 -fieldsplit_phi_pc_gamg_threshold 0. -fieldsplit_phi_gamg_est_ksp_type cg
748: !
749: !TEST*/