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: !     x, r        - solution, residual vectors
 14: !     J           - Jacobian matrix
 15: !
 16:       SNES     snes
 17:       Vec      x,r,lb,ub
 18:       Mat      J
 19:       PetscErrorCode  ierr
 20:       PetscInt i2
 21:       PetscMPIInt size
 22:       PetscScalar,pointer :: xx(:)
 23:       PetscScalar zero,big
 24:       SNESLineSearch ls

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

 29:       external FormFunction, FormJacobian
 30:       external ShashiPostCheck

 32: !
 33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 34: !                 Beginning of program
 35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 37:       PetscCallA(PetscInitialize(ierr))
 38:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
 39:       if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,1,'requires one process'); endif

 41:       big  = 2.88
 42:       big  = PETSC_INFINITY
 43:       zero = 0.0
 44:       i2  = 26
 45: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
 46: !  Create nonlinear solver context
 47: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

 49:       PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))

 51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52: !  Create matrix and vector data structures; set corresponding routines
 53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 55: !  Create vectors for solution and nonlinear function

 57:       PetscCallA(VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr))
 58:       PetscCallA(VecDuplicate(x,r,ierr))

 60: !  Create Jacobian matrix data structure

 62:       PetscCallA(MatCreateDense(PETSC_COMM_SELF,26,26,26,26,PETSC_NULL_SCALAR,J,ierr))

 64: !  Set function evaluation routine and vector

 66:       PetscCallA(SNESSetFunction(snes,r,FormFunction,0,ierr))

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

 70:       PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))

 72:       PetscCallA(VecDuplicate(x,lb,ierr))
 73:       PetscCallA(VecDuplicate(x,ub,ierr))
 74:       PetscCallA(VecSet(lb,zero,ierr))
 75: !      PetscCallA(VecGetArrayF90(lb,xx,ierr))
 76: !      PetscCallA(ShashiLowerBound(xx))
 77: !      PetscCallA(VecRestoreArrayF90(lb,xx,ierr))
 78:       PetscCallA(VecSet(ub,big,ierr))
 79: !      PetscCallA(SNESVISetVariableBounds(snes,lb,ub,ierr))

 81:       PetscCallA(SNESGetLineSearch(snes,ls,ierr))
 82:       PetscCallA(SNESLineSearchSetPostCheck(ls,ShashiPostCheck,0,ierr))
 83:       PetscCallA(SNESSetType(snes,SNESVINEWTONRSLS,ierr))

 85:       PetscCallA(SNESSetFromOptions(snes,ierr))

 87: !     set initial guess

 89:       PetscCallA(VecGetArrayF90(x,xx,ierr))
 90:       PetscCallA(ShashiInitialGuess(xx))
 91:       PetscCallA(VecRestoreArrayF90(x,xx,ierr))

 93:       PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))

 95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96: !  Free work space.  All PETSc objects should be destroyed when they
 97: !  are no longer needed.
 98: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 99:       PetscCallA(VecDestroy(lb,ierr))
100:       PetscCallA(VecDestroy(ub,ierr))
101:       PetscCallA(VecDestroy(x,ierr))
102:       PetscCallA(VecDestroy(r,ierr))
103:       PetscCallA(MatDestroy(J,ierr))
104:       PetscCallA(SNESDestroy(snes,ierr))
105:       PetscCallA(PetscFinalize(ierr))
106:       end
107: !
108: ! ------------------------------------------------------------------------
109: !
110: !  FormFunction - Evaluates nonlinear function, F(x).
111: !
112: !  Input Parameters:
113: !  snes - the SNES context
114: !  x - input vector
115: !  dummy - optional user-defined context (not used here)
116: !
117: !  Output Parameter:
118: !  f - function vector
119: !
120:       subroutine FormFunction(snes,x,f,dummy,ierr)
121: #include <petsc/finclude/petscsnes.h>
122:       use petscsnes
123:       implicit none
124:       SNES     snes
125:       Vec      x,f
126:       PetscErrorCode ierr
127:       integer dummy(*)

129: !  Declarations for use with local arrays

131:       PetscScalar,pointer ::lx_v(:),lf_v(:)

133: !  Get pointers to vector data.
134: !    - For default PETSc vectors, VecGetArray() returns a pointer to
135: !      the data array.  Otherwise, the routine is implementation dependent.
136: !    - You MUST call VecRestoreArray() when you no longer need access to
137: !      the array.
138: !    - Note that the Fortran interface to VecGetArray() differs from the
139: !      C version.  See the Fortran chapter of the users manual for details.

141:       PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
142:       PetscCall(VecGetArrayF90(f,lf_v,ierr))
143:       PetscCall(ShashiFormFunction(lx_v,lf_v))
144:       PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
145:       PetscCall(VecRestoreArrayF90(f,lf_v,ierr))

147:       return
148:       end

150: ! ---------------------------------------------------------------------
151: !
152: !  FormJacobian - Evaluates Jacobian matrix.
153: !
154: !  Input Parameters:
155: !  snes - the SNES context
156: !  x - input vector
157: !  dummy - optional user-defined context (not used here)
158: !
159: !  Output Parameters:
160: !  A - Jacobian matrix
161: !  B - optionally different preconditioning matrix
162: !  flag - flag indicating matrix structure
163: !
164:       subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
165: #include <petsc/finclude/petscsnes.h>
166:       use petscsnes
167:       implicit none
168:       SNES         snes
169:       Vec          X
170:       Mat          jac,B
171:       PetscErrorCode ierr
172:       integer dummy(*)

174: !  Declarations for use with local arrays
175:       PetscScalar,pointer ::lx_v(:),lf_v(:,:)

177: !  Get pointer to vector data

179:       PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
180:       PetscCall(MatDenseGetArrayF90(B,lf_v,ierr))
181:       PetscCall(ShashiFormJacobian(lx_v,lf_v))
182:       PetscCall(MatDenseRestoreArrayF90(B,lf_v,ierr))
183:       PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))

185: !  Assemble matrix

187:       PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
188:       PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))

190:       return
191:       end

193:             subroutine ShashiLowerBound(an_r)
194: !        implicit PetscScalar (a-h,o-z)
195:         implicit none
196:         PetscScalar an_r(26)
197:         PetscInt i

199:         do i=2,26
200:            an_r(i) = 1000.0/6.023D+23
201:         enddo
202:         return
203:         end

205:             subroutine ShashiInitialGuess(an_r)
206: !        implicit PetscScalar (a-h,o-z)
207:         implicit none
208:         PetscScalar an_c_additive
209:         PetscScalar       an_h_additive
210:         PetscScalar an_o_additive
211:         PetscScalar   atom_c_init
212:         PetscScalar atom_h_init
213:         PetscScalar atom_n_init
214:         PetscScalar atom_o_init
215:         PetscScalar h_init
216:         PetscScalar p_init
217:         PetscInt nfuel
218:         PetscScalar temp,pt
219:         PetscScalar an_r(26)
220:         PetscInt an_h(1),an_c(1)

222:         pt = 0.1
223:         atom_c_init =6.7408177364816552D-022
224:         atom_h_init = 2.0
225:         atom_o_init = 1.0
226:         atom_n_init = 3.76
227:         nfuel = 1
228:         an_c(1) = 1
229:         an_h(1) = 4
230:         an_c_additive = 2
231:         an_h_additive = 6
232:         an_o_additive = 1
233:         h_init = 128799.7267952987
234:         p_init = 0.1
235:         temp = 1500

237:       an_r( 1) =     1.66000D-24
238:       an_r( 2) =     1.66030D-22
239:       an_r( 3) =     5.00000D-01
240:       an_r( 4) =     1.66030D-22
241:       an_r( 5) =     1.66030D-22
242:       an_r( 6) =     1.88000D+00
243:       an_r( 7) =     1.66030D-22
244:       an_r( 8) =     1.66030D-22
245:       an_r( 9) =     1.66030D-22
246:       an_r(10) =     1.66030D-22
247:       an_r(11) =     1.66030D-22
248:       an_r(12) =     1.66030D-22
249:       an_r(13) =     1.66030D-22
250:       an_r(14) =     1.00000D+00
251:       an_r(15) =     1.66030D-22
252:       an_r(16) =     1.66030D-22
253:       an_r(17) =     1.66000D-24
254:       an_r(18) =     1.66030D-24
255:       an_r(19) =     1.66030D-24
256:       an_r(20) =     1.66030D-24
257:       an_r(21) =     1.66030D-24
258:       an_r(22) =     1.66030D-24
259:       an_r(23) =     1.66030D-24
260:       an_r(24) =     1.66030D-24
261:       an_r(25) =     1.66030D-24
262:       an_r(26) =     1.66030D-24

264:       an_r = 0
265:       an_r( 3) =     .5
266:       an_r( 6) =     1.88000
267:       an_r(14) =     1.

269: #if defined(solution)
270:       an_r( 1) =      3.802208D-33
271:       an_r( 2) =      1.298287D-29
272:       an_r( 3) =      2.533067D-04
273:       an_r( 4) =      6.865078D-22
274:       an_r( 5) =      9.993125D-01
275:       an_r( 6) =      1.879964D+00
276:       an_r( 7) =      4.449489D-13
277:       an_r( 8) =      3.428687D-07
278:       an_r( 9) =      7.105138D-05
279:       an_r(10) =      1.094368D-04
280:       an_r(11) =      2.362305D-06
281:       an_r(12) =      1.107145D-09
282:       an_r(13) =      1.276162D-24
283:       an_r(14) =      6.315538D-04
284:       an_r(15) =      2.356540D-09
285:       an_r(16) =      2.048248D-09
286:       an_r(17) =      1.966187D-22
287:       an_r(18) =      7.856497D-29
288:       an_r(19) =      1.987840D-36
289:       an_r(20) =      8.182441D-22
290:       an_r(21) =      2.684880D-16
291:       an_r(22) =      2.680473D-16
292:       an_r(23) =      6.594967D-18
293:       an_r(24) =      2.509714D-21
294:       an_r(25) =      3.096459D-21
295:       an_r(26) =      6.149551D-18
296: #endif

298:       return
299:       end

301:       subroutine ShashiFormFunction(an_r,f_eq)
302: !       implicit PetscScalar (a-h,o-z)
303:         implicit none
304:         PetscScalar an_c_additive
305:         PetscScalar       an_h_additive
306:         PetscScalar an_o_additive
307:         PetscScalar   atom_c_init
308:         PetscScalar atom_h_init
309:         PetscScalar atom_n_init
310:         PetscScalar atom_o_init
311:         PetscScalar h_init
312:         PetscScalar p_init
313:         PetscInt nfuel
314:         PetscScalar temp,pt
315:         PetscScalar an_r(26),k_eq(23),f_eq(26)
316:         PetscScalar H_molar(26)
317:         PetscInt an_h(1),an_c(1)
318:         PetscScalar part_p(26),idiff
319:         PetscInt i_cc,i_hh,i_h2o
320:         PetscScalar an_t,sum_h
321:         PetscScalar a_io2
322:         PetscInt i,ip
323:         pt = 0.1
324:         atom_c_init =6.7408177364816552D-022
325:         atom_h_init = 2.0
326:         atom_o_init = 1.0
327:         atom_n_init = 3.76
328:         nfuel = 1
329:         an_c(1) = 1
330:         an_h(1) = 4
331:         an_c_additive = 2
332:         an_h_additive = 6
333:         an_o_additive = 1
334:         h_init = 128799.7267952987
335:         p_init = 0.1
336:         temp = 1500

338:        k_eq( 1) =     1.75149D-05
339:        k_eq( 2) =     4.01405D-06
340:        k_eq( 3) =     6.04663D-14
341:        k_eq( 4) =     2.73612D-01
342:        k_eq( 5) =     3.25592D-03
343:        k_eq( 6) =     5.33568D+05
344:        k_eq( 7) =     2.07479D+05
345:        k_eq( 8) =     1.11841D-02
346:        k_eq( 9) =     1.72684D-03
347:        k_eq(10) =     1.98588D-07
348:        k_eq(11) =     7.23600D+27
349:        k_eq(12) =     5.73926D+49
350:        k_eq(13) =     1.00000D+00
351:        k_eq(14) =     1.64493D+16
352:        k_eq(15) =     2.73837D-29
353:        k_eq(16) =     3.27419D+50
354:        k_eq(17) =     1.72447D-23
355:        k_eq(18) =     4.24657D-06
356:        k_eq(19) =     1.16065D-14
357:        k_eq(20) =     3.28020D+25
358:        k_eq(21) =     1.06291D+00
359:        k_eq(22) =     9.11507D+02
360:        k_eq(23) =     6.02837D+03

362:        H_molar( 1) =     3.26044D+03
363:        H_molar( 2) =    -8.00407D+04
364:        H_molar( 3) =     4.05873D+04
365:        H_molar( 4) =    -3.31849D+05
366:        H_molar( 5) =    -1.93654D+05
367:        H_molar( 6) =     3.84035D+04
368:        H_molar( 7) =     4.97589D+05
369:        H_molar( 8) =     2.74483D+05
370:        H_molar( 9) =     1.30022D+05
371:        H_molar(10) =     7.58429D+04
372:        H_molar(11) =     2.42948D+05
373:        H_molar(12) =     1.44588D+05
374:        H_molar(13) =    -7.16891D+04
375:        H_molar(14) =     3.63075D+04
376:        H_molar(15) =     9.23880D+04
377:        H_molar(16) =     6.50477D+04
378:        H_molar(17) =     3.04310D+05
379:        H_molar(18) =     7.41707D+05
380:        H_molar(19) =     6.32767D+05
381:        H_molar(20) =     8.90624D+05
382:        H_molar(21) =     2.49805D+04
383:        H_molar(22) =     6.43473D+05
384:        H_molar(23) =     1.02861D+06
385:        H_molar(24) =    -6.07503D+03
386:        H_molar(25) =     1.27020D+05
387:        H_molar(26) =    -1.07011D+05
388: !=============
389:       an_t = 0
390:       sum_h = 0

392:       do i = 1,26
393:          an_t = an_t + an_r(i)
394:       enddo

396:         f_eq(1) = atom_h_init                                           &
397:      &          - (an_h(1)*an_r(1) + an_h_additive*an_r(2)              &
398:      &              + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14)      &
399:      &              + an_r(16)+ 2*an_r(17) + an_r(19)                   &
400:      &              +an_r(20) + 3*an_r(22)+an_r(26))

402:         f_eq(2) = atom_o_init                                           &
403:      &          - (an_o_additive*an_r(2) + 2*an_r(3)                    &
404:      &             + 2*an_r(4) + an_r(5)                                &
405:      &             + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
406:      &             + 2*an_r(15) + 2*an_r(16)+ an_r(20) + an_r(22)       &
407:      &             + an_r(23) + 2*an_r(24) + 1*an_r(25)+an_r(26))

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

411:         f_eq(4) = atom_c_init                                           &
412:      &          - (an_c(1)*an_r(1) + an_c_additive * an_r(2)            &
413:      &          + an_r(4) + an_r(13)+ 2*an_r(17) + an_r(18)             &
414:      &          + an_r(19) + an_r(20))

416:         do ip = 1,26
417:            part_p(ip) = (an_r(ip)/an_t) * pt
418:         enddo

420:         i_cc      = an_c(1)
421:         i_hh      = an_h(1)
422:         a_io2      = i_cc + i_hh/4.0
423:         i_h2o     = i_hh/2
424:         idiff     = (i_cc + i_h2o) - (a_io2 + 1)

426:         f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2                       &
427:      &          - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
428: !           write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
429: !          stop
430:         f_eq(6) = atom_n_init                                           &
431:      &          - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12)           &
432:      &              + an_r(15)                                          &
433:      &              + an_r(23))

435:       f_eq( 7) = part_p(11)                                             &
436:      &         - (k_eq( 1) * sqrt(part_p(14)+1d-23))
437:       f_eq( 8) = part_p( 8)                                             &
438:      &         - (k_eq( 2) * sqrt(part_p( 3)+1d-23))

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

443:       f_eq(10) = part_p(10)                                             &
444:      &         - (k_eq( 4) * sqrt(part_p( 3)+1d-23))                    &
445:      &         * sqrt(part_p(14))

447:       f_eq(11) = part_p( 9)                                             &
448:      &         - (k_eq( 5) * sqrt(part_p( 3)+1d-23))                    &
449:      &         * sqrt(part_p( 6)+1d-23)
450:       f_eq(12) = part_p( 5)                                             &
451:      &         - (k_eq( 6) * sqrt(part_p( 3)+1d-23))                    &
452:      &         * (part_p(14))

454:       f_eq(13) = part_p( 4)                                             &
455:      &         - (k_eq( 7) * sqrt(part_p(3)+1.0d-23))                   &
456:      &         * (part_p(13))

458:       f_eq(14) = part_p(15)                                             &
459:      &         - (k_eq( 8) * sqrt(part_p(3)+1.0d-50))                   &
460:      &         * (part_p( 9))
461:       f_eq(15) = part_p(16)                                             &
462:      &         - (k_eq( 9) * part_p( 3))                                &
463:      &         * sqrt(part_p(14)+1d-23)

465:       f_eq(16) = part_p(12)                                             &
466:      &         - (k_eq(10) * sqrt(part_p( 3)+1d-23))                    &
467:      &         * (part_p( 6))

469:       f_eq(17) = part_p(14)*part_p(18)**2                               &
470:      &         - (k_eq(15) * part_p(17))

472:       f_eq(18) = (part_p(13)**2)                                        &
473:      &     - (k_eq(16) * part_p(3)*part_p(18)**2)
474:       print*,f_eq(18),part_p(3),part_p(18),part_p(13),k_eq(16)

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

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

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

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

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

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

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

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

493:              do i = 1,26
494:                  write(44,*)i,f_eq(i)
495:               enddo

497:       return
498:       end

500:       subroutine ShashiFormJacobian(an_r,d_eq)
501: !        implicit PetscScalar (a-h,o-z)
502:         implicit none
503:         PetscScalar an_c_additive
504:         PetscScalar       an_h_additive
505:         PetscScalar an_o_additive
506:         PetscScalar   atom_c_init
507:         PetscScalar atom_h_init
508:         PetscScalar atom_n_init
509:         PetscScalar atom_o_init
510:         PetscScalar h_init
511:         PetscScalar p_init
512:         PetscInt nfuel
513:         PetscScalar temp,pt
514:         PetscScalar an_t,ai_o2
515:         PetscScalar an_tot1_d,an_tot1
516:         PetscScalar an_tot2_d,an_tot2
517:         PetscScalar const5,const3,const2
518:         PetscScalar const_cube
519:         PetscScalar const_five
520:         PetscScalar const_four
521:         PetscScalar const_six
522:         PetscInt jj,jb,ii3,id,ib,i
523:         PetscScalar pt2,pt1
524:         PetscScalar an_r(26),k_eq(23)
525:         PetscScalar d_eq(26,26),H_molar(26)
526:         PetscInt an_h(1),an_c(1)
527:         PetscScalar ai_pwr1,idiff
528:         PetscInt i_cc,i_hh
529:         PetscInt i_h2o,i_pwr2,i_co2_h2o
530:         PetscScalar pt_cube,pt_five
531:         PetscScalar pt_four
532:         PetscScalar pt_val1,pt_val2
533:         PetscInt j

535:         pt = 0.1
536:         atom_c_init =6.7408177364816552D-022
537:         atom_h_init = 2.0
538:         atom_o_init = 1.0
539:         atom_n_init = 3.76
540:         nfuel = 1
541:         an_c(1) = 1
542:         an_h(1) = 4
543:         an_c_additive = 2
544:         an_h_additive = 6
545:         an_o_additive = 1
546:         h_init = 128799.7267952987
547:         p_init = 0.1
548:         temp = 1500

550:        k_eq( 1) =     1.75149D-05
551:        k_eq( 2) =     4.01405D-06
552:        k_eq( 3) =     6.04663D-14
553:        k_eq( 4) =     2.73612D-01
554:        k_eq( 5) =     3.25592D-03
555:        k_eq( 6) =     5.33568D+05
556:        k_eq( 7) =     2.07479D+05
557:        k_eq( 8) =     1.11841D-02
558:        k_eq( 9) =     1.72684D-03
559:        k_eq(10) =     1.98588D-07
560:        k_eq(11) =     7.23600D+27
561:        k_eq(12) =     5.73926D+49
562:        k_eq(13) =     1.00000D+00
563:        k_eq(14) =     1.64493D+16
564:        k_eq(15) =     2.73837D-29
565:        k_eq(16) =     3.27419D+50
566:        k_eq(17) =     1.72447D-23
567:        k_eq(18) =     4.24657D-06
568:        k_eq(19) =     1.16065D-14
569:        k_eq(20) =     3.28020D+25
570:        k_eq(21) =     1.06291D+00
571:        k_eq(22) =     9.11507D+02
572:        k_eq(23) =     6.02837D+03

574:        H_molar( 1) =     3.26044D+03
575:        H_molar( 2) =    -8.00407D+04
576:        H_molar( 3) =     4.05873D+04
577:        H_molar( 4) =    -3.31849D+05
578:        H_molar( 5) =    -1.93654D+05
579:        H_molar( 6) =     3.84035D+04
580:        H_molar( 7) =     4.97589D+05
581:        H_molar( 8) =     2.74483D+05
582:        H_molar( 9) =     1.30022D+05
583:        H_molar(10) =     7.58429D+04
584:        H_molar(11) =     2.42948D+05
585:        H_molar(12) =     1.44588D+05
586:        H_molar(13) =    -7.16891D+04
587:        H_molar(14) =     3.63075D+04
588:        H_molar(15) =     9.23880D+04
589:        H_molar(16) =     6.50477D+04
590:        H_molar(17) =     3.04310D+05
591:        H_molar(18) =     7.41707D+05
592:        H_molar(19) =     6.32767D+05
593:        H_molar(20) =     8.90624D+05
594:        H_molar(21) =     2.49805D+04
595:        H_molar(22) =     6.43473D+05
596:        H_molar(23) =     1.02861D+06
597:        H_molar(24) =    -6.07503D+03
598:        H_molar(25) =     1.27020D+05
599:        H_molar(26) =    -1.07011D+05

601: !
602: !=======
603:       do jb = 1,26
604:          do ib = 1,26
605:             d_eq(ib,jb) = 0.0d0
606:          end do
607:       end do

609:         an_t = 0.0
610:       do id = 1,26
611:          an_t = an_t + an_r(id)
612:       enddo

614: !====
615: !====
616:         d_eq(1,1) = -an_h(1)
617:         d_eq(1,2) = -an_h_additive
618:         d_eq(1,5) = -2
619:         d_eq(1,10) = -1
620:         d_eq(1,11) = -1
621:         d_eq(1,14) = -2
622:         d_eq(1,16) = -1
623:         d_eq(1,17) = -2
624:         d_eq(1,19) = -1
625:         d_eq(1,20) = -1
626:         d_eq(1,22) = -3
627:         d_eq(1,26) = -1

629:         d_eq(2,2) = -1*an_o_additive
630:         d_eq(2,3) = -2
631:         d_eq(2,4) = -2
632:         d_eq(2,5) = -1
633:         d_eq(2,8) = -1
634:         d_eq(2,9) = -1
635:         d_eq(2,10) = -1
636:         d_eq(2,12) = -1
637:         d_eq(2,13) = -1
638:         d_eq(2,15) = -2
639:         d_eq(2,16) = -2
640:         d_eq(2,20) = -1
641:         d_eq(2,22) = -1
642:         d_eq(2,23) = -1
643:         d_eq(2,24) = -2
644:         d_eq(2,25) = -1
645:         d_eq(2,26) = -1

647:         d_eq(6,6) = -2
648:         d_eq(6,7) = -1
649:         d_eq(6,9) = -1
650:         d_eq(6,12) = -2
651:         d_eq(6,15) = -1
652:         d_eq(6,23) = -1

654:         d_eq(4,1) = -an_c(1)
655:         d_eq(4,2) = -an_c_additive
656:         d_eq(4,4) = -1
657:         d_eq(4,13) = -1
658:         d_eq(4,17) = -2
659:         d_eq(4,18) = -1
660:         d_eq(4,19) = -1
661:         d_eq(4,20) = -1

663: !----------
664:         const2 = an_t*an_t
665:         const3 = (an_t)*sqrt(an_t)
666:         const5 = an_t*const3

668:            const_cube =  an_t*an_t*an_t
669:            const_four =  const2*const2
670:            const_five =  const2*const_cube
671:            const_six  = const_cube*const_cube
672:            pt_cube = pt*pt*pt
673:            pt_four = pt_cube*pt
674:            pt_five = pt_cube*pt*pt

676:            i_cc = an_c(1)
677:            i_hh = an_h(1)
678:            ai_o2     = i_cc + float(i_hh)/4.0
679:            i_co2_h2o = i_cc + i_hh/2
680:            i_h2o     = i_hh/2
681:            ai_pwr1  = 1 + i_cc + float(i_hh)/4.0
682:            i_pwr2  = i_cc + i_hh/2
683:            idiff     = (i_cc + i_h2o) - (ai_o2 + 1)

685:            pt1   = pt**(ai_pwr1)
686:            an_tot1 = an_t**(ai_pwr1)
687:            pt_val1 = (pt/an_t)**(ai_pwr1)
688:            an_tot1_d = an_tot1*an_t

690:            pt2   = pt**(i_pwr2)
691:            an_tot2 = an_t**(i_pwr2)
692:            pt_val2 = (pt/an_t)**(i_pwr2)
693:            an_tot2_d = an_tot2*an_t

695:            d_eq(5,1) =                                                  &
696:      &           -(an_r(4)**i_cc)*(an_r(5)**i_h2o)                      &
697:      &           *((pt/an_t)**idiff) *(-idiff/an_t)

699:            do jj = 2,26
700:               d_eq(5,jj) = d_eq(5,1)
701:            enddo

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

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

708:            d_eq(5,4) = d_eq(5,4) - (i_cc*an_r(4)**(i_cc-1))*            &
709:      &                           (an_r(5)**i_h2o)* ((pt/an_t)**idiff)
710:            d_eq(5,5) = d_eq(5,5)                                        &
711:      &               - (i_h2o*(an_r(5)**(i_h2o-1)))                     &
712:      &               * (an_r(4)**i_cc)* ((pt/an_t)**idiff)

714:            d_eq(3,1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
715:            do jj = 2,26
716:               d_eq(3,jj) = d_eq(3,1)
717:            enddo

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

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

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

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

728:            do ii3 = 1,26
729:               d_eq(3,ii3) = 0.0d0
730:            enddo
731:            d_eq(3,2) = 1.0d0

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

736:         do jj = 2,26
737:            d_eq(7,jj) = d_eq(7,1)
738:         enddo

740:         d_eq(7,11) = d_eq(7,11) + pt/an_t
741:         d_eq(7,14) = d_eq(7,14)                                         &
742:      &            - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14)+1d-50)*an_t)))

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

747:         do jj = 2,26
748:            d_eq(8,jj) = d_eq(8,1)
749:         enddo

751:         d_eq(8,3) = d_eq(8,3)                                           &
752:      &            -k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3)+1.0d-50)*an_t)))
753:         d_eq(8,8) = d_eq(8,8) + pt/an_t

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

758:         do jj = 2,26
759:            d_eq(9,jj) = d_eq(9,1)
760:         enddo

762:         d_eq(9,7) = d_eq(9,7) + pt/an_t
763:         d_eq(9,6) = d_eq(9,6)                                           &
764:      &             -k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))

766:         d_eq(10,1) = pt*an_r(10)*(-1.0)/const2                          &
767:      &             -k_eq(4)*(pt)*sqrt((an_r(3)+1.0d-50)                 &
768:      &       *an_r(14))*(-1.0/const2)
769:         do jj = 2,26
770:            d_eq(10,jj) = d_eq(10,1)
771:         enddo

773:         d_eq(10,3) = d_eq(10,3)                                         &
774:      &           -k_eq(4)*(pt)*sqrt(an_r(14))                           &
775:      &           *(0.5/(sqrt(an_r(3)+1.0d-50)*an_t))
776:         d_eq(10,10) = d_eq(10,10) + pt/an_t
777:         d_eq(10,14) = d_eq(10,14)                                       &
778:      &           -k_eq(4)*(pt)*sqrt(an_r(3)+1.0d-50)                    &
779:      &            *(0.5/(sqrt(an_r(14)+1.0d-50)*an_t))

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

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

789:         d_eq(11,3) = d_eq(11,3)                                         &
790:      &            -k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/                     &
791:      &       (sqrt(an_r(3)+1.0d-50)*an_t))
792:         d_eq(11,6) = d_eq(11,6)                                         &
793:      &            -k_eq(5)*(pt)*sqrt(an_r(3)+1.0d-50)                   &
794:      &       *(0.5/(sqrt(an_r(6))*an_t))
795:         d_eq(11,9) = d_eq(11,9) + pt/an_t

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

801:         do jj = 2,26
802:            d_eq(12,jj) = d_eq(12,1)
803:         enddo

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

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

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

817:         do jj = 2,26
818:            d_eq(13,jj) = d_eq(13,1)
819:         enddo

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

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

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

833:         do jj = 2,26
834:            d_eq(14,jj) = d_eq(14,1)
835:         enddo

837:         d_eq(14,3) = d_eq(14,3)                                         &
838:      &            -k_eq(8)*(pt**1.5)*(an_r(9)/const3)                   &
839:      &            *(0.5/sqrt(an_r(3)+1.0d-50))
840:         d_eq(14,9) = d_eq(14,9)                                         &
841:      &            -k_eq(8)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
842:         d_eq(14,15) = d_eq(14,15)+ pt/an_t

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

848:         do jj = 2,26
849:            d_eq(15,jj) = d_eq(15,1)
850:         enddo

852:         d_eq(15,3) = d_eq(15,3)                                         &
853:      &            -k_eq(9)*(pt**1.5)*(sqrt(an_r(14)+1.0d-50)/const3)
854:         d_eq(15,14) = d_eq(15,14)                                       &
855:      &            -k_eq(9)*(pt**1.5)*(an_r(3)/const3)                   &
856:      &            *(0.5/sqrt(an_r(14)+1.0d-50))
857:         d_eq(15,16) = d_eq(15,16) + pt/an_t

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

863:         do jj = 2,26
864:            d_eq(16,jj) = d_eq(16,1)
865:         enddo

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

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

875:         const_cube =  an_t*an_t*an_t
876:         const_four =  const2*const2

878:         d_eq(17,1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
879:      &             - k_eq(15) * an_r(17)*pt * (-1/const2)
880:         do jj = 2,26
881:            d_eq(17,jj) = d_eq(17,1)
882:         enddo
883:         d_eq(17,14) = d_eq(17,14) + an_r(18)*an_r(18)*(pt**3)/const_cube
884:         d_eq(17,17) = d_eq(17,17) - k_eq(15)*pt/an_t
885:         d_eq(17,18) = d_eq(17,18) + 2*an_r(18)*an_r(14)                 &
886:      &                            *(pt**3)/const_cube

888:         d_eq(18,1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube)          &
889:      &             - k_eq(16) * an_r(3)*an_r(18)*an_r(18)               &
890:      &              * (pt*pt*pt) * (-3/const_four)
891:         do jj = 2,26
892:            d_eq(18,jj) = d_eq(18,1)
893:         enddo
894:         d_eq(18,3) = d_eq(18,3)                                         &
895:      &             - k_eq(16) *an_r(18)* an_r(18)*pt*pt*pt /const_cube
896:         d_eq(18,13) = d_eq(18,13)                                       &
897:      &              + 2* an_r(13)*pt*pt /const2
898:         d_eq(18,18) = d_eq(18,18) -k_eq(16)*an_r(3)                     &
899:      &              * 2*an_r(18)*pt*pt*pt/const_cube

901: !====for eq 19

903:         d_eq(19,1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube)           &
904:      &             - k_eq(17)*an_r(13)*an_r(10)*pt*pt * (-2/const_cube)
905:         do jj = 2,26
906:            d_eq(19,jj) = d_eq(19,1)
907:         enddo
908:         d_eq(19,13) = d_eq(19,13)                                       &
909:      &             - k_eq(17) *an_r(10)*pt*pt /const2
910:         d_eq(19,10) = d_eq(19,10)                                       &
911:      &             - k_eq(17) *an_r(13)*pt*pt /const2
912:         d_eq(19,3) = d_eq(19,3) + an_r(19)*pt*pt/const2
913:         d_eq(19,19) = d_eq(19,19) + an_r(3)*pt*pt/const2
914: !====for eq 20

916:         d_eq(20,1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube)          &
917:      &             - k_eq(18) * an_r(19)*an_r(8)*pt*pt * (-2/const_cube)
918:         do jj = 2,26
919:            d_eq(20,jj) = d_eq(20,1)
920:         enddo
921:         d_eq(20,8) = d_eq(20,8)                                         &
922:      &             - k_eq(18) *an_r(19)*pt*pt /const2
923:         d_eq(20,19) = d_eq(20,19)                                       &
924:      &             - k_eq(18) *an_r(8)*pt*pt /const2
925:         d_eq(20,20) = d_eq(20,20) + an_r(21)*pt*pt/const2
926:         d_eq(20,21) = d_eq(20,21) + an_r(20)*pt*pt/const2

928: !========
929: !====for eq 21

931:         d_eq(21,1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube)          &
932:      &             - k_eq(19)*an_r(7)*an_r(8)*pt*pt * (-2/const_cube)
933:         do jj = 2,26
934:            d_eq(21,jj) = d_eq(21,1)
935:         enddo
936:         d_eq(21,7) = d_eq(21,7)                                         &
937:      &             - k_eq(19) *an_r(8)*pt*pt /const2
938:         d_eq(21,8) = d_eq(21,8)                                         &
939:      &             - k_eq(19) *an_r(7)*pt*pt /const2
940:         d_eq(21,21) = d_eq(21,21) + an_r(23)*pt*pt/const2
941:         d_eq(21,23) = d_eq(21,23) + an_r(21)*pt*pt/const2

943: !========
944: !  for 22
945:         d_eq(22,1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube)           &
946:      &         -k_eq(20)*an_r(21)*an_r(22)*pt*pt * (-2/const_cube)
947:         do jj = 2,26
948:            d_eq(22,jj) = d_eq(22,1)
949:         enddo
950:         d_eq(22,21) = d_eq(22,21)                                       &
951:      &             - k_eq(20) *an_r(22)*pt*pt /const2
952:         d_eq(22,22) = d_eq(22,22)                                       &
953:      &             - k_eq(20) *an_r(21)*pt*pt /const2
954:         d_eq(22,11) = d_eq(22,11) + an_r(5)*pt*pt/(const2)
955:         d_eq(22,5) = d_eq(22,5) + an_r(11)*pt*pt/(const2)

957: !========
958: !  for 23

960:         d_eq(23,1) = an_r(24)*(pt)*(-1/const2)                          &
961:      &             - k_eq(21)*an_r(21)*an_r(3)*pt*pt * (-2/const_cube)
962:         do jj = 2,26
963:            d_eq(23,jj) = d_eq(23,1)
964:         enddo
965:         d_eq(23,3) = d_eq(23,3)                                         &
966:      &             - k_eq(21) *an_r(21)*pt*pt /const2
967:         d_eq(23,21) = d_eq(23,21)                                       &
968:      &             - k_eq(21) *an_r(3)*pt*pt /const2
969:         d_eq(23,24) = d_eq(23,24) + pt/(an_t)

971: !========
972: !  for 24
973:         d_eq(24,1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube)           &
974:      &             - k_eq(22)*an_r(24)*an_r(8)*pt*pt * (-2/const_cube)
975:         do jj = 2,26
976:            d_eq(24,jj) = d_eq(24,1)
977:         enddo
978:         d_eq(24,8) = d_eq(24,8)                                         &
979:      &             - k_eq(22) *an_r(24)*pt*pt /const2
980:         d_eq(24,24) = d_eq(24,24)                                       &
981:      &             - k_eq(22) *an_r(8)*pt*pt /const2
982:         d_eq(24,3) = d_eq(24,3) + an_r(25)*pt*pt/const2
983:         d_eq(24,25) = d_eq(24,25) + an_r(3)*pt*pt/const2

985: !========
986: !for 25

988:         d_eq(25,1) = an_r(26)*(pt)*(-1/const2)                          &
989:      &       - k_eq(23)*an_r(21)*an_r(10)*pt*pt * (-2/const_cube)
990:         do jj = 2,26
991:            d_eq(25,jj) = d_eq(25,1)
992:         enddo
993:         d_eq(25,10) = d_eq(25,10)                                       &
994:      &             - k_eq(23) *an_r(21)*pt*pt /const2
995:         d_eq(25,21) = d_eq(25,21)                                       &
996:      &             - k_eq(23) *an_r(10)*pt*pt /const2
997:         d_eq(25,26) = d_eq(25,26) + pt/(an_t)

999: !============
1000: !  for 26
1001:         d_eq(26,20) = -1
1002:         d_eq(26,22) = -1
1003:         d_eq(26,23) = -1
1004:         d_eq(26,21) = 1
1005:         d_eq(26,24) = 1
1006:         d_eq(26,25) = 1
1007:         d_eq(26,26) = 1

1009:            do j = 1,26
1010:          do i = 1,26
1011:                 write(44,*)i,j,d_eq(i,j)
1012:               enddo
1013:            enddo

1015:         return
1016:         end

1018:       subroutine ShashiPostCheck(ls,X,Y,W,c_Y,c_W,dummy)
1019: #include <petsc/finclude/petscsnes.h>
1020:       use petscsnes
1021:       implicit none
1022:       SNESLineSearch ls
1023:       PetscErrorCode ierr
1024:       Vec X,Y,W
1025:       PetscObject dummy
1026:       PetscBool c_Y,c_W
1027:       PetscScalar,pointer :: xx(:)
1028:       PetscInt i
1029:       PetscCall(VecGetArrayF90(W,xx,ierr))
1030:       do i=1,26
1031:          if (xx(i) < 0.0) then
1032:             xx(i) = 0.0
1033:             c_W = PETSC_TRUE
1034:          endif
1035:         if (xx(i) > 3.0) then
1036:            xx(i) = 3.0
1037:         endif
1038:       enddo
1039:       PetscCall(VecRestoreArrayF90(W,xx,ierr))
1040:       return
1041:       end