Actual source code: shashi.F90


  2:       program main
  3: #include <petsc/finclude/petsc.h>
  4:       use petsc
  5:       implicit none

  7: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  8: !                   Variable declarations
  9: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 10: !
 11: !  Variables:
 12: !     snes        - nonlinear solver
 13: !     ksp        - linear solver
 14: !     pc          - preconditioner context
 15: !     ksp         - Krylov subspace method context
 16: !     x, r        - solution, residual vectors
 17: !     J           - Jacobian matrix
 18: !     its         - iterations for convergence
 19: !
 20:       SNES     snes
 21:       PC       pc
 22:       KSP      ksp
 23:       Vec      x,r,lb,ub
 24:       Mat      J
 25:       SNESLineSearch linesearch
 26:       PetscErrorCode  ierr
 27:       PetscInt its,i2,i20
 28:       PetscMPIInt size
 29:       PetscScalar   pfive
 30:       PetscReal   tol
 31:       PetscBool   setls
 32:       PetscScalar,pointer :: xx(:)
 33:       PetscScalar zero,big
 34:       SNESLineSearch ls

 36: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 37: !  MUST be declared as external.

 39:       external FormFunction, FormJacobian
 40:       external ShashiPostCheck

 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !                   Macro definitions
 44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45: !
 46: !  Macros to make clearer the process of setting values in vectors and
 47: !  getting values from vectors.  These vectors are used in the routines
 48: !  FormFunction() and FormJacobian().
 49: !   - The element lx_a(ib) is element ib in the vector x
 50: !
 51: #define lx_a(ib) lx_v(lx_i + (ib))
 52: #define lf_a(ib) lf_v(lf_i + (ib))
 53: !
 54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55: !                 Beginning of program
 56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 58:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 59:       if (ierr .ne. 0) then
 60:          print*,'Unable to initialize PETSc'
 61:          stop
 62:       endif
 63:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 64:       if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,1,'requires one process'); endif

 66:       big  = 2.88
 67:       big  = PETSC_INFINITY
 68:       zero = 0.0
 69:       i2  = 26
 70: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 71: !  Create nonlinear solver context
 72: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

 74:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

 76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77: !  Create matrix and vector data structures; set corresponding routines
 78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 80: !  Create vectors for solution and nonlinear function

 82:       call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
 83:       call VecDuplicate(x,r,ierr)

 85: !  Create Jacobian matrix data structure

 87:       call MatCreateDense(PETSC_COMM_SELF,26,26,26,26,                          &
 88:      &                    PETSC_NULL_SCALAR,J,ierr)

 90: !  Set function evaluation routine and vector

 92:       call SNESSetFunction(snes,r,FormFunction,0,ierr)

 94: !  Set Jacobian matrix data structure and Jacobian evaluation routine

 96:       call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)

 98:       call VecDuplicate(x,lb,ierr)
 99:       call VecDuplicate(x,ub,ierr)
100:       call VecSet(lb,zero,ierr)
101: !      call VecGetArrayF90(lb,xx,ierr)
102: !      call ShashiLowerBound(xx)
103: !      call VecRestoreArrayF90(lb,xx,ierr)
104:       call VecSet(ub,big,ierr)
105: !      call SNESVISetVariableBounds(snes,lb,ub,ierr)

107:       call SNESGetLineSearch(snes,ls,ierr)
108:       call SNESLineSearchSetPostCheck(ls,ShashiPostCheck,                 &
109:      &                                0,ierr)
110:       call SNESSetType(snes,SNESVINEWTONRSLS,ierr)

112:       call SNESSetFromOptions(snes,ierr)

114: !     set initial guess

116:       call VecGetArrayF90(x,xx,ierr)
117:       call ShashiInitialGuess(xx)
118:       call VecRestoreArrayF90(x,xx,ierr)

120:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)

122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: !  Free work space.  All PETSc objects should be destroyed when they
124: !  are no longer needed.
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:       call VecDestroy(lb,ierr)
127:       call VecDestroy(ub,ierr)
128:       call VecDestroy(x,ierr)
129:       call VecDestroy(r,ierr)
130:       call MatDestroy(J,ierr)
131:       call SNESDestroy(snes,ierr)
132:       call PetscFinalize(ierr)
133:       end
134: !
135: ! ------------------------------------------------------------------------
136: !
137: !  FormFunction - Evaluates nonlinear function, F(x).
138: !
139: !  Input Parameters:
140: !  snes - the SNES context
141: !  x - input vector
142: !  dummy - optional user-defined context (not used here)
143: !
144: !  Output Parameter:
145: !  f - function vector
146: !
147:       subroutine FormFunction(snes,x,f,dummy,ierr)
148: #include <petsc/finclude/petscsnes.h>
149:       use petscsnes
150:       implicit none
151:       SNES     snes
152:       Vec      x,f
153:       PetscErrorCode ierr
154:       integer dummy(*)

156: !  Declarations for use with local arrays

158:       PetscScalar  lx_v(2),lf_v(2)
159:       PetscOffset  lx_i,lf_i

161: !  Get pointers to vector data.
162: !    - For default PETSc vectors, VecGetArray() returns a pointer to
163: !      the data array.  Otherwise, the routine is implementation dependent.
164: !    - You MUST call VecRestoreArray() when you no longer need access to
165: !      the array.
166: !    - Note that the Fortran interface to VecGetArray() differs from the
167: !      C version.  See the Fortran chapter of the users manual for details.

169:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
170:       call VecGetArray(f,lf_v,lf_i,ierr)
171:       call ShashiFormFunction(lx_a(1),lf_a(1))
172:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
173:       call VecRestoreArray(f,lf_v,lf_i,ierr)

175:       return
176:       end

178: ! ---------------------------------------------------------------------
179: !
180: !  FormJacobian - Evaluates Jacobian matrix.
181: !
182: !  Input Parameters:
183: !  snes - the SNES context
184: !  x - input vector
185: !  dummy - optional user-defined context (not used here)
186: !
187: !  Output Parameters:
188: !  A - Jacobian matrix
189: !  B - optionally different preconditioning matrix
190: !  flag - flag indicating matrix structure
191: !
192:       subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
193: #include <petsc/finclude/petscsnes.h>
194:       use petscsnes
195:       implicit none
196:       SNES         snes
197:       Vec          X
198:       Mat          jac,B
199:       PetscScalar  A(4)
200:       PetscErrorCode ierr
201:       PetscInt idx(2),i2
202:       integer dummy(*)

204: !  Declarations for use with local arrays

206:       PetscScalar lx_v(1),lf_v(1)
207:       PetscOffset lx_i,lf_i

209: !  Get pointer to vector data

211:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
212:       call MatDenseGetArray(B,lf_v,lf_i,ierr)
213:       call ShashiFormJacobian(lx_a(1),lf_a(1))
214:       call MatDenseRestoreArray(B,lf_v,lf_i,ierr)
215:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

217: !  Assemble matrix

219:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
220:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)

222:       return
223:       end

225:             subroutine ShashiLowerBound(an_r)
226: !        implicit PetscScalar (a-h,o-z)
227:         implicit none
228:         PetscScalar an_r(26)
229:         PetscInt i

231:         do i=2,26
232:            an_r(i) = 1000.0/6.023D+23
233:         enddo
234:         return
235:         end

237:             subroutine ShashiInitialGuess(an_r)
238: !        implicit PetscScalar (a-h,o-z)
239:         implicit none
240:         PetscScalar an_c_additive
241:         PetscScalar       an_h_additive
242:         PetscScalar an_o_additive
243:         PetscScalar   atom_c_init
244:         PetscScalar atom_h_init
245:         PetscScalar atom_n_init
246:         PetscScalar atom_o_init
247:         PetscScalar h_init
248:         PetscScalar p_init
249:         PetscInt nfuel
250:         PetscScalar temp,pt
251:         PetscScalar an_r(26),k_eq(23),f_eq(26)
252:         PetscScalar d_eq(26,26),H_molar(26)
253:         PetscInt an_h(1),an_c(1)
254:         PetscScalar part_p(26)
255:         PetscInt i_cc,i_hh,i_h2o
256:         PetscInt i_pwr2,i_co2_h2o

258:         pt = 0.1
259:         atom_c_init =6.7408177364816552D-022
260:         atom_h_init = 2.0
261:         atom_o_init = 1.0
262:         atom_n_init = 3.76
263:         nfuel = 1
264:         an_c(1) = 1
265:         an_h(1) = 4
266:         an_c_additive = 2
267:         an_h_additive = 6
268:         an_o_additive = 1
269:         h_init = 128799.7267952987
270:         p_init = 0.1
271:         temp = 1500

273:       an_r( 1) =     1.66000D-24
274:       an_r( 2) =     1.66030D-22
275:       an_r( 3) =     5.00000D-01
276:       an_r( 4) =     1.66030D-22
277:       an_r( 5) =     1.66030D-22
278:       an_r( 6) =     1.88000D+00
279:       an_r( 7) =     1.66030D-22
280:       an_r( 8) =     1.66030D-22
281:       an_r( 9) =     1.66030D-22
282:       an_r(10) =     1.66030D-22
283:       an_r(11) =     1.66030D-22
284:       an_r(12) =     1.66030D-22
285:       an_r(13) =     1.66030D-22
286:       an_r(14) =     1.00000D+00
287:       an_r(15) =     1.66030D-22
288:       an_r(16) =     1.66030D-22
289:       an_r(17) =     1.66000D-24
290:       an_r(18) =     1.66030D-24
291:       an_r(19) =     1.66030D-24
292:       an_r(20) =     1.66030D-24
293:       an_r(21) =     1.66030D-24
294:       an_r(22) =     1.66030D-24
295:       an_r(23) =     1.66030D-24
296:       an_r(24) =     1.66030D-24
297:       an_r(25) =     1.66030D-24
298:       an_r(26) =     1.66030D-24

300:       an_r = 0
301:       an_r( 3) =     .5
302:       an_r( 6) =     1.88000
303:       an_r(14) =     1.

305: #if defined(solution)
306:       an_r( 1) =      3.802208D-33
307:       an_r( 2) =      1.298287D-29
308:       an_r( 3) =      2.533067D-04
309:       an_r( 4) =      6.865078D-22
310:       an_r( 5) =      9.993125D-01
311:       an_r( 6) =      1.879964D+00
312:       an_r( 7) =      4.449489D-13
313:       an_r( 8) =      3.428687D-07
314:       an_r( 9) =      7.105138D-05
315:       an_r(10) =      1.094368D-04
316:       an_r(11) =      2.362305D-06
317:       an_r(12) =      1.107145D-09
318:       an_r(13) =      1.276162D-24
319:       an_r(14) =      6.315538D-04
320:       an_r(15) =      2.356540D-09
321:       an_r(16) =      2.048248D-09
322:       an_r(17) =      1.966187D-22
323:       an_r(18) =      7.856497D-29
324:       an_r(19) =      1.987840D-36
325:       an_r(20) =      8.182441D-22
326:       an_r(21) =      2.684880D-16
327:       an_r(22) =      2.680473D-16
328:       an_r(23) =      6.594967D-18
329:       an_r(24) =      2.509714D-21
330:       an_r(25) =      3.096459D-21
331:       an_r(26) =      6.149551D-18
332: #endif

334:       return
335:       end

337:       subroutine ShashiFormFunction(an_r,f_eq)
338: !       implicit PetscScalar (a-h,o-z)
339:         implicit none
340:         PetscScalar an_c_additive
341:         PetscScalar       an_h_additive
342:         PetscScalar an_o_additive
343:         PetscScalar   atom_c_init
344:         PetscScalar atom_h_init
345:         PetscScalar atom_n_init
346:         PetscScalar atom_o_init
347:         PetscScalar h_init
348:         PetscScalar p_init
349:         PetscInt nfuel
350:         PetscScalar temp,pt
351:        PetscScalar an_r(26),k_eq(23),f_eq(26)
352:        PetscScalar d_eq(26,26),H_molar(26)
353:        PetscInt an_h(1),an_c(1)
354:        PetscScalar part_p(26),idiff
355:         PetscInt i_cc,i_hh,i_h2o
356:         PetscInt i_pwr2,i_co2_h2o
357:         PetscScalar an_t,sum_h,pt_cubed,pt_five
358:         PetscScalar pt_four,pt_val1,pt_val2
359:         PetscScalar a_io2
360:         PetscInt i,ip
361:         pt = 0.1
362:         atom_c_init =6.7408177364816552D-022
363:         atom_h_init = 2.0
364:         atom_o_init = 1.0
365:         atom_n_init = 3.76
366:         nfuel = 1
367:         an_c(1) = 1
368:         an_h(1) = 4
369:         an_c_additive = 2
370:         an_h_additive = 6
371:         an_o_additive = 1
372:         h_init = 128799.7267952987
373:         p_init = 0.1
374:         temp = 1500

376:        k_eq( 1) =     1.75149D-05
377:        k_eq( 2) =     4.01405D-06
378:        k_eq( 3) =     6.04663D-14
379:        k_eq( 4) =     2.73612D-01
380:        k_eq( 5) =     3.25592D-03
381:        k_eq( 6) =     5.33568D+05
382:        k_eq( 7) =     2.07479D+05
383:        k_eq( 8) =     1.11841D-02
384:        k_eq( 9) =     1.72684D-03
385:        k_eq(10) =     1.98588D-07
386:        k_eq(11) =     7.23600D+27
387:        k_eq(12) =     5.73926D+49
388:        k_eq(13) =     1.00000D+00
389:        k_eq(14) =     1.64493D+16
390:        k_eq(15) =     2.73837D-29
391:        k_eq(16) =     3.27419D+50
392:        k_eq(17) =     1.72447D-23
393:        k_eq(18) =     4.24657D-06
394:        k_eq(19) =     1.16065D-14
395:        k_eq(20) =     3.28020D+25
396:        k_eq(21) =     1.06291D+00
397:        k_eq(22) =     9.11507D+02
398:        k_eq(23) =     6.02837D+03

400:        H_molar( 1) =     3.26044D+03
401:        H_molar( 2) =    -8.00407D+04
402:        H_molar( 3) =     4.05873D+04
403:        H_molar( 4) =    -3.31849D+05
404:        H_molar( 5) =    -1.93654D+05
405:        H_molar( 6) =     3.84035D+04
406:        H_molar( 7) =     4.97589D+05
407:        H_molar( 8) =     2.74483D+05
408:        H_molar( 9) =     1.30022D+05
409:        H_molar(10) =     7.58429D+04
410:        H_molar(11) =     2.42948D+05
411:        H_molar(12) =     1.44588D+05
412:        H_molar(13) =    -7.16891D+04
413:        H_molar(14) =     3.63075D+04
414:        H_molar(15) =     9.23880D+04
415:        H_molar(16) =     6.50477D+04
416:        H_molar(17) =     3.04310D+05
417:        H_molar(18) =     7.41707D+05
418:        H_molar(19) =     6.32767D+05
419:        H_molar(20) =     8.90624D+05
420:        H_molar(21) =     2.49805D+04
421:        H_molar(22) =     6.43473D+05
422:        H_molar(23) =     1.02861D+06
423:        H_molar(24) =    -6.07503D+03
424:        H_molar(25) =     1.27020D+05
425:        H_molar(26) =    -1.07011D+05
426: !=============
427:       an_t = 0
428:       sum_h = 0

430:       do i = 1,26
431:          an_t = an_t + an_r(i)
432:       enddo

434:         f_eq(1) = atom_h_init                                           &
435:      &          - (an_h(1)*an_r(1) + an_h_additive*an_r(2)              &
436:      &              + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14)      &
437:      &              + an_r(16)+ 2*an_r(17) + an_r(19)                   &
438:      &              +an_r(20) + 3*an_r(22)+an_r(26))

440:         f_eq(2) = atom_o_init                                           &
441:      &          - (an_o_additive*an_r(2) + 2*an_r(3)                    &
442:      &             + 2*an_r(4) + an_r(5)                                &
443:      &             + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
444:      &             + 2*an_r(15) + 2*an_r(16)+ an_r(20) + an_r(22)       &
445:      &             + an_r(23) + 2*an_r(24) + 1*an_r(25)+an_r(26))

447:         f_eq(3) = an_r(2)-1.0d-150

449:         f_eq(4) = atom_c_init                                           &
450:      &          - (an_c(1)*an_r(1) + an_c_additive * an_r(2)            &
451:      &          + an_r(4) + an_r(13)+ 2*an_r(17) + an_r(18)             &
452:      &          + an_r(19) + an_r(20))

454:         do ip = 1,26
455:            part_p(ip) = (an_r(ip)/an_t) * pt
456:         enddo

458:         i_cc      = an_c(1)
459:         i_hh      = an_h(1)
460:         a_io2      = i_cc + i_hh/4.0
461:         i_h2o     = i_hh/2
462:         idiff     = (i_cc + i_h2o) - (a_io2 + 1)

464:         f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2                       &
465:      &          - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
466: !           write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
467: !          stop
468:         f_eq(6) = atom_n_init                                           &
469:      &          - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12)           &
470:      &              + an_r(15)                                          &
471:      &              + an_r(23))

473:       f_eq( 7) = part_p(11)                                             &
474:      &         - (k_eq( 1) * sqrt(part_p(14)+1d-23))
475:       f_eq( 8) = part_p( 8)                                             &
476:      &         - (k_eq( 2) * sqrt(part_p( 3)+1d-23))

478:       f_eq( 9) = part_p( 7)                                             &
479:      &         - (k_eq( 3) * sqrt(part_p( 6)+1d-23))

481:       f_eq(10) = part_p(10)                                             &
482:      &         - (k_eq( 4) * sqrt(part_p( 3)+1d-23))                    &
483:      &         * sqrt(part_p(14))

485:       f_eq(11) = part_p( 9)                                             &
486:      &         - (k_eq( 5) * sqrt(part_p( 3)+1d-23))                    &
487:      &         * sqrt(part_p( 6)+1d-23)
488:       f_eq(12) = part_p( 5)                                             &
489:      &         - (k_eq( 6) * sqrt(part_p( 3)+1d-23))                    &
490:      &         * (part_p(14))

492:       f_eq(13) = part_p( 4)                                             &
493:      &         - (k_eq( 7) * sqrt(part_p(3)+1.0d-23))                   &
494:      &         * (part_p(13))

496:       f_eq(14) = part_p(15)                                             &
497:      &         - (k_eq( 8) * sqrt(part_p(3)+1.0d-50))                   &
498:      &         * (part_p( 9))
499:       f_eq(15) = part_p(16)                                             &
500:      &         - (k_eq( 9) * part_p( 3))                                &
501:      &         * sqrt(part_p(14)+1d-23)

503:       f_eq(16) = part_p(12)                                             &
504:      &         - (k_eq(10) * sqrt(part_p( 3)+1d-23))                    &
505:      &         * (part_p( 6))

507:       f_eq(17) = part_p(14)*part_p(18)**2                               &
508:      &         - (k_eq(15) * part_p(17))

510:       f_eq(18) = (part_p(13)**2)                                        &
511:      &     - (k_eq(16) * part_p(3)*part_p(18)**2)
512:       print*,f_eq(18),part_p(3),part_p(18),part_p(13),k_eq(16)

514:       f_eq(19) = part_p(19)*part_p(3) - k_eq(17)*part_p(13)*part_p(10)

516:       f_eq(20) = part_p(21)*part_p(20) - k_eq(18)*part_p(19)*part_p(8)

518:       f_eq(21) = part_p(21)*part_p(23) - k_eq(19)*part_p(7)*part_p(8)

520:       f_eq(22) = part_p(5)*part_p(11) - k_eq(20)*part_p(21)*part_p(22)

522:       f_eq(23) = part_p(24) - k_eq(21)*part_p(21)*part_p(3)

524:       f_eq(24) =  part_p(3)*part_p(25) - k_eq(22)*part_p(24)*part_p(8)

526:       f_eq(25) =  part_p(26) - k_eq(23)*part_p(21)*part_p(10)

528:       f_eq(26) = -(an_r(20) + an_r(22) + an_r(23))                      &
529:      &          +(an_r(21) + an_r(24) + an_r(25) + an_r(26))

531:              do i = 1,26
532:                  write(44,*)i,f_eq(i)
533:               enddo

535:       return
536:       end

538:       subroutine ShashiFormJacobian(an_r,d_eq)
539: !        implicit PetscScalar (a-h,o-z)
540:         implicit none
541:         PetscScalar an_c_additive
542:         PetscScalar       an_h_additive
543:         PetscScalar an_o_additive
544:         PetscScalar   atom_c_init
545:         PetscScalar atom_h_init
546:         PetscScalar atom_n_init
547:         PetscScalar atom_o_init
548:         PetscScalar h_init
549:         PetscScalar p_init
550:         PetscInt nfuel
551:         PetscScalar temp,pt
552:         PetscScalar an_t,ai_o2,sum_h
553:         PetscScalar an_tot1_d,an_tot1
554:         PetscScalar an_tot2_d,an_tot2
555:         PetscScalar const5,const3,const2
556:         PetscScalar   const_cube
557:         PetscScalar   const_five
558:         PetscScalar   const_four
559:         PetscScalar   const_six
560:         PetscInt jj,jb,ii3,id,ib,ip,i
561:         PetscScalar   pt2,pt1
562:         PetscScalar an_r(26),k_eq(23),f_eq(26)
563:         PetscScalar d_eq(26,26),H_molar(26)
564:         PetscInt an_h(1),an_c(1)
565:         PetscScalar ai_pwr1,part_p(26),idiff
566:         PetscInt i_cc,i_hh
567:         PetscInt i_h2o,i_pwr2,i_co2_h2o
568:         PetscScalar pt_cube,pt_five
569:         PetscScalar pt_four
570:         PetscScalar pt_val1,pt_val2,a_io2
571:         PetscInt j

573:         pt = 0.1
574:         atom_c_init =6.7408177364816552D-022
575:         atom_h_init = 2.0
576:         atom_o_init = 1.0
577:         atom_n_init = 3.76
578:         nfuel = 1
579:         an_c(1) = 1
580:         an_h(1) = 4
581:         an_c_additive = 2
582:         an_h_additive = 6
583:         an_o_additive = 1
584:         h_init = 128799.7267952987
585:         p_init = 0.1
586:         temp = 1500

588:        k_eq( 1) =     1.75149D-05
589:        k_eq( 2) =     4.01405D-06
590:        k_eq( 3) =     6.04663D-14
591:        k_eq( 4) =     2.73612D-01
592:        k_eq( 5) =     3.25592D-03
593:        k_eq( 6) =     5.33568D+05
594:        k_eq( 7) =     2.07479D+05
595:        k_eq( 8) =     1.11841D-02
596:        k_eq( 9) =     1.72684D-03
597:        k_eq(10) =     1.98588D-07
598:        k_eq(11) =     7.23600D+27
599:        k_eq(12) =     5.73926D+49
600:        k_eq(13) =     1.00000D+00
601:        k_eq(14) =     1.64493D+16
602:        k_eq(15) =     2.73837D-29
603:        k_eq(16) =     3.27419D+50
604:        k_eq(17) =     1.72447D-23
605:        k_eq(18) =     4.24657D-06
606:        k_eq(19) =     1.16065D-14
607:        k_eq(20) =     3.28020D+25
608:        k_eq(21) =     1.06291D+00
609:        k_eq(22) =     9.11507D+02
610:        k_eq(23) =     6.02837D+03

612:        H_molar( 1) =     3.26044D+03
613:        H_molar( 2) =    -8.00407D+04
614:        H_molar( 3) =     4.05873D+04
615:        H_molar( 4) =    -3.31849D+05
616:        H_molar( 5) =    -1.93654D+05
617:        H_molar( 6) =     3.84035D+04
618:        H_molar( 7) =     4.97589D+05
619:        H_molar( 8) =     2.74483D+05
620:        H_molar( 9) =     1.30022D+05
621:        H_molar(10) =     7.58429D+04
622:        H_molar(11) =     2.42948D+05
623:        H_molar(12) =     1.44588D+05
624:        H_molar(13) =    -7.16891D+04
625:        H_molar(14) =     3.63075D+04
626:        H_molar(15) =     9.23880D+04
627:        H_molar(16) =     6.50477D+04
628:        H_molar(17) =     3.04310D+05
629:        H_molar(18) =     7.41707D+05
630:        H_molar(19) =     6.32767D+05
631:        H_molar(20) =     8.90624D+05
632:        H_molar(21) =     2.49805D+04
633:        H_molar(22) =     6.43473D+05
634:        H_molar(23) =     1.02861D+06
635:        H_molar(24) =    -6.07503D+03
636:        H_molar(25) =     1.27020D+05
637:        H_molar(26) =    -1.07011D+05

639: !
640: !=======
641:       do jb = 1,26
642:          do ib = 1,26
643:             d_eq(ib,jb) = 0.0d0
644:          end do
645:       end do

647:         an_t = 0.0
648:       do id = 1,26
649:          an_t = an_t + an_r(id)
650:       enddo

652: !====
653: !====
654:         d_eq(1,1) = -an_h(1)
655:         d_eq(1,2) = -an_h_additive
656:         d_eq(1,5) = -2
657:         d_eq(1,10) = -1
658:         d_eq(1,11) = -1
659:         d_eq(1,14) = -2
660:         d_eq(1,16) = -1
661:         d_eq(1,17) = -2
662:         d_eq(1,19) = -1
663:         d_eq(1,20) = -1
664:         d_eq(1,22) = -3
665:         d_eq(1,26) = -1

667:         d_eq(2,2) = -1*an_o_additive
668:         d_eq(2,3) = -2
669:         d_eq(2,4) = -2
670:         d_eq(2,5) = -1
671:         d_eq(2,8) = -1
672:         d_eq(2,9) = -1
673:         d_eq(2,10) = -1
674:         d_eq(2,12) = -1
675:         d_eq(2,13) = -1
676:         d_eq(2,15) = -2
677:         d_eq(2,16) = -2
678:         d_eq(2,20) = -1
679:         d_eq(2,22) = -1
680:         d_eq(2,23) = -1
681:         d_eq(2,24) = -2
682:         d_eq(2,25) = -1
683:         d_eq(2,26) = -1

685:         d_eq(6,6) = -2
686:         d_eq(6,7) = -1
687:         d_eq(6,9) = -1
688:         d_eq(6,12) = -2
689:         d_eq(6,15) = -1
690:         d_eq(6,23) = -1

692:         d_eq(4,1) = -an_c(1)
693:         d_eq(4,2) = -an_c_additive
694:         d_eq(4,4) = -1
695:         d_eq(4,13) = -1
696:         d_eq(4,17) = -2
697:         d_eq(4,18) = -1
698:         d_eq(4,19) = -1
699:         d_eq(4,20) = -1

701: !----------
702:         const2 = an_t*an_t
703:         const3 = (an_t)*sqrt(an_t)
704:         const5 = an_t*const3

706:            const_cube =  an_t*an_t*an_t
707:            const_four =  const2*const2
708:            const_five =  const2*const_cube
709:            const_six  = const_cube*const_cube
710:            pt_cube = pt*pt*pt
711:            pt_four = pt_cube*pt
712:            pt_five = pt_cube*pt*pt

714:            i_cc = an_c(1)
715:            i_hh = an_h(1)
716:            ai_o2     = i_cc + float(i_hh)/4.0
717:            i_co2_h2o = i_cc + i_hh/2
718:            i_h2o     = i_hh/2
719:            ai_pwr1  = 1 + i_cc + float(i_hh)/4.0
720:            i_pwr2  = i_cc + i_hh/2
721:            idiff     = (i_cc + i_h2o) - (ai_o2 + 1)

723:            pt1   = pt**(ai_pwr1)
724:            an_tot1 = an_t**(ai_pwr1)
725:            pt_val1 = (pt/an_t)**(ai_pwr1)
726:            an_tot1_d = an_tot1*an_t

728:            pt2   = pt**(i_pwr2)
729:            an_tot2 = an_t**(i_pwr2)
730:            pt_val2 = (pt/an_t)**(i_pwr2)
731:            an_tot2_d = an_tot2*an_t

733:            d_eq(5,1) =                                                  &
734:      &           -(an_r(4)**i_cc)*(an_r(5)**i_h2o)                      &
735:      &           *((pt/an_t)**idiff) *(-idiff/an_t)

737:            do jj = 2,26
738:               d_eq(5,jj) = d_eq(5,1)
739:            enddo

741:            d_eq(5,1) = d_eq(5,1) + k_eq(11)* (an_r(3) **ai_o2)

743:            d_eq(5,3) = d_eq(5,3) + k_eq(11)* (ai_o2*an_r(3)**(ai_o2-1)) &
744:      &                           * an_r(1)

746:            d_eq(5,4) = d_eq(5,4) - (i_cc*an_r(4)**(i_cc-1))*            &
747:      &                           (an_r(5)**i_h2o)* ((pt/an_t)**idiff)
748:            d_eq(5,5) = d_eq(5,5)                                        &
749:      &               - (i_h2o*(an_r(5)**(i_h2o-1)))                     &
750:      &               * (an_r(4)**i_cc)* ((pt/an_t)**idiff)

752:            d_eq(3,1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
753:            do jj = 2,26
754:               d_eq(3,jj) = d_eq(3,1)
755:            enddo

757:            d_eq(3,2) = d_eq(3,2) + k_eq(12)* (an_r(3)**3)

759:            d_eq(3,3) = d_eq(3,3) + k_eq(12)* (3*an_r(3)**2)*an_r(2)

761:            d_eq(3,4) = d_eq(3,4) - 2*an_r(4)*(an_r(5)**3)*(pt/an_t)

763:            d_eq(3,5) = d_eq(3,5) - 3*(an_r(5)**2)*(an_r(4)**2)*(pt/an_t)
764: !     &                           *(pt_five/const_five)

766:            do ii3 = 1,26
767:               d_eq(3,ii3) = 0.0d0
768:            enddo
769:            d_eq(3,2) = 1.0d0

771:         d_eq(7,1) = pt*an_r(11)*(-1.0)/const2                           &
772:      &            -k_eq(1)*sqrt(pt)*sqrt(an_r(14)+1d-50)*(-0.5/const3)

774:         do jj = 2,26
775:            d_eq(7,jj) = d_eq(7,1)
776:         enddo

778:         d_eq(7,11) = d_eq(7,11) + pt/an_t
779:         d_eq(7,14) = d_eq(7,14)                                         &
780:      &            - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14)+1d-50)*an_t)))

782:         d_eq(8,1) = pt*an_r(8)*(-1.0)/const2                            &
783:      &            -k_eq(2)*sqrt(pt)*sqrt(an_r(3)+1.0d-50)*(-0.5/const3)

785:         do jj = 2,26
786:            d_eq(8,jj) = d_eq(8,1)
787:         enddo

789:         d_eq(8,3) = d_eq(8,3)                                           &
790:      &            -k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3)+1.0d-50)*an_t)))
791:         d_eq(8,8) = d_eq(8,8) + pt/an_t

793:         d_eq(9,1) = pt*an_r(7)*(-1.0)/const2                            &
794:      &            -k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3)

796:         do jj = 2,26
797:            d_eq(9,jj) = d_eq(9,1)
798:         enddo

800:         d_eq(9,7) = d_eq(9,7) + pt/an_t
801:         d_eq(9,6) = d_eq(9,6)                                           &
802:      &             -k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))

804:         d_eq(10,1) = pt*an_r(10)*(-1.0)/const2                          &
805:      &             -k_eq(4)*(pt)*sqrt((an_r(3)+1.0d-50)                 &
806:      &       *an_r(14))*(-1.0/const2)
807:         do jj = 2,26
808:            d_eq(10,jj) = d_eq(10,1)
809:         enddo

811:         d_eq(10,3) = d_eq(10,3)                                         &
812:      &           -k_eq(4)*(pt)*sqrt(an_r(14))                           &
813:      &           *(0.5/(sqrt(an_r(3)+1.0d-50)*an_t))
814:         d_eq(10,10) = d_eq(10,10) + pt/an_t
815:         d_eq(10,14) = d_eq(10,14)                                       &
816:      &           -k_eq(4)*(pt)*sqrt(an_r(3)+1.0d-50)                    &
817:      &            *(0.5/(sqrt(an_r(14)+1.0d-50)*an_t))

819:         d_eq(11,1) = pt*an_r(9)*(-1.0)/const2                           &
820:      &             -k_eq(5)*(pt)*sqrt((an_r(3)+1.0d-50)*an_r(6))        &
821:      &             *(-1.0/const2)

823:         do jj = 2,26
824:            d_eq(11,jj) = d_eq(11,1)
825:         enddo

827:         d_eq(11,3) = d_eq(11,3)                                         &
828:      &            -k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/                     &
829:      &       (sqrt(an_r(3)+1.0d-50)*an_t))
830:         d_eq(11,6) = d_eq(11,6)                                         &
831:      &            -k_eq(5)*(pt)*sqrt(an_r(3)+1.0d-50)                   &
832:      &       *(0.5/(sqrt(an_r(6))*an_t))
833:         d_eq(11,9) = d_eq(11,9) + pt/an_t

835:         d_eq(12,1) = pt*an_r(5)*(-1.0)/const2                           &
836:      &             -k_eq(6)*(pt**1.5)*sqrt(an_r(3)+1.0d-50)             &
837:      &             *(an_r(14))*(-1.5/const5)

839:         do jj = 2,26
840:            d_eq(12,jj) = d_eq(12,1)
841:         enddo

843:         d_eq(12,3) = d_eq(12,3)                                         &
844:      &            -k_eq(6)*(pt**1.5)*((an_r(14)+1.0d-50)/const3)        &
845:      &            *(0.5/sqrt(an_r(3)+1.0d-50))

847:         d_eq(12,5) = d_eq(12,5) + pt/an_t
848:         d_eq(12,14) = d_eq(12,14)                                       &
849:      &            -k_eq(6)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)

851:         d_eq(13,1) = pt*an_r(4)*(-1.0)/const2                           &
852:      &             -k_eq(7)*(pt**1.5)*sqrt(an_r(3)+1.0d-50)             &
853:      &             *(an_r(13))*(-1.5/const5)

855:         do jj = 2,26
856:            d_eq(13,jj) = d_eq(13,1)
857:         enddo

859:         d_eq(13,3) = d_eq(13,3)                                         &
860:      &            -k_eq(7)*(pt**1.5)*(an_r(13)/const3)                  &
861:      &            *(0.5/sqrt(an_r(3)+1.0d-50))

863:         d_eq(13,4) = d_eq(13,4) + pt/an_t
864:         d_eq(13,13) = d_eq(13,13)                                       &
865:      &            -k_eq(7)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)

867:         d_eq(14,1) = pt*an_r(15)*(-1.0)/const2                          &
868:      &             -k_eq(8)*(pt**1.5)*sqrt(an_r(3)+1.0d-50)             &
869:      &             *(an_r(9))*(-1.5/const5)

871:         do jj = 2,26
872:            d_eq(14,jj) = d_eq(14,1)
873:         enddo

875:         d_eq(14,3) = d_eq(14,3)                                         &
876:      &            -k_eq(8)*(pt**1.5)*(an_r(9)/const3)                   &
877:      &            *(0.5/sqrt(an_r(3)+1.0d-50))
878:         d_eq(14,9) = d_eq(14,9)                                         &
879:      &            -k_eq(8)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
880:         d_eq(14,15) = d_eq(14,15)+ pt/an_t

882:         d_eq(15,1) = pt*an_r(16)*(-1.0)/const2                          &
883:      &             -k_eq(9)*(pt**1.5)*sqrt(an_r(14)+1.0d-50)            &
884:      &             *(an_r(3))*(-1.5/const5)

886:         do jj = 2,26
887:            d_eq(15,jj) = d_eq(15,1)
888:         enddo

890:         d_eq(15,3) = d_eq(15,3)                                         &
891:      &            -k_eq(9)*(pt**1.5)*(sqrt(an_r(14)+1.0d-50)/const3)
892:         d_eq(15,14) = d_eq(15,14)                                       &
893:      &            -k_eq(9)*(pt**1.5)*(an_r(3)/const3)                   &
894:      &            *(0.5/sqrt(an_r(14)+1.0d-50))
895:         d_eq(15,16) = d_eq(15,16) + pt/an_t

897:         d_eq(16,1) = pt*an_r(12)*(-1.0)/const2                          &
898:      &             -k_eq(10)*(pt**1.5)*sqrt(an_r(3)+1.0d-50)            &
899:      &             *(an_r(6))*(-1.5/const5)

901:         do jj = 2,26
902:            d_eq(16,jj) = d_eq(16,1)
903:         enddo

905:         d_eq(16,3) = d_eq(16,3)                                         &
906:      &             -k_eq(10)*(pt**1.5)*(an_r(6)/const3)                 &
907:      &             *(0.5/sqrt(an_r(3)+1.0d-50))

909:         d_eq(16,6) = d_eq(16,6)                                         &
910:      &             -k_eq(10)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
911:         d_eq(16,12) = d_eq(16,12) + pt/an_t

913:         const_cube =  an_t*an_t*an_t
914:         const_four =  const2*const2

916:         d_eq(17,1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
917:      &             - k_eq(15) * an_r(17)*pt * (-1/const2)
918:         do jj = 2,26
919:            d_eq(17,jj) = d_eq(17,1)
920:         enddo
921:         d_eq(17,14) = d_eq(17,14) + an_r(18)*an_r(18)*(pt**3)/const_cube
922:         d_eq(17,17) = d_eq(17,17) - k_eq(15)*pt/an_t
923:         d_eq(17,18) = d_eq(17,18) + 2*an_r(18)*an_r(14)                 &
924:      &                            *(pt**3)/const_cube

926:         d_eq(18,1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube)          &
927:      &             - k_eq(16) * an_r(3)*an_r(18)*an_r(18)               &
928:      &              * (pt*pt*pt) * (-3/const_four)
929:         do jj = 2,26
930:            d_eq(18,jj) = d_eq(18,1)
931:         enddo
932:         d_eq(18,3) = d_eq(18,3)                                         &
933:      &             - k_eq(16) *an_r(18)* an_r(18)*pt*pt*pt /const_cube
934:         d_eq(18,13) = d_eq(18,13)                                       &
935:      &              + 2* an_r(13)*pt*pt /const2
936:         d_eq(18,18) = d_eq(18,18) -k_eq(16)*an_r(3)                     &
937:      &              * 2*an_r(18)*pt*pt*pt/const_cube

939: !====for eq 19

941:         d_eq(19,1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube)           &
942:      &             - k_eq(17)*an_r(13)*an_r(10)*pt*pt * (-2/const_cube)
943:         do jj = 2,26
944:            d_eq(19,jj) = d_eq(19,1)
945:         enddo
946:         d_eq(19,13) = d_eq(19,13)                                       &
947:      &             - k_eq(17) *an_r(10)*pt*pt /const2
948:         d_eq(19,10) = d_eq(19,10)                                       &
949:      &             - k_eq(17) *an_r(13)*pt*pt /const2
950:         d_eq(19,3) = d_eq(19,3) + an_r(19)*pt*pt/const2
951:         d_eq(19,19) = d_eq(19,19) + an_r(3)*pt*pt/const2
952: !====for eq 20

954:         d_eq(20,1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube)          &
955:      &             - k_eq(18) * an_r(19)*an_r(8)*pt*pt * (-2/const_cube)
956:         do jj = 2,26
957:            d_eq(20,jj) = d_eq(20,1)
958:         enddo
959:         d_eq(20,8) = d_eq(20,8)                                         &
960:      &             - k_eq(18) *an_r(19)*pt*pt /const2
961:         d_eq(20,19) = d_eq(20,19)                                       &
962:      &             - k_eq(18) *an_r(8)*pt*pt /const2
963:         d_eq(20,20) = d_eq(20,20) + an_r(21)*pt*pt/const2
964:         d_eq(20,21) = d_eq(20,21) + an_r(20)*pt*pt/const2

966: !========
967: !====for eq 21

969:         d_eq(21,1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube)          &
970:      &             - k_eq(19)*an_r(7)*an_r(8)*pt*pt * (-2/const_cube)
971:         do jj = 2,26
972:            d_eq(21,jj) = d_eq(21,1)
973:         enddo
974:         d_eq(21,7) = d_eq(21,7)                                         &
975:      &             - k_eq(19) *an_r(8)*pt*pt /const2
976:         d_eq(21,8) = d_eq(21,8)                                         &
977:      &             - k_eq(19) *an_r(7)*pt*pt /const2
978:         d_eq(21,21) = d_eq(21,21) + an_r(23)*pt*pt/const2
979:         d_eq(21,23) = d_eq(21,23) + an_r(21)*pt*pt/const2

981: !========
982: !  for 22
983:         d_eq(22,1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube)           &
984:      &         -k_eq(20)*an_r(21)*an_r(22)*pt*pt * (-2/const_cube)
985:         do jj = 2,26
986:            d_eq(22,jj) = d_eq(22,1)
987:         enddo
988:         d_eq(22,21) = d_eq(22,21)                                       &
989:      &             - k_eq(20) *an_r(22)*pt*pt /const2
990:         d_eq(22,22) = d_eq(22,22)                                       &
991:      &             - k_eq(20) *an_r(21)*pt*pt /const2
992:         d_eq(22,11) = d_eq(22,11) + an_r(5)*pt*pt/(const2)
993:         d_eq(22,5) = d_eq(22,5) + an_r(11)*pt*pt/(const2)

995: !========
996: !  for 23

998:         d_eq(23,1) = an_r(24)*(pt)*(-1/const2)                          &
999:      &             - k_eq(21)*an_r(21)*an_r(3)*pt*pt * (-2/const_cube)
1000:         do jj = 2,26
1001:            d_eq(23,jj) = d_eq(23,1)
1002:         enddo
1003:         d_eq(23,3) = d_eq(23,3)                                         &
1004:      &             - k_eq(21) *an_r(21)*pt*pt /const2
1005:         d_eq(23,21) = d_eq(23,21)                                       &
1006:      &             - k_eq(21) *an_r(3)*pt*pt /const2
1007:         d_eq(23,24) = d_eq(23,24) + pt/(an_t)

1009: !========
1010: !  for 24
1011:         d_eq(24,1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube)           &
1012:      &             - k_eq(22)*an_r(24)*an_r(8)*pt*pt * (-2/const_cube)
1013:         do jj = 2,26
1014:            d_eq(24,jj) = d_eq(24,1)
1015:         enddo
1016:         d_eq(24,8) = d_eq(24,8)                                         &
1017:      &             - k_eq(22) *an_r(24)*pt*pt /const2
1018:         d_eq(24,24) = d_eq(24,24)                                       &
1019:      &             - k_eq(22) *an_r(8)*pt*pt /const2
1020:         d_eq(24,3) = d_eq(24,3) + an_r(25)*pt*pt/const2
1021:         d_eq(24,25) = d_eq(24,25) + an_r(3)*pt*pt/const2

1023: !========
1024: !for 25

1026:         d_eq(25,1) = an_r(26)*(pt)*(-1/const2)                          &
1027:      &       - k_eq(23)*an_r(21)*an_r(10)*pt*pt * (-2/const_cube)
1028:         do jj = 2,26
1029:            d_eq(25,jj) = d_eq(25,1)
1030:         enddo
1031:         d_eq(25,10) = d_eq(25,10)                                       &
1032:      &             - k_eq(23) *an_r(21)*pt*pt /const2
1033:         d_eq(25,21) = d_eq(25,21)                                       &
1034:      &             - k_eq(23) *an_r(10)*pt*pt /const2
1035:         d_eq(25,26) = d_eq(25,26) + pt/(an_t)

1037: !============
1038: !  for 26
1039:         d_eq(26,20) = -1
1040:         d_eq(26,22) = -1
1041:         d_eq(26,23) = -1
1042:         d_eq(26,21) = 1
1043:         d_eq(26,24) = 1
1044:         d_eq(26,25) = 1
1045:         d_eq(26,26) = 1

1047:            do j = 1,26
1048:          do i = 1,26
1049:                 write(44,*)i,j,d_eq(i,j)
1050:               enddo
1051:            enddo

1053:         return
1054:         end

1056:       subroutine ShashiPostCheck(ls,X,Y,W,c_Y,c_W,dummy)
1057: #include <petsc/finclude/petscsnes.h>
1058:       use petscsnes
1059:       implicit none
1060:       SNESLineSearch ls
1061:       PetscErrorCode ierr
1062:       Vec X,Y,W
1063:       PetscObject dummy
1064:       PetscBool c_Y,c_W
1065:       PetscScalar,pointer :: xx(:)
1066:       PetscInt i
1067:       call VecGetArrayF90(W,xx,ierr)
1068:       do i=1,26
1069:          if (xx(i) < 0.0) then
1070:             xx(i) = 0.0
1071:             c_W = PETSC_TRUE
1072:          endif
1073:         if (xx(i) > 3.0) then
1074:            xx(i) = 3.0
1075:         endif
1076:       enddo
1077:       call VecRestoreArrayF90(W,xx,ierr)
1078:       return
1079:       end