Actual source code: shashi.F90
1: program main
2: #include <petsc/finclude/petsc.h>
3: use petsc
4: implicit none
6: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
7: ! Variable declarations
8: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: !
10: ! Variables:
11: ! snes - nonlinear solver
12: ! x, r - solution, residual vectors
13: ! J - Jacobian matrix
14: !
15: SNES snes
16: Vec x, r, lb, ub
17: Mat J
18: PetscErrorCode ierr
19: PetscInt i2
20: PetscMPIInt size
21: PetscScalar, pointer :: xx(:)
22: PetscScalar zero, big
23: SNESLineSearch ls
25: ! Note: Any user-defined Fortran routines (such as FormJacobian)
26: ! MUST be declared as external.
28: external FormFunction, FormJacobian
29: external ShashiPostCheck
31: !
32: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
33: ! Beginning of program
34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: PetscCallA(PetscInitialize(ierr))
37: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
38: PetscCheckA(size == 1, PETSC_COMM_WORLD, 1, 'requires one process')
40: big = 2.88
41: big = PETSC_INFINITY
42: zero = 0.0
43: i2 = 26
44: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Create nonlinear solver context
46: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
48: PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Create matrix and vector data structures; set corresponding routines
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: ! Create vectors for solution and nonlinear function
56: PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i2, x, ierr))
57: PetscCallA(VecDuplicate(x, r, ierr))
59: ! Create Jacobian matrix data structure
61: PetscCallA(MatCreateDense(PETSC_COMM_SELF, 26, 26, 26, 26, PETSC_NULL_SCALAR, J, ierr))
63: ! Set function evaluation routine and vector
65: PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr))
67: ! Set Jacobian matrix data structure and Jacobian evaluation routine
69: PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr))
71: PetscCallA(VecDuplicate(x, lb, ierr))
72: PetscCallA(VecDuplicate(x, ub, ierr))
73: PetscCallA(VecSet(lb, zero, ierr))
74: ! PetscCallA(VecGetArray(lb,xx,ierr))
75: ! PetscCallA(ShashiLowerBound(xx))
76: ! PetscCallA(VecRestoreArray(lb,xx,ierr))
77: PetscCallA(VecSet(ub, big, ierr))
78: ! PetscCallA(SNESVISetVariableBounds(snes,lb,ub,ierr))
80: PetscCallA(SNESGetLineSearch(snes, ls, ierr))
81: PetscCallA(SNESLineSearchSetPostCheck(ls, ShashiPostCheck, 0, ierr))
82: PetscCallA(SNESSetType(snes, SNESVINEWTONRSLS, ierr))
84: PetscCallA(SNESSetFromOptions(snes, ierr))
86: ! set initial guess
88: PetscCallA(VecGetArray(x, xx, ierr))
89: PetscCallA(ShashiInitialGuess(xx))
90: PetscCallA(VecRestoreArray(x, xx, ierr))
92: PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: ! Free work space. All PETSc objects should be destroyed when they
96: ! are no longer needed.
97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: PetscCallA(VecDestroy(lb, ierr))
99: PetscCallA(VecDestroy(ub, ierr))
100: PetscCallA(VecDestroy(x, ierr))
101: PetscCallA(VecDestroy(r, ierr))
102: PetscCallA(MatDestroy(J, ierr))
103: PetscCallA(SNESDestroy(snes, ierr))
104: PetscCallA(PetscFinalize(ierr))
105: end
106: !
107: ! ------------------------------------------------------------------------
108: !
109: ! FormFunction - Evaluates nonlinear function, F(x).
110: !
111: ! Input Parameters:
112: ! snes - the SNES context
113: ! x - input vector
114: ! dummy - optional user-defined context (not used here)
115: !
116: ! Output Parameter:
117: ! f - function vector
118: !
119: subroutine FormFunction(snes, x, f, dummy, ierr)
120: #include <petsc/finclude/petscsnes.h>
121: use petscsnes
122: implicit none
123: SNES snes
124: Vec x, f
125: PetscErrorCode ierr
126: integer dummy(*)
128: ! Declarations for use with local arrays
130: PetscScalar, pointer ::lx_v(:), lf_v(:)
132: ! Get pointers to vector data.
133: ! - For default PETSc vectors, VecGetArray() returns a pointer to
134: ! the data array. Otherwise, the routine is implementation dependent.
135: ! - You MUST call VecRestoreArray() when you no longer need access to
136: ! the array.
137: ! - Note that the Fortran interface to VecGetArray() differs from the
138: ! C version. See the Fortran chapter of the users manual for details.
140: PetscCall(VecGetArrayRead(x, lx_v, ierr))
141: PetscCall(VecGetArray(f, lf_v, ierr))
142: PetscCall(ShashiFormFunction(lx_v, lf_v))
143: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
144: PetscCall(VecRestoreArray(f, lf_v, ierr))
146: end
148: ! ---------------------------------------------------------------------
149: !
150: ! FormJacobian - Evaluates Jacobian matrix.
151: !
152: ! Input Parameters:
153: ! snes - the SNES context
154: ! x - input vector
155: ! dummy - optional user-defined context (not used here)
156: !
157: ! Output Parameters:
158: ! A - Jacobian matrix
159: ! B - optionally different matrix used to construct the preconditioner
160: !
161: subroutine FormJacobian(snes, X, jac, B, dummy, ierr)
162: #include <petsc/finclude/petscsnes.h>
163: use petscsnes
164: implicit none
165: SNES snes
166: Vec X
167: Mat jac, B
168: PetscErrorCode ierr
169: integer dummy(*)
171: ! Declarations for use with local arrays
172: PetscScalar, pointer ::lx_v(:), lf_v(:, :)
174: ! Get pointer to vector data
176: PetscCall(VecGetArrayRead(x, lx_v, ierr))
177: PetscCall(MatDenseGetArray(B, lf_v, ierr))
178: PetscCall(ShashiFormJacobian(lx_v, lf_v))
179: PetscCall(MatDenseRestoreArray(B, lf_v, ierr))
180: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
182: ! Assemble matrix
184: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
185: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
187: end
189: subroutine ShashiLowerBound(an_r)
190: ! implicit PetscScalar (a-h,o-z)
191: implicit none
192: PetscScalar an_r(26)
193: PetscInt i
195: do i = 2, 26
196: an_r(i) = 1000.0/6.023D+23
197: end do
198: end
200: subroutine ShashiInitialGuess(an_r)
201: ! implicit PetscScalar (a-h,o-z)
202: implicit none
203: PetscScalar an_c_additive
204: PetscScalar an_h_additive
205: PetscScalar an_o_additive
206: PetscScalar atom_c_init
207: PetscScalar atom_h_init
208: PetscScalar atom_n_init
209: PetscScalar atom_o_init
210: PetscScalar h_init
211: PetscScalar p_init
212: PetscInt nfuel
213: PetscScalar temp, pt
214: PetscScalar an_r(26)
215: PetscInt an_h(1), an_c(1)
217: pt = 0.1
218: atom_c_init = 6.7408177364816552D-022
219: atom_h_init = 2.0
220: atom_o_init = 1.0
221: atom_n_init = 3.76
222: nfuel = 1
223: an_c(1) = 1
224: an_h(1) = 4
225: an_c_additive = 2
226: an_h_additive = 6
227: an_o_additive = 1
228: h_init = 128799.7267952987
229: p_init = 0.1
230: temp = 1500
232: an_r(1) = 1.66000D-24
233: an_r(2) = 1.66030D-22
234: an_r(3) = 5.00000D-01
235: an_r(4) = 1.66030D-22
236: an_r(5) = 1.66030D-22
237: an_r(6) = 1.88000D+00
238: an_r(7) = 1.66030D-22
239: an_r(8) = 1.66030D-22
240: an_r(9) = 1.66030D-22
241: an_r(10) = 1.66030D-22
242: an_r(11) = 1.66030D-22
243: an_r(12) = 1.66030D-22
244: an_r(13) = 1.66030D-22
245: an_r(14) = 1.00000D+00
246: an_r(15) = 1.66030D-22
247: an_r(16) = 1.66030D-22
248: an_r(17) = 1.66000D-24
249: an_r(18) = 1.66030D-24
250: an_r(19) = 1.66030D-24
251: an_r(20) = 1.66030D-24
252: an_r(21) = 1.66030D-24
253: an_r(22) = 1.66030D-24
254: an_r(23) = 1.66030D-24
255: an_r(24) = 1.66030D-24
256: an_r(25) = 1.66030D-24
257: an_r(26) = 1.66030D-24
259: an_r = 0
260: an_r(3) = .5
261: an_r(6) = 1.88000
262: an_r(14) = 1.
264: #if defined(solution)
265: an_r(1) = 3.802208D-33
266: an_r(2) = 1.298287D-29
267: an_r(3) = 2.533067D-04
268: an_r(4) = 6.865078D-22
269: an_r(5) = 9.993125D-01
270: an_r(6) = 1.879964D+00
271: an_r(7) = 4.449489D-13
272: an_r(8) = 3.428687D-07
273: an_r(9) = 7.105138D-05
274: an_r(10) = 1.094368D-04
275: an_r(11) = 2.362305D-06
276: an_r(12) = 1.107145D-09
277: an_r(13) = 1.276162D-24
278: an_r(14) = 6.315538D-04
279: an_r(15) = 2.356540D-09
280: an_r(16) = 2.048248D-09
281: an_r(17) = 1.966187D-22
282: an_r(18) = 7.856497D-29
283: an_r(19) = 1.987840D-36
284: an_r(20) = 8.182441D-22
285: an_r(21) = 2.684880D-16
286: an_r(22) = 2.680473D-16
287: an_r(23) = 6.594967D-18
288: an_r(24) = 2.509714D-21
289: an_r(25) = 3.096459D-21
290: an_r(26) = 6.149551D-18
291: #endif
293: end
295: subroutine ShashiFormFunction(an_r, f_eq)
296: ! implicit PetscScalar (a-h,o-z)
297: implicit none
298: PetscScalar an_c_additive
299: PetscScalar an_h_additive
300: PetscScalar an_o_additive
301: PetscScalar atom_c_init
302: PetscScalar atom_h_init
303: PetscScalar atom_n_init
304: PetscScalar atom_o_init
305: PetscScalar h_init
306: PetscScalar p_init
307: PetscInt nfuel
308: PetscScalar temp, pt
309: PetscScalar an_r(26), k_eq(23), f_eq(26)
310: PetscScalar H_molar(26)
311: PetscInt an_h(1), an_c(1)
312: PetscScalar part_p(26), idiff
313: PetscInt i_cc, i_hh, i_h2o
314: PetscScalar an_t, sum_h
315: PetscScalar a_io2
316: PetscInt i, ip
317: pt = 0.1
318: atom_c_init = 6.7408177364816552D-022
319: atom_h_init = 2.0
320: atom_o_init = 1.0
321: atom_n_init = 3.76
322: nfuel = 1
323: an_c(1) = 1
324: an_h(1) = 4
325: an_c_additive = 2
326: an_h_additive = 6
327: an_o_additive = 1
328: h_init = 128799.7267952987
329: p_init = 0.1
330: temp = 1500
332: k_eq(1) = 1.75149D-05
333: k_eq(2) = 4.01405D-06
334: k_eq(3) = 6.04663D-14
335: k_eq(4) = 2.73612D-01
336: k_eq(5) = 3.25592D-03
337: k_eq(6) = 5.33568D+05
338: k_eq(7) = 2.07479D+05
339: k_eq(8) = 1.11841D-02
340: k_eq(9) = 1.72684D-03
341: k_eq(10) = 1.98588D-07
342: k_eq(11) = 7.23600D+27
343: k_eq(12) = 5.73926D+49
344: k_eq(13) = 1.00000D+00
345: k_eq(14) = 1.64493D+16
346: k_eq(15) = 2.73837D-29
347: k_eq(16) = 3.27419D+50
348: k_eq(17) = 1.72447D-23
349: k_eq(18) = 4.24657D-06
350: k_eq(19) = 1.16065D-14
351: k_eq(20) = 3.28020D+25
352: k_eq(21) = 1.06291D+00
353: k_eq(22) = 9.11507D+02
354: k_eq(23) = 6.02837D+03
356: H_molar(1) = 3.26044D+03
357: H_molar(2) = -8.00407D+04
358: H_molar(3) = 4.05873D+04
359: H_molar(4) = -3.31849D+05
360: H_molar(5) = -1.93654D+05
361: H_molar(6) = 3.84035D+04
362: H_molar(7) = 4.97589D+05
363: H_molar(8) = 2.74483D+05
364: H_molar(9) = 1.30022D+05
365: H_molar(10) = 7.58429D+04
366: H_molar(11) = 2.42948D+05
367: H_molar(12) = 1.44588D+05
368: H_molar(13) = -7.16891D+04
369: H_molar(14) = 3.63075D+04
370: H_molar(15) = 9.23880D+04
371: H_molar(16) = 6.50477D+04
372: H_molar(17) = 3.04310D+05
373: H_molar(18) = 7.41707D+05
374: H_molar(19) = 6.32767D+05
375: H_molar(20) = 8.90624D+05
376: H_molar(21) = 2.49805D+04
377: H_molar(22) = 6.43473D+05
378: H_molar(23) = 1.02861D+06
379: H_molar(24) = -6.07503D+03
380: H_molar(25) = 1.27020D+05
381: H_molar(26) = -1.07011D+05
382: !=============
383: an_t = 0
384: sum_h = 0
386: do i = 1, 26
387: an_t = an_t + an_r(i)
388: end do
390: f_eq(1) = atom_h_init &
391: & - (an_h(1)*an_r(1) + an_h_additive*an_r(2) &
392: & + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) &
393: & + an_r(16) + 2*an_r(17) + an_r(19) &
394: & + an_r(20) + 3*an_r(22) + an_r(26))
396: f_eq(2) = atom_o_init &
397: & - (an_o_additive*an_r(2) + 2*an_r(3) &
398: & + 2*an_r(4) + an_r(5) &
399: & + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
400: & + 2*an_r(15) + 2*an_r(16) + an_r(20) + an_r(22) &
401: & + an_r(23) + 2*an_r(24) + 1*an_r(25) + an_r(26))
403: f_eq(3) = an_r(2) - 1.0d-150
405: f_eq(4) = atom_c_init &
406: & - (an_c(1)*an_r(1) + an_c_additive*an_r(2) &
407: & + an_r(4) + an_r(13) + 2*an_r(17) + an_r(18) &
408: & + an_r(19) + an_r(20))
410: do ip = 1, 26
411: part_p(ip) = (an_r(ip)/an_t)*pt
412: end do
414: i_cc = an_c(1)
415: i_hh = an_h(1)
416: a_io2 = i_cc + i_hh/4.0
417: i_h2o = i_hh/2
418: idiff = (i_cc + i_h2o) - (a_io2 + 1)
420: f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2 &
421: & - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
422: ! write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
423: ! stop
424: f_eq(6) = atom_n_init &
425: & - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) &
426: & + an_r(15) &
427: & + an_r(23))
429: f_eq(7) = part_p(11) &
430: & - (k_eq(1)*sqrt(part_p(14) + 1d-23))
431: f_eq(8) = part_p(8) &
432: & - (k_eq(2)*sqrt(part_p(3) + 1d-23))
434: f_eq(9) = part_p(7) &
435: & - (k_eq(3)*sqrt(part_p(6) + 1d-23))
437: f_eq(10) = part_p(10) &
438: & - (k_eq(4)*sqrt(part_p(3) + 1d-23)) &
439: & *sqrt(part_p(14))
441: f_eq(11) = part_p(9) &
442: & - (k_eq(5)*sqrt(part_p(3) + 1d-23)) &
443: & *sqrt(part_p(6) + 1d-23)
444: f_eq(12) = part_p(5) &
445: & - (k_eq(6)*sqrt(part_p(3) + 1d-23)) &
446: & *(part_p(14))
448: f_eq(13) = part_p(4) &
449: & - (k_eq(7)*sqrt(part_p(3) + 1.0d-23)) &
450: & *(part_p(13))
452: f_eq(14) = part_p(15) &
453: & - (k_eq(8)*sqrt(part_p(3) + 1.0d-50)) &
454: & *(part_p(9))
455: f_eq(15) = part_p(16) &
456: & - (k_eq(9)*part_p(3)) &
457: & *sqrt(part_p(14) + 1d-23)
459: f_eq(16) = part_p(12) &
460: & - (k_eq(10)*sqrt(part_p(3) + 1d-23)) &
461: & *(part_p(6))
463: f_eq(17) = part_p(14)*part_p(18)**2 &
464: & - (k_eq(15)*part_p(17))
466: f_eq(18) = (part_p(13)**2) &
467: & - (k_eq(16)*part_p(3)*part_p(18)**2)
468: print *, f_eq(18), part_p(3), part_p(18), part_p(13), k_eq(16)
470: f_eq(19) = part_p(19)*part_p(3) - k_eq(17)*part_p(13)*part_p(10)
472: f_eq(20) = part_p(21)*part_p(20) - k_eq(18)*part_p(19)*part_p(8)
474: f_eq(21) = part_p(21)*part_p(23) - k_eq(19)*part_p(7)*part_p(8)
476: f_eq(22) = part_p(5)*part_p(11) - k_eq(20)*part_p(21)*part_p(22)
478: f_eq(23) = part_p(24) - k_eq(21)*part_p(21)*part_p(3)
480: f_eq(24) = part_p(3)*part_p(25) - k_eq(22)*part_p(24)*part_p(8)
482: f_eq(25) = part_p(26) - k_eq(23)*part_p(21)*part_p(10)
484: f_eq(26) = -(an_r(20) + an_r(22) + an_r(23)) &
485: & + (an_r(21) + an_r(24) + an_r(25) + an_r(26))
487: do i = 1, 26
488: write (44, *) i, f_eq(i)
489: end do
491: end
493: subroutine ShashiFormJacobian(an_r, d_eq)
494: ! implicit PetscScalar (a-h,o-z)
495: implicit none
496: PetscScalar an_c_additive
497: PetscScalar an_h_additive
498: PetscScalar an_o_additive
499: PetscScalar atom_c_init
500: PetscScalar atom_h_init
501: PetscScalar atom_n_init
502: PetscScalar atom_o_init
503: PetscScalar h_init
504: PetscScalar p_init
505: PetscInt nfuel
506: PetscScalar temp, pt
507: PetscScalar an_t, ai_o2
508: PetscScalar an_tot1_d, an_tot1
509: PetscScalar an_tot2_d, an_tot2
510: PetscScalar const5, const3, const2
511: PetscScalar const_cube
512: PetscScalar const_five
513: PetscScalar const_four
514: PetscScalar const_six
515: PetscInt jj, jb, ii3, id, ib, i
516: PetscScalar pt2, pt1
517: PetscScalar an_r(26), k_eq(23)
518: PetscScalar d_eq(26, 26), H_molar(26)
519: PetscInt an_h(1), an_c(1)
520: PetscScalar ai_pwr1, idiff
521: PetscInt i_cc, i_hh
522: PetscInt i_h2o, i_pwr2, i_co2_h2o
523: PetscScalar pt_cube, pt_five
524: PetscScalar pt_four
525: PetscScalar pt_val1, pt_val2
526: PetscInt j
528: pt = 0.1
529: atom_c_init = 6.7408177364816552D-022
530: atom_h_init = 2.0
531: atom_o_init = 1.0
532: atom_n_init = 3.76
533: nfuel = 1
534: an_c(1) = 1
535: an_h(1) = 4
536: an_c_additive = 2
537: an_h_additive = 6
538: an_o_additive = 1
539: h_init = 128799.7267952987
540: p_init = 0.1
541: temp = 1500
543: k_eq(1) = 1.75149D-05
544: k_eq(2) = 4.01405D-06
545: k_eq(3) = 6.04663D-14
546: k_eq(4) = 2.73612D-01
547: k_eq(5) = 3.25592D-03
548: k_eq(6) = 5.33568D+05
549: k_eq(7) = 2.07479D+05
550: k_eq(8) = 1.11841D-02
551: k_eq(9) = 1.72684D-03
552: k_eq(10) = 1.98588D-07
553: k_eq(11) = 7.23600D+27
554: k_eq(12) = 5.73926D+49
555: k_eq(13) = 1.00000D+00
556: k_eq(14) = 1.64493D+16
557: k_eq(15) = 2.73837D-29
558: k_eq(16) = 3.27419D+50
559: k_eq(17) = 1.72447D-23
560: k_eq(18) = 4.24657D-06
561: k_eq(19) = 1.16065D-14
562: k_eq(20) = 3.28020D+25
563: k_eq(21) = 1.06291D+00
564: k_eq(22) = 9.11507D+02
565: k_eq(23) = 6.02837D+03
567: H_molar(1) = 3.26044D+03
568: H_molar(2) = -8.00407D+04
569: H_molar(3) = 4.05873D+04
570: H_molar(4) = -3.31849D+05
571: H_molar(5) = -1.93654D+05
572: H_molar(6) = 3.84035D+04
573: H_molar(7) = 4.97589D+05
574: H_molar(8) = 2.74483D+05
575: H_molar(9) = 1.30022D+05
576: H_molar(10) = 7.58429D+04
577: H_molar(11) = 2.42948D+05
578: H_molar(12) = 1.44588D+05
579: H_molar(13) = -7.16891D+04
580: H_molar(14) = 3.63075D+04
581: H_molar(15) = 9.23880D+04
582: H_molar(16) = 6.50477D+04
583: H_molar(17) = 3.04310D+05
584: H_molar(18) = 7.41707D+05
585: H_molar(19) = 6.32767D+05
586: H_molar(20) = 8.90624D+05
587: H_molar(21) = 2.49805D+04
588: H_molar(22) = 6.43473D+05
589: H_molar(23) = 1.02861D+06
590: H_molar(24) = -6.07503D+03
591: H_molar(25) = 1.27020D+05
592: H_molar(26) = -1.07011D+05
594: !
595: !=======
596: do jb = 1, 26
597: do ib = 1, 26
598: d_eq(ib, jb) = 0.0d0
599: end do
600: end do
602: an_t = 0.0
603: do id = 1, 26
604: an_t = an_t + an_r(id)
605: end do
607: !====
608: !====
609: d_eq(1, 1) = -an_h(1)
610: d_eq(1, 2) = -an_h_additive
611: d_eq(1, 5) = -2
612: d_eq(1, 10) = -1
613: d_eq(1, 11) = -1
614: d_eq(1, 14) = -2
615: d_eq(1, 16) = -1
616: d_eq(1, 17) = -2
617: d_eq(1, 19) = -1
618: d_eq(1, 20) = -1
619: d_eq(1, 22) = -3
620: d_eq(1, 26) = -1
622: d_eq(2, 2) = -1*an_o_additive
623: d_eq(2, 3) = -2
624: d_eq(2, 4) = -2
625: d_eq(2, 5) = -1
626: d_eq(2, 8) = -1
627: d_eq(2, 9) = -1
628: d_eq(2, 10) = -1
629: d_eq(2, 12) = -1
630: d_eq(2, 13) = -1
631: d_eq(2, 15) = -2
632: d_eq(2, 16) = -2
633: d_eq(2, 20) = -1
634: d_eq(2, 22) = -1
635: d_eq(2, 23) = -1
636: d_eq(2, 24) = -2
637: d_eq(2, 25) = -1
638: d_eq(2, 26) = -1
640: d_eq(6, 6) = -2
641: d_eq(6, 7) = -1
642: d_eq(6, 9) = -1
643: d_eq(6, 12) = -2
644: d_eq(6, 15) = -1
645: d_eq(6, 23) = -1
647: d_eq(4, 1) = -an_c(1)
648: d_eq(4, 2) = -an_c_additive
649: d_eq(4, 4) = -1
650: d_eq(4, 13) = -1
651: d_eq(4, 17) = -2
652: d_eq(4, 18) = -1
653: d_eq(4, 19) = -1
654: d_eq(4, 20) = -1
656: !----------
657: const2 = an_t*an_t
658: const3 = (an_t)*sqrt(an_t)
659: const5 = an_t*const3
661: const_cube = an_t*an_t*an_t
662: const_four = const2*const2
663: const_five = const2*const_cube
664: const_six = const_cube*const_cube
665: pt_cube = pt*pt*pt
666: pt_four = pt_cube*pt
667: pt_five = pt_cube*pt*pt
669: i_cc = an_c(1)
670: i_hh = an_h(1)
671: ai_o2 = i_cc + float(i_hh)/4.0
672: i_co2_h2o = i_cc + i_hh/2
673: i_h2o = i_hh/2
674: ai_pwr1 = 1 + i_cc + float(i_hh)/4.0
675: i_pwr2 = i_cc + i_hh/2
676: idiff = (i_cc + i_h2o) - (ai_o2 + 1)
678: pt1 = pt**(ai_pwr1)
679: an_tot1 = an_t**(ai_pwr1)
680: pt_val1 = (pt/an_t)**(ai_pwr1)
681: an_tot1_d = an_tot1*an_t
683: pt2 = pt**(i_pwr2)
684: an_tot2 = an_t**(i_pwr2)
685: pt_val2 = (pt/an_t)**(i_pwr2)
686: an_tot2_d = an_tot2*an_t
688: d_eq(5, 1) = &
689: & -(an_r(4)**i_cc)*(an_r(5)**i_h2o) &
690: & *((pt/an_t)**idiff)*(-idiff/an_t)
692: do jj = 2, 26
693: d_eq(5, jj) = d_eq(5, 1)
694: end do
696: d_eq(5, 1) = d_eq(5, 1) + k_eq(11)*(an_r(3)**ai_o2)
698: d_eq(5, 3) = d_eq(5, 3) + k_eq(11)*(ai_o2*an_r(3)**(ai_o2 - 1)) &
699: & *an_r(1)
701: d_eq(5, 4) = d_eq(5, 4) - (i_cc*an_r(4)**(i_cc - 1))* &
702: & (an_r(5)**i_h2o)*((pt/an_t)**idiff)
703: d_eq(5, 5) = d_eq(5, 5) &
704: & - (i_h2o*(an_r(5)**(i_h2o - 1))) &
705: & *(an_r(4)**i_cc)*((pt/an_t)**idiff)
707: d_eq(3, 1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
708: do jj = 2, 26
709: d_eq(3, jj) = d_eq(3, 1)
710: end do
712: d_eq(3, 2) = d_eq(3, 2) + k_eq(12)*(an_r(3)**3)
714: d_eq(3, 3) = d_eq(3, 3) + k_eq(12)*(3*an_r(3)**2)*an_r(2)
716: d_eq(3, 4) = d_eq(3, 4) - 2*an_r(4)*(an_r(5)**3)*(pt/an_t)
718: d_eq(3, 5) = d_eq(3, 5) - 3*(an_r(5)**2)*(an_r(4)**2)*(pt/an_t)
719: ! & *(pt_five/const_five)
721: do ii3 = 1, 26
722: d_eq(3, ii3) = 0.0d0
723: end do
724: d_eq(3, 2) = 1.0d0
726: d_eq(7, 1) = pt*an_r(11)*(-1.0)/const2 &
727: & - k_eq(1)*sqrt(pt)*sqrt(an_r(14) + 1d-50)*(-0.5/const3)
729: do jj = 2, 26
730: d_eq(7, jj) = d_eq(7, 1)
731: end do
733: d_eq(7, 11) = d_eq(7, 11) + pt/an_t
734: d_eq(7, 14) = d_eq(7, 14) &
735: & - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14) + 1d-50)*an_t)))
737: d_eq(8, 1) = pt*an_r(8)*(-1.0)/const2 &
738: & - k_eq(2)*sqrt(pt)*sqrt(an_r(3) + 1.0d-50)*(-0.5/const3)
740: do jj = 2, 26
741: d_eq(8, jj) = d_eq(8, 1)
742: end do
744: d_eq(8, 3) = d_eq(8, 3) &
745: & - k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3) + 1.0d-50)*an_t)))
746: d_eq(8, 8) = d_eq(8, 8) + pt/an_t
748: d_eq(9, 1) = pt*an_r(7)*(-1.0)/const2 &
749: & - k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3)
751: do jj = 2, 26
752: d_eq(9, jj) = d_eq(9, 1)
753: end do
755: d_eq(9, 7) = d_eq(9, 7) + pt/an_t
756: d_eq(9, 6) = d_eq(9, 6) &
757: & - k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))
759: d_eq(10, 1) = pt*an_r(10)*(-1.0)/const2 &
760: & - k_eq(4)*(pt)*sqrt((an_r(3) + 1.0d-50) &
761: & *an_r(14))*(-1.0/const2)
762: do jj = 2, 26
763: d_eq(10, jj) = d_eq(10, 1)
764: end do
766: d_eq(10, 3) = d_eq(10, 3) &
767: & - k_eq(4)*(pt)*sqrt(an_r(14)) &
768: & *(0.5/(sqrt(an_r(3) + 1.0d-50)*an_t))
769: d_eq(10, 10) = d_eq(10, 10) + pt/an_t
770: d_eq(10, 14) = d_eq(10, 14) &
771: & - k_eq(4)*(pt)*sqrt(an_r(3) + 1.0d-50) &
772: & *(0.5/(sqrt(an_r(14) + 1.0d-50)*an_t))
774: d_eq(11, 1) = pt*an_r(9)*(-1.0)/const2 &
775: & - k_eq(5)*(pt)*sqrt((an_r(3) + 1.0d-50)*an_r(6)) &
776: & *(-1.0/const2)
778: do jj = 2, 26
779: d_eq(11, jj) = d_eq(11, 1)
780: end do
782: d_eq(11, 3) = d_eq(11, 3) &
783: & - k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ &
784: & (sqrt(an_r(3) + 1.0d-50)*an_t))
785: d_eq(11, 6) = d_eq(11, 6) &
786: & - k_eq(5)*(pt)*sqrt(an_r(3) + 1.0d-50) &
787: & *(0.5/(sqrt(an_r(6))*an_t))
788: d_eq(11, 9) = d_eq(11, 9) + pt/an_t
790: d_eq(12, 1) = pt*an_r(5)*(-1.0)/const2 &
791: & - k_eq(6)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
792: & *(an_r(14))*(-1.5/const5)
794: do jj = 2, 26
795: d_eq(12, jj) = d_eq(12, 1)
796: end do
798: d_eq(12, 3) = d_eq(12, 3) &
799: & - k_eq(6)*(pt**1.5)*((an_r(14) + 1.0d-50)/const3) &
800: & *(0.5/sqrt(an_r(3) + 1.0d-50))
802: d_eq(12, 5) = d_eq(12, 5) + pt/an_t
803: d_eq(12, 14) = d_eq(12, 14) &
804: & - k_eq(6)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
806: d_eq(13, 1) = pt*an_r(4)*(-1.0)/const2 &
807: & - k_eq(7)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
808: & *(an_r(13))*(-1.5/const5)
810: do jj = 2, 26
811: d_eq(13, jj) = d_eq(13, 1)
812: end do
814: d_eq(13, 3) = d_eq(13, 3) &
815: & - k_eq(7)*(pt**1.5)*(an_r(13)/const3) &
816: & *(0.5/sqrt(an_r(3) + 1.0d-50))
818: d_eq(13, 4) = d_eq(13, 4) + pt/an_t
819: d_eq(13, 13) = d_eq(13, 13) &
820: & - k_eq(7)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
822: d_eq(14, 1) = pt*an_r(15)*(-1.0)/const2 &
823: & - k_eq(8)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
824: & *(an_r(9))*(-1.5/const5)
826: do jj = 2, 26
827: d_eq(14, jj) = d_eq(14, 1)
828: end do
830: d_eq(14, 3) = d_eq(14, 3) &
831: & - k_eq(8)*(pt**1.5)*(an_r(9)/const3) &
832: & *(0.5/sqrt(an_r(3) + 1.0d-50))
833: d_eq(14, 9) = d_eq(14, 9) &
834: & - k_eq(8)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
835: d_eq(14, 15) = d_eq(14, 15) + pt/an_t
837: d_eq(15, 1) = pt*an_r(16)*(-1.0)/const2 &
838: & - k_eq(9)*(pt**1.5)*sqrt(an_r(14) + 1.0d-50) &
839: & *(an_r(3))*(-1.5/const5)
841: do jj = 2, 26
842: d_eq(15, jj) = d_eq(15, 1)
843: end do
845: d_eq(15, 3) = d_eq(15, 3) &
846: & - k_eq(9)*(pt**1.5)*(sqrt(an_r(14) + 1.0d-50)/const3)
847: d_eq(15, 14) = d_eq(15, 14) &
848: & - k_eq(9)*(pt**1.5)*(an_r(3)/const3) &
849: & *(0.5/sqrt(an_r(14) + 1.0d-50))
850: d_eq(15, 16) = d_eq(15, 16) + pt/an_t
852: d_eq(16, 1) = pt*an_r(12)*(-1.0)/const2 &
853: & - k_eq(10)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
854: & *(an_r(6))*(-1.5/const5)
856: do jj = 2, 26
857: d_eq(16, jj) = d_eq(16, 1)
858: end do
860: d_eq(16, 3) = d_eq(16, 3) &
861: & - k_eq(10)*(pt**1.5)*(an_r(6)/const3) &
862: & *(0.5/sqrt(an_r(3) + 1.0d-50))
864: d_eq(16, 6) = d_eq(16, 6) &
865: & - k_eq(10)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
866: d_eq(16, 12) = d_eq(16, 12) + pt/an_t
868: const_cube = an_t*an_t*an_t
869: const_four = const2*const2
871: d_eq(17, 1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
872: & - k_eq(15)*an_r(17)*pt*(-1/const2)
873: do jj = 2, 26
874: d_eq(17, jj) = d_eq(17, 1)
875: end do
876: d_eq(17, 14) = d_eq(17, 14) + an_r(18)*an_r(18)*(pt**3)/const_cube
877: d_eq(17, 17) = d_eq(17, 17) - k_eq(15)*pt/an_t
878: d_eq(17, 18) = d_eq(17, 18) + 2*an_r(18)*an_r(14) &
879: & *(pt**3)/const_cube
881: d_eq(18, 1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube) &
882: & - k_eq(16)*an_r(3)*an_r(18)*an_r(18) &
883: & *(pt*pt*pt)*(-3/const_four)
884: do jj = 2, 26
885: d_eq(18, jj) = d_eq(18, 1)
886: end do
887: d_eq(18, 3) = d_eq(18, 3) &
888: & - k_eq(16)*an_r(18)*an_r(18)*pt*pt*pt/const_cube
889: d_eq(18, 13) = d_eq(18, 13) &
890: & + 2*an_r(13)*pt*pt/const2
891: d_eq(18, 18) = d_eq(18, 18) - k_eq(16)*an_r(3) &
892: & *2*an_r(18)*pt*pt*pt/const_cube
894: !====for eq 19
896: d_eq(19, 1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube) &
897: & - k_eq(17)*an_r(13)*an_r(10)*pt*pt*(-2/const_cube)
898: do jj = 2, 26
899: d_eq(19, jj) = d_eq(19, 1)
900: end do
901: d_eq(19, 13) = d_eq(19, 13) &
902: & - k_eq(17)*an_r(10)*pt*pt/const2
903: d_eq(19, 10) = d_eq(19, 10) &
904: & - k_eq(17)*an_r(13)*pt*pt/const2
905: d_eq(19, 3) = d_eq(19, 3) + an_r(19)*pt*pt/const2
906: d_eq(19, 19) = d_eq(19, 19) + an_r(3)*pt*pt/const2
907: !====for eq 20
909: d_eq(20, 1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube) &
910: & - k_eq(18)*an_r(19)*an_r(8)*pt*pt*(-2/const_cube)
911: do jj = 2, 26
912: d_eq(20, jj) = d_eq(20, 1)
913: end do
914: d_eq(20, 8) = d_eq(20, 8) &
915: & - k_eq(18)*an_r(19)*pt*pt/const2
916: d_eq(20, 19) = d_eq(20, 19) &
917: & - k_eq(18)*an_r(8)*pt*pt/const2
918: d_eq(20, 20) = d_eq(20, 20) + an_r(21)*pt*pt/const2
919: d_eq(20, 21) = d_eq(20, 21) + an_r(20)*pt*pt/const2
921: !========
922: !====for eq 21
924: d_eq(21, 1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube) &
925: & - k_eq(19)*an_r(7)*an_r(8)*pt*pt*(-2/const_cube)
926: do jj = 2, 26
927: d_eq(21, jj) = d_eq(21, 1)
928: end do
929: d_eq(21, 7) = d_eq(21, 7) &
930: & - k_eq(19)*an_r(8)*pt*pt/const2
931: d_eq(21, 8) = d_eq(21, 8) &
932: & - k_eq(19)*an_r(7)*pt*pt/const2
933: d_eq(21, 21) = d_eq(21, 21) + an_r(23)*pt*pt/const2
934: d_eq(21, 23) = d_eq(21, 23) + an_r(21)*pt*pt/const2
936: !========
937: ! for 22
938: d_eq(22, 1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube) &
939: & - k_eq(20)*an_r(21)*an_r(22)*pt*pt*(-2/const_cube)
940: do jj = 2, 26
941: d_eq(22, jj) = d_eq(22, 1)
942: end do
943: d_eq(22, 21) = d_eq(22, 21) &
944: & - k_eq(20)*an_r(22)*pt*pt/const2
945: d_eq(22, 22) = d_eq(22, 22) &
946: & - k_eq(20)*an_r(21)*pt*pt/const2
947: d_eq(22, 11) = d_eq(22, 11) + an_r(5)*pt*pt/(const2)
948: d_eq(22, 5) = d_eq(22, 5) + an_r(11)*pt*pt/(const2)
950: !========
951: ! for 23
953: d_eq(23, 1) = an_r(24)*(pt)*(-1/const2) &
954: & - k_eq(21)*an_r(21)*an_r(3)*pt*pt*(-2/const_cube)
955: do jj = 2, 26
956: d_eq(23, jj) = d_eq(23, 1)
957: end do
958: d_eq(23, 3) = d_eq(23, 3) &
959: & - k_eq(21)*an_r(21)*pt*pt/const2
960: d_eq(23, 21) = d_eq(23, 21) &
961: & - k_eq(21)*an_r(3)*pt*pt/const2
962: d_eq(23, 24) = d_eq(23, 24) + pt/(an_t)
964: !========
965: ! for 24
966: d_eq(24, 1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube) &
967: & - k_eq(22)*an_r(24)*an_r(8)*pt*pt*(-2/const_cube)
968: do jj = 2, 26
969: d_eq(24, jj) = d_eq(24, 1)
970: end do
971: d_eq(24, 8) = d_eq(24, 8) &
972: & - k_eq(22)*an_r(24)*pt*pt/const2
973: d_eq(24, 24) = d_eq(24, 24) &
974: & - k_eq(22)*an_r(8)*pt*pt/const2
975: d_eq(24, 3) = d_eq(24, 3) + an_r(25)*pt*pt/const2
976: d_eq(24, 25) = d_eq(24, 25) + an_r(3)*pt*pt/const2
978: !========
979: !for 25
981: d_eq(25, 1) = an_r(26)*(pt)*(-1/const2) &
982: & - k_eq(23)*an_r(21)*an_r(10)*pt*pt*(-2/const_cube)
983: do jj = 2, 26
984: d_eq(25, jj) = d_eq(25, 1)
985: end do
986: d_eq(25, 10) = d_eq(25, 10) &
987: & - k_eq(23)*an_r(21)*pt*pt/const2
988: d_eq(25, 21) = d_eq(25, 21) &
989: & - k_eq(23)*an_r(10)*pt*pt/const2
990: d_eq(25, 26) = d_eq(25, 26) + pt/(an_t)
992: !============
993: ! for 26
994: d_eq(26, 20) = -1
995: d_eq(26, 22) = -1
996: d_eq(26, 23) = -1
997: d_eq(26, 21) = 1
998: d_eq(26, 24) = 1
999: d_eq(26, 25) = 1
1000: d_eq(26, 26) = 1
1002: do j = 1, 26
1003: do i = 1, 26
1004: write (44, *) i, j, d_eq(i, j)
1005: end do
1006: end do
1008: end
1010: subroutine ShashiPostCheck(ls, X, Y, W, c_Y, c_W, dummy)
1011: #include <petsc/finclude/petscsnes.h>
1012: use petscsnes
1013: implicit none
1014: SNESLineSearch ls
1015: PetscErrorCode ierr
1016: Vec X, Y, W
1017: PetscObject dummy
1018: PetscBool c_Y, c_W
1019: PetscScalar, pointer :: xx(:)
1020: PetscInt i
1021: PetscCall(VecGetArray(W, xx, ierr))
1022: do i = 1, 26
1023: if (xx(i) < 0.0) then
1024: xx(i) = 0.0
1025: c_W = PETSC_TRUE
1026: end if
1027: if (xx(i) > 3.0) then
1028: xx(i) = 3.0
1029: end if
1030: end do
1031: PetscCall(VecRestoreArray(W, xx, ierr))
1032: end