Actual source code: shashi.F90

petsc-3.13.4 2020-08-01
Report Typos and Errors


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

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

 40:       external FormFunction, FormJacobian
 41:       external ShashiPostCheck

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

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

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

 75:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

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

 81: !  Create vectors for solution and nonlinear function

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


 87: !  Create Jacobian matrix data structure

 89:       call MatCreateDense(PETSC_COMM_SELF,26,26,26,26,                          &
 90:      &                    PETSC_NULL_SCALAR,J,ierr)

 92: !  Set function evaluation routine and vector

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

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

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

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

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

114:       call SNESSetFromOptions(snes,ierr)

116: !     set initial guess

118:       call VecGetArrayF90(x,xx,ierr)
119:       call ShashiInitialGuess(xx)
120:       call VecRestoreArrayF90(x,xx,ierr)

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

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

158: !  Declarations for use with local arrays

160:       PetscScalar  lx_v(2),lf_v(2)
161:       PetscOffset  lx_i,lf_i

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

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

177:       return
178:       end

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

208:       PetscScalar lx_v(1),lf_v(1)
209:       PetscOffset lx_i,lf_i

211: !  Get pointer to vector data

213:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
214:       call MatDenseGetArray(B,lf_v,lf_i,ierr)
215:       call ShashiFormJacobian(lx_a(1),lf_a(1))
216:       call MatDenseRestoreArray(B,lf_v,lf_i,ierr)
217:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
218:       
219: !  Assemble matrix

221:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
222:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
223:       
224:       return
225:       end

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

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

239:             subroutine ShashiInitialGuess(an_r)
240: !        implicit PetscScalar (a-h,o-z)
241:         implicit none
242:         PetscScalar an_c_additive
243:         PetscScalar       an_h_additive
244:         PetscScalar an_o_additive
245:         PetscScalar   atom_c_init
246:         PetscScalar atom_h_init
247:         PetscScalar atom_n_init
248:         PetscScalar atom_o_init
249:         PetscScalar h_init
250:         PetscScalar p_init
251:         PetscInt nfuel
252:         PetscScalar temp,pt
253:         PetscScalar an_r(26),k_eq(23),f_eq(26)
254:         PetscScalar d_eq(26,26),H_molar(26)
255:         PetscInt an_h(1),an_c(1)
256:         PetscScalar part_p(26)
257:         PetscInt i_cc,i_hh,i_h2o
258:         PetscInt i_pwr2,i_co2_h2o 
259:         
260:         pt = 0.1
261:         atom_c_init =6.7408177364816552D-022 
262:         atom_h_init = 2.0
263:         atom_o_init = 1.0
264:         atom_n_init = 3.76
265:         nfuel = 1
266:         an_c(1) = 1
267:         an_h(1) = 4
268:         an_c_additive = 2
269:         an_h_additive = 6
270:         an_o_additive = 1
271:         h_init = 128799.7267952987
272:         p_init = 0.1
273:         temp = 1500


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

303:       an_r = 0
304:       an_r( 3) =     .5
305:       an_r( 6) =     1.88000
306:       an_r(14) =     1.

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



343:       subroutine ShashiFormFunction(an_r,f_eq)
344: !       implicit PetscScalar (a-h,o-z)
345:         implicit none
346:         PetscScalar an_c_additive
347:         PetscScalar       an_h_additive
348:         PetscScalar an_o_additive
349:         PetscScalar   atom_c_init
350:         PetscScalar atom_h_init
351:         PetscScalar atom_n_init
352:         PetscScalar atom_o_init
353:         PetscScalar h_init
354:         PetscScalar p_init
355:         PetscInt nfuel
356:         PetscScalar temp,pt
357:        PetscScalar an_r(26),k_eq(23),f_eq(26)
358:        PetscScalar d_eq(26,26),H_molar(26)
359:        PetscInt an_h(1),an_c(1)
360:        PetscScalar part_p(26),idiff
361:         PetscInt i_cc,i_hh,i_h2o
362:         PetscInt i_pwr2,i_co2_h2o 
363:         PetscScalar an_t,sum_h,pt_cubed,pt_five
364:         PetscScalar pt_four,pt_val1,pt_val2
365:         PetscScalar a_io2
366:         PetscInt i,ip
367:         pt = 0.1
368:         atom_c_init =6.7408177364816552D-022 
369:         atom_h_init = 2.0
370:         atom_o_init = 1.0
371:         atom_n_init = 3.76
372:         nfuel = 1
373:         an_c(1) = 1
374:         an_h(1) = 4
375:         an_c_additive = 2
376:         an_h_additive = 6
377:         an_o_additive = 1
378:         h_init = 128799.7267952987
379:         p_init = 0.1
380:         temp = 1500
381:         
382:        k_eq( 1) =     1.75149D-05
383:        k_eq( 2) =     4.01405D-06
384:        k_eq( 3) =     6.04663D-14
385:        k_eq( 4) =     2.73612D-01
386:        k_eq( 5) =     3.25592D-03
387:        k_eq( 6) =     5.33568D+05
388:        k_eq( 7) =     2.07479D+05
389:        k_eq( 8) =     1.11841D-02
390:        k_eq( 9) =     1.72684D-03
391:        k_eq(10) =     1.98588D-07
392:        k_eq(11) =     7.23600D+27
393:        k_eq(12) =     5.73926D+49
394:        k_eq(13) =     1.00000D+00
395:        k_eq(14) =     1.64493D+16
396:        k_eq(15) =     2.73837D-29
397:        k_eq(16) =     3.27419D+50
398:        k_eq(17) =     1.72447D-23
399:        k_eq(18) =     4.24657D-06
400:        k_eq(19) =     1.16065D-14
401:        k_eq(20) =     3.28020D+25
402:        k_eq(21) =     1.06291D+00
403:        k_eq(22) =     9.11507D+02
404:        k_eq(23) =     6.02837D+03

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

436:       do i = 1,26
437:          an_t = an_t + an_r(i)
438:       enddo

440:         f_eq(1) = atom_h_init                                           &
441:      &          - (an_h(1)*an_r(1) + an_h_additive*an_r(2)              &
442:      &              + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14)      &
443:      &              + an_r(16)+ 2*an_r(17) + an_r(19)                   &
444:      &              +an_r(20) + 3*an_r(22)+an_r(26)) 


447:         f_eq(2) = atom_o_init                                           &
448:      &          - (an_o_additive*an_r(2) + 2*an_r(3)                    &
449:      &             + 2*an_r(4) + an_r(5)                                &
450:      &             + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
451:      &             + 2*an_r(15) + 2*an_r(16)+ an_r(20) + an_r(22)       &
452:      &             + an_r(23) + 2*an_r(24) + 1*an_r(25)+an_r(26))


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

457:         f_eq(4) = atom_c_init                                           &
458:      &          - (an_c(1)*an_r(1) + an_c_additive * an_r(2)            &
459:      &          + an_r(4) + an_r(13)+ 2*an_r(17) + an_r(18)             &
460:      &          + an_r(19) + an_r(20))



464:         
465:         do ip = 1,26
466:            part_p(ip) = (an_r(ip)/an_t) * pt
467:         enddo

469:         i_cc      = an_c(1)
470:         i_hh      = an_h(1)
471:         a_io2      = i_cc + i_hh/4.0
472:         i_h2o     = i_hh/2
473:         idiff     = (i_cc + i_h2o) - (a_io2 + 1)

475:         f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2                       &
476:      &          - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
477: !           write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
478: !          stop
479:         f_eq(6) = atom_n_init                                           &
480:      &          - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12)           &
481:      &              + an_r(15)                                          &
482:      &              + an_r(23))




487:       f_eq( 7) = part_p(11)                                             &
488:      &         - (k_eq( 1) * sqrt(part_p(14)+1d-23))
489:       f_eq( 8) = part_p( 8)                                             &
490:      &         - (k_eq( 2) * sqrt(part_p( 3)+1d-23))

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

495:       f_eq(10) = part_p(10)                                             &
496:      &         - (k_eq( 4) * sqrt(part_p( 3)+1d-23))                    &
497:      &         * sqrt(part_p(14))

499:       f_eq(11) = part_p( 9)                                             &
500:      &         - (k_eq( 5) * sqrt(part_p( 3)+1d-23))                    &
501:      &         * sqrt(part_p( 6)+1d-23)
502:       f_eq(12) = part_p( 5)                                             &
503:      &         - (k_eq( 6) * sqrt(part_p( 3)+1d-23))                    &
504:      &         * (part_p(14))


507:       f_eq(13) = part_p( 4)                                             &
508:      &         - (k_eq( 7) * sqrt(part_p(3)+1.0d-23))                   &
509:      &         * (part_p(13))

511:       f_eq(14) = part_p(15)                                             &
512:      &         - (k_eq( 8) * sqrt(part_p(3)+1.0d-50))                   &
513:      &         * (part_p( 9))
514:       f_eq(15) = part_p(16)                                             &
515:      &         - (k_eq( 9) * part_p( 3))                                &
516:      &         * sqrt(part_p(14)+1d-23)

518:       f_eq(16) = part_p(12)                                             &
519:      &         - (k_eq(10) * sqrt(part_p( 3)+1d-23))                    &
520:      &         * (part_p( 6))

522:       f_eq(17) = part_p(14)*part_p(18)**2                               &
523:      &         - (k_eq(15) * part_p(17))

525:       f_eq(18) = (part_p(13)**2)                                        &
526:      &     - (k_eq(16) * part_p(3)*part_p(18)**2)
527:       print*,f_eq(18),part_p(3),part_p(18),part_p(13),k_eq(16)

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

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

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


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


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

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

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

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

548:              do i = 1,26
549:                  write(44,*)i,f_eq(i)
550:               enddo
551:               
552:       return
553:       end

555:       subroutine ShashiFormJacobian(an_r,d_eq)
556: !        implicit PetscScalar (a-h,o-z)
557:         implicit none
558:         PetscScalar an_c_additive
559:         PetscScalar       an_h_additive
560:         PetscScalar an_o_additive
561:         PetscScalar   atom_c_init
562:         PetscScalar atom_h_init
563:         PetscScalar atom_n_init
564:         PetscScalar atom_o_init
565:         PetscScalar h_init
566:         PetscScalar p_init
567:         PetscInt nfuel
568:         PetscScalar temp,pt
569:         PetscScalar an_t,ai_o2,sum_h
570:         PetscScalar an_tot1_d,an_tot1
571:         PetscScalar an_tot2_d,an_tot2
572:         PetscScalar const5,const3,const2
573:         PetscScalar   const_cube
574:         PetscScalar   const_five
575:         PetscScalar   const_four
576:         PetscScalar   const_six
577:         PetscInt jj,jb,ii3,id,ib,ip,i
578:         PetscScalar   pt2,pt1
579:         PetscScalar an_r(26),k_eq(23),f_eq(26)
580:         PetscScalar d_eq(26,26),H_molar(26)
581:         PetscInt an_h(1),an_c(1)
582:         PetscScalar ai_pwr1,part_p(26),idiff
583:         PetscInt i_cc,i_hh
584:         PetscInt i_h2o,i_pwr2,i_co2_h2o
585:         PetscScalar pt_cube,pt_five
586:         PetscScalar pt_four
587:         PetscScalar pt_val1,pt_val2,a_io2
588:         PetscInt j

590:         pt = 0.1
591:         atom_c_init =6.7408177364816552D-022 
592:         atom_h_init = 2.0
593:         atom_o_init = 1.0
594:         atom_n_init = 3.76
595:         nfuel = 1
596:         an_c(1) = 1
597:         an_h(1) = 4
598:         an_c_additive = 2
599:         an_h_additive = 6
600:         an_o_additive = 1
601:         h_init = 128799.7267952987
602:         p_init = 0.1
603:         temp = 1500

605:        k_eq( 1) =     1.75149D-05
606:        k_eq( 2) =     4.01405D-06
607:        k_eq( 3) =     6.04663D-14
608:        k_eq( 4) =     2.73612D-01
609:        k_eq( 5) =     3.25592D-03
610:        k_eq( 6) =     5.33568D+05
611:        k_eq( 7) =     2.07479D+05
612:        k_eq( 8) =     1.11841D-02
613:        k_eq( 9) =     1.72684D-03
614:        k_eq(10) =     1.98588D-07
615:        k_eq(11) =     7.23600D+27
616:        k_eq(12) =     5.73926D+49
617:        k_eq(13) =     1.00000D+00
618:        k_eq(14) =     1.64493D+16
619:        k_eq(15) =     2.73837D-29
620:        k_eq(16) =     3.27419D+50
621:        k_eq(17) =     1.72447D-23
622:        k_eq(18) =     4.24657D-06
623:        k_eq(19) =     1.16065D-14
624:        k_eq(20) =     3.28020D+25
625:        k_eq(21) =     1.06291D+00
626:        k_eq(22) =     9.11507D+02
627:        k_eq(23) =     6.02837D+03

629:        H_molar( 1) =     3.26044D+03
630:        H_molar( 2) =    -8.00407D+04
631:        H_molar( 3) =     4.05873D+04
632:        H_molar( 4) =    -3.31849D+05
633:        H_molar( 5) =    -1.93654D+05
634:        H_molar( 6) =     3.84035D+04
635:        H_molar( 7) =     4.97589D+05
636:        H_molar( 8) =     2.74483D+05
637:        H_molar( 9) =     1.30022D+05
638:        H_molar(10) =     7.58429D+04
639:        H_molar(11) =     2.42948D+05
640:        H_molar(12) =     1.44588D+05
641:        H_molar(13) =    -7.16891D+04
642:        H_molar(14) =     3.63075D+04
643:        H_molar(15) =     9.23880D+04
644:        H_molar(16) =     6.50477D+04
645:        H_molar(17) =     3.04310D+05
646:        H_molar(18) =     7.41707D+05
647:        H_molar(19) =     6.32767D+05
648:        H_molar(20) =     8.90624D+05
649:        H_molar(21) =     2.49805D+04
650:        H_molar(22) =     6.43473D+05
651:        H_molar(23) =     1.02861D+06
652:        H_molar(24) =    -6.07503D+03
653:        H_molar(25) =     1.27020D+05
654:        H_molar(26) =    -1.07011D+05
655:        
656: !
657: !=======
658:       do jb = 1,26
659:          do ib = 1,26
660:             d_eq(ib,jb) = 0.0d0
661:          end do
662:       end do

664:         an_t = 0.0
665:       do id = 1,26
666:          an_t = an_t + an_r(id)
667:       enddo

669: !====
670: !====
671:         d_eq(1,1) = -an_h(1)
672:         d_eq(1,2) = -an_h_additive
673:         d_eq(1,5) = -2    
674:         d_eq(1,10) = -1    
675:         d_eq(1,11) = -1    
676:         d_eq(1,14) = -2  
677:         d_eq(1,16) = -1  
678:         d_eq(1,17) = -2
679:         d_eq(1,19) = -1  
680:         d_eq(1,20) = -1  
681:         d_eq(1,22) = -3
682:         d_eq(1,26) = -1

684:         d_eq(2,2) = -1*an_o_additive
685:         d_eq(2,3) = -2 
686:         d_eq(2,4) = -2
687:         d_eq(2,5) = -1    
688:         d_eq(2,8) = -1    
689:         d_eq(2,9) = -1    
690:         d_eq(2,10) = -1    
691:         d_eq(2,12) = -1
692:         d_eq(2,13) = -1
693:         d_eq(2,15) = -2
694:         d_eq(2,16) = -2
695:         d_eq(2,20) = -1
696:         d_eq(2,22) = -1
697:         d_eq(2,23) = -1
698:         d_eq(2,24) = -2
699:         d_eq(2,25) = -1
700:         d_eq(2,26) = -1



704:         d_eq(6,6) = -2
705:         d_eq(6,7) = -1
706:         d_eq(6,9) = -1
707:         d_eq(6,12) = -2
708:         d_eq(6,15) = -1
709:         d_eq(6,23) = -1



713:         d_eq(4,1) = -an_c(1)
714:         d_eq(4,2) = -an_c_additive
715:         d_eq(4,4) = -1
716:         d_eq(4,13) = -1
717:         d_eq(4,17) = -2
718:         d_eq(4,18) = -1
719:         d_eq(4,19) = -1
720:         d_eq(4,20) = -1


723: !----------
724:         const2 = an_t*an_t
725:         const3 = (an_t)*sqrt(an_t)
726:         const5 = an_t*const3


729:            const_cube =  an_t*an_t*an_t
730:            const_four =  const2*const2
731:            const_five =  const2*const_cube
732:            const_six  = const_cube*const_cube
733:            pt_cube = pt*pt*pt
734:            pt_four = pt_cube*pt
735:            pt_five = pt_cube*pt*pt

737:            i_cc = an_c(1)
738:            i_hh = an_h(1)
739:            ai_o2     = i_cc + float(i_hh)/4.0
740:            i_co2_h2o = i_cc + i_hh/2
741:            i_h2o     = i_hh/2
742:            ai_pwr1  = 1 + i_cc + float(i_hh)/4.0
743:            i_pwr2  = i_cc + i_hh/2
744:            idiff     = (i_cc + i_h2o) - (ai_o2 + 1)

746:            pt1   = pt**(ai_pwr1)
747:            an_tot1 = an_t**(ai_pwr1)
748:            pt_val1 = (pt/an_t)**(ai_pwr1)
749:            an_tot1_d = an_tot1*an_t

751:            pt2   = pt**(i_pwr2)
752:            an_tot2 = an_t**(i_pwr2)
753:            pt_val2 = (pt/an_t)**(i_pwr2)
754:            an_tot2_d = an_tot2*an_t


757:            d_eq(5,1) =                                                  &
758:      &           -(an_r(4)**i_cc)*(an_r(5)**i_h2o)                      &
759:      &           *((pt/an_t)**idiff) *(-idiff/an_t)


762:            do jj = 2,26
763:               d_eq(5,jj) = d_eq(5,1)
764:            enddo

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

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

771:            d_eq(5,4) = d_eq(5,4) - (i_cc*an_r(4)**(i_cc-1))*            &
772:      &                           (an_r(5)**i_h2o)* ((pt/an_t)**idiff)
773:            d_eq(5,5) = d_eq(5,5)                                        &
774:      &               - (i_h2o*(an_r(5)**(i_h2o-1)))                     &
775:      &               * (an_r(4)**i_cc)* ((pt/an_t)**idiff)



779:            d_eq(3,1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
780:            do jj = 2,26
781:               d_eq(3,jj) = d_eq(3,1)
782:            enddo


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



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

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


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

797:            do ii3 = 1,26
798:               d_eq(3,ii3) = 0.0d0
799:            enddo
800:            d_eq(3,2) = 1.0d0



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

807:         do jj = 2,26
808:            d_eq(7,jj) = d_eq(7,1)
809:         enddo

811:         d_eq(7,11) = d_eq(7,11) + pt/an_t
812:         d_eq(7,14) = d_eq(7,14)                                         &
813:      &            - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14)+1d-50)*an_t)))


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

819:         do jj = 2,26
820:            d_eq(8,jj) = d_eq(8,1)
821:         enddo

823:         d_eq(8,3) = d_eq(8,3)                                           &
824:      &            -k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3)+1.0d-50)*an_t)))
825:         d_eq(8,8) = d_eq(8,8) + pt/an_t


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

831:         do jj = 2,26
832:            d_eq(9,jj) = d_eq(9,1)
833:         enddo

835:         d_eq(9,7) = d_eq(9,7) + pt/an_t
836:         d_eq(9,6) = d_eq(9,6)                                           &
837:      &             -k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))


840:         d_eq(10,1) = pt*an_r(10)*(-1.0)/const2                          &
841:      &             -k_eq(4)*(pt)*sqrt((an_r(3)+1.0d-50)                 &
842:      &       *an_r(14))*(-1.0/const2)
843:         do jj = 2,26
844:            d_eq(10,jj) = d_eq(10,1)
845:         enddo

847:         d_eq(10,3) = d_eq(10,3)                                         &
848:      &           -k_eq(4)*(pt)*sqrt(an_r(14))                           &
849:      &           *(0.5/(sqrt(an_r(3)+1.0d-50)*an_t))
850:         d_eq(10,10) = d_eq(10,10) + pt/an_t
851:         d_eq(10,14) = d_eq(10,14)                                       &
852:      &           -k_eq(4)*(pt)*sqrt(an_r(3)+1.0d-50)                    &
853:      &            *(0.5/(sqrt(an_r(14)+1.0d-50)*an_t))

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

859:         do jj = 2,26
860:            d_eq(11,jj) = d_eq(11,1)
861:         enddo

863:         d_eq(11,3) = d_eq(11,3)                                         &
864:      &            -k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/                     &
865:      &       (sqrt(an_r(3)+1.0d-50)*an_t))
866:         d_eq(11,6) = d_eq(11,6)                                         &
867:      &            -k_eq(5)*(pt)*sqrt(an_r(3)+1.0d-50)                   &
868:      &       *(0.5/(sqrt(an_r(6))*an_t))
869:         d_eq(11,9) = d_eq(11,9) + pt/an_t



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


878:         do jj = 2,26
879:            d_eq(12,jj) = d_eq(12,1)
880:         enddo

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

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


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

895:         do jj = 2,26
896:            d_eq(13,jj) = d_eq(13,1)
897:         enddo

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

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

907:  


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

914:         do jj = 2,26
915:            d_eq(14,jj) = d_eq(14,1)
916:         enddo

918:         d_eq(14,3) = d_eq(14,3)                                         &
919:      &            -k_eq(8)*(pt**1.5)*(an_r(9)/const3)                   &
920:      &            *(0.5/sqrt(an_r(3)+1.0d-50))
921:         d_eq(14,9) = d_eq(14,9)                                         &
922:      &            -k_eq(8)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
923:         d_eq(14,15) = d_eq(14,15)+ pt/an_t



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

931:         do jj = 2,26
932:            d_eq(15,jj) = d_eq(15,1)
933:         enddo

935:         d_eq(15,3) = d_eq(15,3)                                         &
936:      &            -k_eq(9)*(pt**1.5)*(sqrt(an_r(14)+1.0d-50)/const3)
937:         d_eq(15,14) = d_eq(15,14)                                       &
938:      &            -k_eq(9)*(pt**1.5)*(an_r(3)/const3)                   &
939:      &            *(0.5/sqrt(an_r(14)+1.0d-50))  
940:         d_eq(15,16) = d_eq(15,16) + pt/an_t  


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

947:         do jj = 2,26
948:            d_eq(16,jj) = d_eq(16,1)
949:         enddo

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

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





963:         const_cube =  an_t*an_t*an_t
964:         const_four =  const2*const2


967:         d_eq(17,1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
968:      &             - k_eq(15) * an_r(17)*pt * (-1/const2)
969:         do jj = 2,26
970:            d_eq(17,jj) = d_eq(17,1)
971:         enddo
972:         d_eq(17,14) = d_eq(17,14) + an_r(18)*an_r(18)*(pt**3)/const_cube
973:         d_eq(17,17) = d_eq(17,17) - k_eq(15)*pt/an_t
974:         d_eq(17,18) = d_eq(17,18) + 2*an_r(18)*an_r(14)                 &
975:      &                            *(pt**3)/const_cube


978:         d_eq(18,1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube)          &
979:      &             - k_eq(16) * an_r(3)*an_r(18)*an_r(18)               &
980:      &              * (pt*pt*pt) * (-3/const_four)
981:         do jj = 2,26
982:            d_eq(18,jj) = d_eq(18,1)
983:         enddo
984:         d_eq(18,3) = d_eq(18,3)                                         &
985:      &             - k_eq(16) *an_r(18)* an_r(18)*pt*pt*pt /const_cube
986:         d_eq(18,13) = d_eq(18,13)                                       &
987:      &              + 2* an_r(13)*pt*pt /const2
988:         d_eq(18,18) = d_eq(18,18) -k_eq(16)*an_r(3)                     &
989:      &              * 2*an_r(18)*pt*pt*pt/const_cube



993: !====for eq 19 

995:         d_eq(19,1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube)           &
996:      &             - k_eq(17)*an_r(13)*an_r(10)*pt*pt * (-2/const_cube)
997:         do jj = 2,26
998:            d_eq(19,jj) = d_eq(19,1)
999:         enddo
1000:         d_eq(19,13) = d_eq(19,13)                                       &
1001:      &             - k_eq(17) *an_r(10)*pt*pt /const2
1002:         d_eq(19,10) = d_eq(19,10)                                       &
1003:      &             - k_eq(17) *an_r(13)*pt*pt /const2
1004:         d_eq(19,3) = d_eq(19,3) + an_r(19)*pt*pt/const2
1005:         d_eq(19,19) = d_eq(19,19) + an_r(3)*pt*pt/const2
1006: !====for eq 20 

1008:         d_eq(20,1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube)          &
1009:      &             - k_eq(18) * an_r(19)*an_r(8)*pt*pt * (-2/const_cube)
1010:         do jj = 2,26
1011:            d_eq(20,jj) = d_eq(20,1)
1012:         enddo
1013:         d_eq(20,8) = d_eq(20,8)                                         &
1014:      &             - k_eq(18) *an_r(19)*pt*pt /const2
1015:         d_eq(20,19) = d_eq(20,19)                                       &
1016:      &             - k_eq(18) *an_r(8)*pt*pt /const2 
1017:         d_eq(20,20) = d_eq(20,20) + an_r(21)*pt*pt/const2
1018:         d_eq(20,21) = d_eq(20,21) + an_r(20)*pt*pt/const2

1020: !========
1021: !====for eq 21 

1023:         d_eq(21,1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube)          &
1024:      &             - k_eq(19)*an_r(7)*an_r(8)*pt*pt * (-2/const_cube)
1025:         do jj = 2,26
1026:            d_eq(21,jj) = d_eq(21,1)
1027:         enddo
1028:         d_eq(21,7) = d_eq(21,7)                                         &
1029:      &             - k_eq(19) *an_r(8)*pt*pt /const2
1030:         d_eq(21,8) = d_eq(21,8)                                         &
1031:      &             - k_eq(19) *an_r(7)*pt*pt /const2
1032:         d_eq(21,21) = d_eq(21,21) + an_r(23)*pt*pt/const2
1033:         d_eq(21,23) = d_eq(21,23) + an_r(21)*pt*pt/const2

1035: !======== 
1036: !  for 22  
1037:         d_eq(22,1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube)           &
1038:      &         -k_eq(20)*an_r(21)*an_r(22)*pt*pt * (-2/const_cube)
1039:         do jj = 2,26
1040:            d_eq(22,jj) = d_eq(22,1)
1041:         enddo
1042:         d_eq(22,21) = d_eq(22,21)                                       &
1043:      &             - k_eq(20) *an_r(22)*pt*pt /const2
1044:         d_eq(22,22) = d_eq(22,22)                                       &
1045:      &             - k_eq(20) *an_r(21)*pt*pt /const2
1046:         d_eq(22,11) = d_eq(22,11) + an_r(5)*pt*pt/(const2)
1047:         d_eq(22,5) = d_eq(22,5) + an_r(11)*pt*pt/(const2)



1051: !======== 
1052: !  for 23 

1054:         d_eq(23,1) = an_r(24)*(pt)*(-1/const2)                          &
1055:      &             - k_eq(21)*an_r(21)*an_r(3)*pt*pt * (-2/const_cube)
1056:         do jj = 2,26
1057:            d_eq(23,jj) = d_eq(23,1)
1058:         enddo
1059:         d_eq(23,3) = d_eq(23,3)                                         &
1060:      &             - k_eq(21) *an_r(21)*pt*pt /const2
1061:         d_eq(23,21) = d_eq(23,21)                                       &
1062:      &             - k_eq(21) *an_r(3)*pt*pt /const2
1063:         d_eq(23,24) = d_eq(23,24) + pt/(an_t)

1065: !======== 
1066: !  for 24 
1067:         d_eq(24,1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube)           &
1068:      &             - k_eq(22)*an_r(24)*an_r(8)*pt*pt * (-2/const_cube)
1069:         do jj = 2,26
1070:            d_eq(24,jj) = d_eq(24,1)
1071:         enddo
1072:         d_eq(24,8) = d_eq(24,8)                                         &
1073:      &             - k_eq(22) *an_r(24)*pt*pt /const2
1074:         d_eq(24,24) = d_eq(24,24)                                       &
1075:      &             - k_eq(22) *an_r(8)*pt*pt /const2
1076:         d_eq(24,3) = d_eq(24,3) + an_r(25)*pt*pt/const2
1077:         d_eq(24,25) = d_eq(24,25) + an_r(3)*pt*pt/const2

1079: !======== 
1080: !for 25 
1081:         
1082:         d_eq(25,1) = an_r(26)*(pt)*(-1/const2)                          &
1083:      &       - k_eq(23)*an_r(21)*an_r(10)*pt*pt * (-2/const_cube)
1084:         do jj = 2,26
1085:            d_eq(25,jj) = d_eq(25,1)
1086:         enddo
1087:         d_eq(25,10) = d_eq(25,10)                                       &
1088:      &             - k_eq(23) *an_r(21)*pt*pt /const2
1089:         d_eq(25,21) = d_eq(25,21)                                       &
1090:      &             - k_eq(23) *an_r(10)*pt*pt /const2
1091:         d_eq(25,26) = d_eq(25,26) + pt/(an_t)

1093: !============
1094: !  for 26 
1095:         d_eq(26,20) = -1
1096:         d_eq(26,22) = -1 
1097:         d_eq(26,23) = -1
1098:         d_eq(26,21) = 1    
1099:         d_eq(26,24) = 1    
1100:         d_eq(26,25) = 1   
1101:         d_eq(26,26) = 1


1104:            do j = 1,26
1105:          do i = 1,26
1106:                 write(44,*)i,j,d_eq(i,j)
1107:               enddo
1108:            enddo

1110:            
1111:         return
1112:         end

1114:       subroutine ShashiPostCheck(ls,X,Y,W,c_Y,c_W,dummy)
1115:  #include <petsc/finclude/petscsnes.h>
1116:       use petscsnes
1117:       implicit none
1118:       SNESLineSearch ls
1119:       PetscErrorCode ierr
1120:       Vec X,Y,W
1121:       PetscObject dummy
1122:       PetscBool c_Y,c_W
1123:       PetscScalar,pointer :: xx(:)
1124:       PetscInt i
1125:       call VecGetArrayF90(W,xx,ierr)
1126:       do i=1,26
1127:          if (xx(i) < 0.0) then
1128:             xx(i) = 0.0
1129:             c_W = PETSC_TRUE
1130:          endif
1131:         if (xx(i) > 3.0) then
1132:            xx(i) = 3.0
1133:         endif
1134:       enddo
1135:       call VecRestoreArrayF90(W,xx,ierr)
1136:       return
1137:       end
1138: