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: ! Macro definitions
34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: !
36: ! Macros to make clearer the process of setting values in vectors and
37: ! getting values from vectors. These vectors are used in the routines
38: ! FormFunction() and FormJacobian().
39: ! - The element lx_a(ib) is element ib in the vector x
40: !
41: #define lx_a(ib) lx_v(lx_i + (ib))
42: #define lf_a(ib) lf_v(lf_i + (ib))
43: !
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Beginning of program
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: PetscCallA(PetscInitialize(ierr))
49: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
50: if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,1,'requires one process'); endif
52: big = 2.88
53: big = PETSC_INFINITY
54: zero = 0.0
55: i2 = 26
56: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
57: ! Create nonlinear solver context
58: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
60: PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: ! Create matrix and vector data structures; set corresponding routines
64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: ! Create vectors for solution and nonlinear function
68: PetscCallA(VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr))
69: PetscCallA(VecDuplicate(x,r,ierr))
71: ! Create Jacobian matrix data structure
73: PetscCallA(MatCreateDense(PETSC_COMM_SELF,26,26,26,26,PETSC_NULL_SCALAR,J,ierr))
75: ! Set function evaluation routine and vector
77: PetscCallA(SNESSetFunction(snes,r,FormFunction,0,ierr))
79: ! Set Jacobian matrix data structure and Jacobian evaluation routine
81: PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,0,ierr))
83: PetscCallA(VecDuplicate(x,lb,ierr))
84: PetscCallA(VecDuplicate(x,ub,ierr))
85: PetscCallA(VecSet(lb,zero,ierr))
86: ! PetscCallA(VecGetArrayF90(lb,xx,ierr))
87: ! PetscCallA(ShashiLowerBound(xx))
88: ! PetscCallA(VecRestoreArrayF90(lb,xx,ierr))
89: PetscCallA(VecSet(ub,big,ierr))
90: ! PetscCallA(SNESVISetVariableBounds(snes,lb,ub,ierr))
92: PetscCallA(SNESGetLineSearch(snes,ls,ierr))
93: PetscCallA(SNESLineSearchSetPostCheck(ls,ShashiPostCheck,0,ierr))
94: PetscCallA(SNESSetType(snes,SNESVINEWTONRSLS,ierr))
96: PetscCallA(SNESSetFromOptions(snes,ierr))
98: ! set initial guess
100: PetscCallA(VecGetArrayF90(x,xx,ierr))
101: PetscCallA(ShashiInitialGuess(xx))
102: PetscCallA(VecRestoreArrayF90(x,xx,ierr))
104: PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: ! Free work space. All PETSc objects should be destroyed when they
108: ! are no longer needed.
109: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: PetscCallA(VecDestroy(lb,ierr))
111: PetscCallA(VecDestroy(ub,ierr))
112: PetscCallA(VecDestroy(x,ierr))
113: PetscCallA(VecDestroy(r,ierr))
114: PetscCallA(MatDestroy(J,ierr))
115: PetscCallA(SNESDestroy(snes,ierr))
116: PetscCallA(PetscFinalize(ierr))
117: end
118: !
119: ! ------------------------------------------------------------------------
120: !
121: ! FormFunction - Evaluates nonlinear function, F(x).
122: !
123: ! Input Parameters:
124: ! snes - the SNES context
125: ! x - input vector
126: ! dummy - optional user-defined context (not used here)
127: !
128: ! Output Parameter:
129: ! f - function vector
130: !
131: subroutine FormFunction(snes,x,f,dummy,ierr)
132: #include <petsc/finclude/petscsnes.h>
133: use petscsnes
134: implicit none
135: SNES snes
136: Vec x,f
137: PetscErrorCode ierr
138: integer dummy(*)
140: ! Declarations for use with local arrays
142: PetscScalar lx_v(2),lf_v(2)
143: PetscOffset lx_i,lf_i
145: ! Get pointers to vector data.
146: ! - For default PETSc vectors, VecGetArray() returns a pointer to
147: ! the data array. Otherwise, the routine is implementation dependent.
148: ! - You MUST call VecRestoreArray() when you no longer need access to
149: ! the array.
150: ! - Note that the Fortran interface to VecGetArray() differs from the
151: ! C version. See the Fortran chapter of the users manual for details.
153: PetscCall(VecGetArrayRead(x,lx_v,lx_i,ierr))
154: PetscCall(VecGetArray(f,lf_v,lf_i,ierr))
155: PetscCall(ShashiFormFunction(lx_a(1),lf_a(1)))
156: PetscCall(VecRestoreArrayRead(x,lx_v,lx_i,ierr))
157: PetscCall(VecRestoreArray(f,lf_v,lf_i,ierr))
159: return
160: end
162: ! ---------------------------------------------------------------------
163: !
164: ! FormJacobian - Evaluates Jacobian matrix.
165: !
166: ! Input Parameters:
167: ! snes - the SNES context
168: ! x - input vector
169: ! dummy - optional user-defined context (not used here)
170: !
171: ! Output Parameters:
172: ! A - Jacobian matrix
173: ! B - optionally different preconditioning matrix
174: ! flag - flag indicating matrix structure
175: !
176: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
177: #include <petsc/finclude/petscsnes.h>
178: use petscsnes
179: implicit none
180: SNES snes
181: Vec X
182: Mat jac,B
183: PetscErrorCode ierr
184: integer dummy(*)
186: ! Declarations for use with local arrays
188: PetscScalar lx_v(1),lf_v(1)
189: PetscOffset lx_i,lf_i
191: ! Get pointer to vector data
193: PetscCall(VecGetArrayRead(x,lx_v,lx_i,ierr))
194: PetscCall(MatDenseGetArray(B,lf_v,lf_i,ierr))
195: PetscCall(ShashiFormJacobian(lx_a(1),lf_a(1)))
196: PetscCall(MatDenseRestoreArray(B,lf_v,lf_i,ierr))
197: PetscCall(VecRestoreArrayRead(x,lx_v,lx_i,ierr))
199: ! Assemble matrix
201: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
202: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
204: return
205: end
207: subroutine ShashiLowerBound(an_r)
208: ! implicit PetscScalar (a-h,o-z)
209: implicit none
210: PetscScalar an_r(26)
211: PetscInt i
213: do i=2,26
214: an_r(i) = 1000.0/6.023D+23
215: enddo
216: return
217: end
219: subroutine ShashiInitialGuess(an_r)
220: ! implicit PetscScalar (a-h,o-z)
221: implicit none
222: PetscScalar an_c_additive
223: PetscScalar an_h_additive
224: PetscScalar an_o_additive
225: PetscScalar atom_c_init
226: PetscScalar atom_h_init
227: PetscScalar atom_n_init
228: PetscScalar atom_o_init
229: PetscScalar h_init
230: PetscScalar p_init
231: PetscInt nfuel
232: PetscScalar temp,pt
233: PetscScalar an_r(26)
234: PetscInt an_h(1),an_c(1)
236: pt = 0.1
237: atom_c_init =6.7408177364816552D-022
238: atom_h_init = 2.0
239: atom_o_init = 1.0
240: atom_n_init = 3.76
241: nfuel = 1
242: an_c(1) = 1
243: an_h(1) = 4
244: an_c_additive = 2
245: an_h_additive = 6
246: an_o_additive = 1
247: h_init = 128799.7267952987
248: p_init = 0.1
249: temp = 1500
251: an_r( 1) = 1.66000D-24
252: an_r( 2) = 1.66030D-22
253: an_r( 3) = 5.00000D-01
254: an_r( 4) = 1.66030D-22
255: an_r( 5) = 1.66030D-22
256: an_r( 6) = 1.88000D+00
257: an_r( 7) = 1.66030D-22
258: an_r( 8) = 1.66030D-22
259: an_r( 9) = 1.66030D-22
260: an_r(10) = 1.66030D-22
261: an_r(11) = 1.66030D-22
262: an_r(12) = 1.66030D-22
263: an_r(13) = 1.66030D-22
264: an_r(14) = 1.00000D+00
265: an_r(15) = 1.66030D-22
266: an_r(16) = 1.66030D-22
267: an_r(17) = 1.66000D-24
268: an_r(18) = 1.66030D-24
269: an_r(19) = 1.66030D-24
270: an_r(20) = 1.66030D-24
271: an_r(21) = 1.66030D-24
272: an_r(22) = 1.66030D-24
273: an_r(23) = 1.66030D-24
274: an_r(24) = 1.66030D-24
275: an_r(25) = 1.66030D-24
276: an_r(26) = 1.66030D-24
278: an_r = 0
279: an_r( 3) = .5
280: an_r( 6) = 1.88000
281: an_r(14) = 1.
283: #if defined(solution)
284: an_r( 1) = 3.802208D-33
285: an_r( 2) = 1.298287D-29
286: an_r( 3) = 2.533067D-04
287: an_r( 4) = 6.865078D-22
288: an_r( 5) = 9.993125D-01
289: an_r( 6) = 1.879964D+00
290: an_r( 7) = 4.449489D-13
291: an_r( 8) = 3.428687D-07
292: an_r( 9) = 7.105138D-05
293: an_r(10) = 1.094368D-04
294: an_r(11) = 2.362305D-06
295: an_r(12) = 1.107145D-09
296: an_r(13) = 1.276162D-24
297: an_r(14) = 6.315538D-04
298: an_r(15) = 2.356540D-09
299: an_r(16) = 2.048248D-09
300: an_r(17) = 1.966187D-22
301: an_r(18) = 7.856497D-29
302: an_r(19) = 1.987840D-36
303: an_r(20) = 8.182441D-22
304: an_r(21) = 2.684880D-16
305: an_r(22) = 2.680473D-16
306: an_r(23) = 6.594967D-18
307: an_r(24) = 2.509714D-21
308: an_r(25) = 3.096459D-21
309: an_r(26) = 6.149551D-18
310: #endif
312: return
313: end
315: subroutine ShashiFormFunction(an_r,f_eq)
316: ! implicit PetscScalar (a-h,o-z)
317: implicit none
318: PetscScalar an_c_additive
319: PetscScalar an_h_additive
320: PetscScalar an_o_additive
321: PetscScalar atom_c_init
322: PetscScalar atom_h_init
323: PetscScalar atom_n_init
324: PetscScalar atom_o_init
325: PetscScalar h_init
326: PetscScalar p_init
327: PetscInt nfuel
328: PetscScalar temp,pt
329: PetscScalar an_r(26),k_eq(23),f_eq(26)
330: PetscScalar H_molar(26)
331: PetscInt an_h(1),an_c(1)
332: PetscScalar part_p(26),idiff
333: PetscInt i_cc,i_hh,i_h2o
334: PetscScalar an_t,sum_h
335: PetscScalar a_io2
336: PetscInt i,ip
337: pt = 0.1
338: atom_c_init =6.7408177364816552D-022
339: atom_h_init = 2.0
340: atom_o_init = 1.0
341: atom_n_init = 3.76
342: nfuel = 1
343: an_c(1) = 1
344: an_h(1) = 4
345: an_c_additive = 2
346: an_h_additive = 6
347: an_o_additive = 1
348: h_init = 128799.7267952987
349: p_init = 0.1
350: temp = 1500
352: k_eq( 1) = 1.75149D-05
353: k_eq( 2) = 4.01405D-06
354: k_eq( 3) = 6.04663D-14
355: k_eq( 4) = 2.73612D-01
356: k_eq( 5) = 3.25592D-03
357: k_eq( 6) = 5.33568D+05
358: k_eq( 7) = 2.07479D+05
359: k_eq( 8) = 1.11841D-02
360: k_eq( 9) = 1.72684D-03
361: k_eq(10) = 1.98588D-07
362: k_eq(11) = 7.23600D+27
363: k_eq(12) = 5.73926D+49
364: k_eq(13) = 1.00000D+00
365: k_eq(14) = 1.64493D+16
366: k_eq(15) = 2.73837D-29
367: k_eq(16) = 3.27419D+50
368: k_eq(17) = 1.72447D-23
369: k_eq(18) = 4.24657D-06
370: k_eq(19) = 1.16065D-14
371: k_eq(20) = 3.28020D+25
372: k_eq(21) = 1.06291D+00
373: k_eq(22) = 9.11507D+02
374: k_eq(23) = 6.02837D+03
376: H_molar( 1) = 3.26044D+03
377: H_molar( 2) = -8.00407D+04
378: H_molar( 3) = 4.05873D+04
379: H_molar( 4) = -3.31849D+05
380: H_molar( 5) = -1.93654D+05
381: H_molar( 6) = 3.84035D+04
382: H_molar( 7) = 4.97589D+05
383: H_molar( 8) = 2.74483D+05
384: H_molar( 9) = 1.30022D+05
385: H_molar(10) = 7.58429D+04
386: H_molar(11) = 2.42948D+05
387: H_molar(12) = 1.44588D+05
388: H_molar(13) = -7.16891D+04
389: H_molar(14) = 3.63075D+04
390: H_molar(15) = 9.23880D+04
391: H_molar(16) = 6.50477D+04
392: H_molar(17) = 3.04310D+05
393: H_molar(18) = 7.41707D+05
394: H_molar(19) = 6.32767D+05
395: H_molar(20) = 8.90624D+05
396: H_molar(21) = 2.49805D+04
397: H_molar(22) = 6.43473D+05
398: H_molar(23) = 1.02861D+06
399: H_molar(24) = -6.07503D+03
400: H_molar(25) = 1.27020D+05
401: H_molar(26) = -1.07011D+05
402: !=============
403: an_t = 0
404: sum_h = 0
406: do i = 1,26
407: an_t = an_t + an_r(i)
408: enddo
410: f_eq(1) = atom_h_init &
411: & - (an_h(1)*an_r(1) + an_h_additive*an_r(2) &
412: & + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) &
413: & + an_r(16)+ 2*an_r(17) + an_r(19) &
414: & +an_r(20) + 3*an_r(22)+an_r(26))
416: f_eq(2) = atom_o_init &
417: & - (an_o_additive*an_r(2) + 2*an_r(3) &
418: & + 2*an_r(4) + an_r(5) &
419: & + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
420: & + 2*an_r(15) + 2*an_r(16)+ an_r(20) + an_r(22) &
421: & + an_r(23) + 2*an_r(24) + 1*an_r(25)+an_r(26))
423: f_eq(3) = an_r(2)-1.0d-150
425: f_eq(4) = atom_c_init &
426: & - (an_c(1)*an_r(1) + an_c_additive * an_r(2) &
427: & + an_r(4) + an_r(13)+ 2*an_r(17) + an_r(18) &
428: & + an_r(19) + an_r(20))
430: do ip = 1,26
431: part_p(ip) = (an_r(ip)/an_t) * pt
432: enddo
434: i_cc = an_c(1)
435: i_hh = an_h(1)
436: a_io2 = i_cc + i_hh/4.0
437: i_h2o = i_hh/2
438: idiff = (i_cc + i_h2o) - (a_io2 + 1)
440: f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2 &
441: & - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
442: ! write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
443: ! stop
444: f_eq(6) = atom_n_init &
445: & - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) &
446: & + an_r(15) &
447: & + an_r(23))
449: f_eq( 7) = part_p(11) &
450: & - (k_eq( 1) * sqrt(part_p(14)+1d-23))
451: f_eq( 8) = part_p( 8) &
452: & - (k_eq( 2) * sqrt(part_p( 3)+1d-23))
454: f_eq( 9) = part_p( 7) &
455: & - (k_eq( 3) * sqrt(part_p( 6)+1d-23))
457: f_eq(10) = part_p(10) &
458: & - (k_eq( 4) * sqrt(part_p( 3)+1d-23)) &
459: & * sqrt(part_p(14))
461: f_eq(11) = part_p( 9) &
462: & - (k_eq( 5) * sqrt(part_p( 3)+1d-23)) &
463: & * sqrt(part_p( 6)+1d-23)
464: f_eq(12) = part_p( 5) &
465: & - (k_eq( 6) * sqrt(part_p( 3)+1d-23)) &
466: & * (part_p(14))
468: f_eq(13) = part_p( 4) &
469: & - (k_eq( 7) * sqrt(part_p(3)+1.0d-23)) &
470: & * (part_p(13))
472: f_eq(14) = part_p(15) &
473: & - (k_eq( 8) * sqrt(part_p(3)+1.0d-50)) &
474: & * (part_p( 9))
475: f_eq(15) = part_p(16) &
476: & - (k_eq( 9) * part_p( 3)) &
477: & * sqrt(part_p(14)+1d-23)
479: f_eq(16) = part_p(12) &
480: & - (k_eq(10) * sqrt(part_p( 3)+1d-23)) &
481: & * (part_p( 6))
483: f_eq(17) = part_p(14)*part_p(18)**2 &
484: & - (k_eq(15) * part_p(17))
486: f_eq(18) = (part_p(13)**2) &
487: & - (k_eq(16) * part_p(3)*part_p(18)**2)
488: print*,f_eq(18),part_p(3),part_p(18),part_p(13),k_eq(16)
490: f_eq(19) = part_p(19)*part_p(3) - k_eq(17)*part_p(13)*part_p(10)
492: f_eq(20) = part_p(21)*part_p(20) - k_eq(18)*part_p(19)*part_p(8)
494: f_eq(21) = part_p(21)*part_p(23) - k_eq(19)*part_p(7)*part_p(8)
496: f_eq(22) = part_p(5)*part_p(11) - k_eq(20)*part_p(21)*part_p(22)
498: f_eq(23) = part_p(24) - k_eq(21)*part_p(21)*part_p(3)
500: f_eq(24) = part_p(3)*part_p(25) - k_eq(22)*part_p(24)*part_p(8)
502: f_eq(25) = part_p(26) - k_eq(23)*part_p(21)*part_p(10)
504: f_eq(26) = -(an_r(20) + an_r(22) + an_r(23)) &
505: & +(an_r(21) + an_r(24) + an_r(25) + an_r(26))
507: do i = 1,26
508: write(44,*)i,f_eq(i)
509: enddo
511: return
512: end
514: subroutine ShashiFormJacobian(an_r,d_eq)
515: ! implicit PetscScalar (a-h,o-z)
516: implicit none
517: PetscScalar an_c_additive
518: PetscScalar an_h_additive
519: PetscScalar an_o_additive
520: PetscScalar atom_c_init
521: PetscScalar atom_h_init
522: PetscScalar atom_n_init
523: PetscScalar atom_o_init
524: PetscScalar h_init
525: PetscScalar p_init
526: PetscInt nfuel
527: PetscScalar temp,pt
528: PetscScalar an_t,ai_o2
529: PetscScalar an_tot1_d,an_tot1
530: PetscScalar an_tot2_d,an_tot2
531: PetscScalar const5,const3,const2
532: PetscScalar const_cube
533: PetscScalar const_five
534: PetscScalar const_four
535: PetscScalar const_six
536: PetscInt jj,jb,ii3,id,ib,i
537: PetscScalar pt2,pt1
538: PetscScalar an_r(26),k_eq(23)
539: PetscScalar d_eq(26,26),H_molar(26)
540: PetscInt an_h(1),an_c(1)
541: PetscScalar ai_pwr1,idiff
542: PetscInt i_cc,i_hh
543: PetscInt i_h2o,i_pwr2,i_co2_h2o
544: PetscScalar pt_cube,pt_five
545: PetscScalar pt_four
546: PetscScalar pt_val1,pt_val2
547: PetscInt j
549: pt = 0.1
550: atom_c_init =6.7408177364816552D-022
551: atom_h_init = 2.0
552: atom_o_init = 1.0
553: atom_n_init = 3.76
554: nfuel = 1
555: an_c(1) = 1
556: an_h(1) = 4
557: an_c_additive = 2
558: an_h_additive = 6
559: an_o_additive = 1
560: h_init = 128799.7267952987
561: p_init = 0.1
562: temp = 1500
564: k_eq( 1) = 1.75149D-05
565: k_eq( 2) = 4.01405D-06
566: k_eq( 3) = 6.04663D-14
567: k_eq( 4) = 2.73612D-01
568: k_eq( 5) = 3.25592D-03
569: k_eq( 6) = 5.33568D+05
570: k_eq( 7) = 2.07479D+05
571: k_eq( 8) = 1.11841D-02
572: k_eq( 9) = 1.72684D-03
573: k_eq(10) = 1.98588D-07
574: k_eq(11) = 7.23600D+27
575: k_eq(12) = 5.73926D+49
576: k_eq(13) = 1.00000D+00
577: k_eq(14) = 1.64493D+16
578: k_eq(15) = 2.73837D-29
579: k_eq(16) = 3.27419D+50
580: k_eq(17) = 1.72447D-23
581: k_eq(18) = 4.24657D-06
582: k_eq(19) = 1.16065D-14
583: k_eq(20) = 3.28020D+25
584: k_eq(21) = 1.06291D+00
585: k_eq(22) = 9.11507D+02
586: k_eq(23) = 6.02837D+03
588: H_molar( 1) = 3.26044D+03
589: H_molar( 2) = -8.00407D+04
590: H_molar( 3) = 4.05873D+04
591: H_molar( 4) = -3.31849D+05
592: H_molar( 5) = -1.93654D+05
593: H_molar( 6) = 3.84035D+04
594: H_molar( 7) = 4.97589D+05
595: H_molar( 8) = 2.74483D+05
596: H_molar( 9) = 1.30022D+05
597: H_molar(10) = 7.58429D+04
598: H_molar(11) = 2.42948D+05
599: H_molar(12) = 1.44588D+05
600: H_molar(13) = -7.16891D+04
601: H_molar(14) = 3.63075D+04
602: H_molar(15) = 9.23880D+04
603: H_molar(16) = 6.50477D+04
604: H_molar(17) = 3.04310D+05
605: H_molar(18) = 7.41707D+05
606: H_molar(19) = 6.32767D+05
607: H_molar(20) = 8.90624D+05
608: H_molar(21) = 2.49805D+04
609: H_molar(22) = 6.43473D+05
610: H_molar(23) = 1.02861D+06
611: H_molar(24) = -6.07503D+03
612: H_molar(25) = 1.27020D+05
613: H_molar(26) = -1.07011D+05
615: !
616: !=======
617: do jb = 1,26
618: do ib = 1,26
619: d_eq(ib,jb) = 0.0d0
620: end do
621: end do
623: an_t = 0.0
624: do id = 1,26
625: an_t = an_t + an_r(id)
626: enddo
628: !====
629: !====
630: d_eq(1,1) = -an_h(1)
631: d_eq(1,2) = -an_h_additive
632: d_eq(1,5) = -2
633: d_eq(1,10) = -1
634: d_eq(1,11) = -1
635: d_eq(1,14) = -2
636: d_eq(1,16) = -1
637: d_eq(1,17) = -2
638: d_eq(1,19) = -1
639: d_eq(1,20) = -1
640: d_eq(1,22) = -3
641: d_eq(1,26) = -1
643: d_eq(2,2) = -1*an_o_additive
644: d_eq(2,3) = -2
645: d_eq(2,4) = -2
646: d_eq(2,5) = -1
647: d_eq(2,8) = -1
648: d_eq(2,9) = -1
649: d_eq(2,10) = -1
650: d_eq(2,12) = -1
651: d_eq(2,13) = -1
652: d_eq(2,15) = -2
653: d_eq(2,16) = -2
654: d_eq(2,20) = -1
655: d_eq(2,22) = -1
656: d_eq(2,23) = -1
657: d_eq(2,24) = -2
658: d_eq(2,25) = -1
659: d_eq(2,26) = -1
661: d_eq(6,6) = -2
662: d_eq(6,7) = -1
663: d_eq(6,9) = -1
664: d_eq(6,12) = -2
665: d_eq(6,15) = -1
666: d_eq(6,23) = -1
668: d_eq(4,1) = -an_c(1)
669: d_eq(4,2) = -an_c_additive
670: d_eq(4,4) = -1
671: d_eq(4,13) = -1
672: d_eq(4,17) = -2
673: d_eq(4,18) = -1
674: d_eq(4,19) = -1
675: d_eq(4,20) = -1
677: !----------
678: const2 = an_t*an_t
679: const3 = (an_t)*sqrt(an_t)
680: const5 = an_t*const3
682: const_cube = an_t*an_t*an_t
683: const_four = const2*const2
684: const_five = const2*const_cube
685: const_six = const_cube*const_cube
686: pt_cube = pt*pt*pt
687: pt_four = pt_cube*pt
688: pt_five = pt_cube*pt*pt
690: i_cc = an_c(1)
691: i_hh = an_h(1)
692: ai_o2 = i_cc + float(i_hh)/4.0
693: i_co2_h2o = i_cc + i_hh/2
694: i_h2o = i_hh/2
695: ai_pwr1 = 1 + i_cc + float(i_hh)/4.0
696: i_pwr2 = i_cc + i_hh/2
697: idiff = (i_cc + i_h2o) - (ai_o2 + 1)
699: pt1 = pt**(ai_pwr1)
700: an_tot1 = an_t**(ai_pwr1)
701: pt_val1 = (pt/an_t)**(ai_pwr1)
702: an_tot1_d = an_tot1*an_t
704: pt2 = pt**(i_pwr2)
705: an_tot2 = an_t**(i_pwr2)
706: pt_val2 = (pt/an_t)**(i_pwr2)
707: an_tot2_d = an_tot2*an_t
709: d_eq(5,1) = &
710: & -(an_r(4)**i_cc)*(an_r(5)**i_h2o) &
711: & *((pt/an_t)**idiff) *(-idiff/an_t)
713: do jj = 2,26
714: d_eq(5,jj) = d_eq(5,1)
715: enddo
717: d_eq(5,1) = d_eq(5,1) + k_eq(11)* (an_r(3) **ai_o2)
719: d_eq(5,3) = d_eq(5,3) + k_eq(11)* (ai_o2*an_r(3)**(ai_o2-1)) &
720: & * an_r(1)
722: d_eq(5,4) = d_eq(5,4) - (i_cc*an_r(4)**(i_cc-1))* &
723: & (an_r(5)**i_h2o)* ((pt/an_t)**idiff)
724: d_eq(5,5) = d_eq(5,5) &
725: & - (i_h2o*(an_r(5)**(i_h2o-1))) &
726: & * (an_r(4)**i_cc)* ((pt/an_t)**idiff)
728: d_eq(3,1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
729: do jj = 2,26
730: d_eq(3,jj) = d_eq(3,1)
731: enddo
733: d_eq(3,2) = d_eq(3,2) + k_eq(12)* (an_r(3)**3)
735: d_eq(3,3) = d_eq(3,3) + k_eq(12)* (3*an_r(3)**2)*an_r(2)
737: d_eq(3,4) = d_eq(3,4) - 2*an_r(4)*(an_r(5)**3)*(pt/an_t)
739: d_eq(3,5) = d_eq(3,5) - 3*(an_r(5)**2)*(an_r(4)**2)*(pt/an_t)
740: ! & *(pt_five/const_five)
742: do ii3 = 1,26
743: d_eq(3,ii3) = 0.0d0
744: enddo
745: d_eq(3,2) = 1.0d0
747: d_eq(7,1) = pt*an_r(11)*(-1.0)/const2 &
748: & -k_eq(1)*sqrt(pt)*sqrt(an_r(14)+1d-50)*(-0.5/const3)
750: do jj = 2,26
751: d_eq(7,jj) = d_eq(7,1)
752: enddo
754: d_eq(7,11) = d_eq(7,11) + pt/an_t
755: d_eq(7,14) = d_eq(7,14) &
756: & - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14)+1d-50)*an_t)))
758: d_eq(8,1) = pt*an_r(8)*(-1.0)/const2 &
759: & -k_eq(2)*sqrt(pt)*sqrt(an_r(3)+1.0d-50)*(-0.5/const3)
761: do jj = 2,26
762: d_eq(8,jj) = d_eq(8,1)
763: enddo
765: d_eq(8,3) = d_eq(8,3) &
766: & -k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3)+1.0d-50)*an_t)))
767: d_eq(8,8) = d_eq(8,8) + pt/an_t
769: d_eq(9,1) = pt*an_r(7)*(-1.0)/const2 &
770: & -k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3)
772: do jj = 2,26
773: d_eq(9,jj) = d_eq(9,1)
774: enddo
776: d_eq(9,7) = d_eq(9,7) + pt/an_t
777: d_eq(9,6) = d_eq(9,6) &
778: & -k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))
780: d_eq(10,1) = pt*an_r(10)*(-1.0)/const2 &
781: & -k_eq(4)*(pt)*sqrt((an_r(3)+1.0d-50) &
782: & *an_r(14))*(-1.0/const2)
783: do jj = 2,26
784: d_eq(10,jj) = d_eq(10,1)
785: enddo
787: d_eq(10,3) = d_eq(10,3) &
788: & -k_eq(4)*(pt)*sqrt(an_r(14)) &
789: & *(0.5/(sqrt(an_r(3)+1.0d-50)*an_t))
790: d_eq(10,10) = d_eq(10,10) + pt/an_t
791: d_eq(10,14) = d_eq(10,14) &
792: & -k_eq(4)*(pt)*sqrt(an_r(3)+1.0d-50) &
793: & *(0.5/(sqrt(an_r(14)+1.0d-50)*an_t))
795: d_eq(11,1) = pt*an_r(9)*(-1.0)/const2 &
796: & -k_eq(5)*(pt)*sqrt((an_r(3)+1.0d-50)*an_r(6)) &
797: & *(-1.0/const2)
799: do jj = 2,26
800: d_eq(11,jj) = d_eq(11,1)
801: enddo
803: d_eq(11,3) = d_eq(11,3) &
804: & -k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ &
805: & (sqrt(an_r(3)+1.0d-50)*an_t))
806: d_eq(11,6) = d_eq(11,6) &
807: & -k_eq(5)*(pt)*sqrt(an_r(3)+1.0d-50) &
808: & *(0.5/(sqrt(an_r(6))*an_t))
809: d_eq(11,9) = d_eq(11,9) + pt/an_t
811: d_eq(12,1) = pt*an_r(5)*(-1.0)/const2 &
812: & -k_eq(6)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
813: & *(an_r(14))*(-1.5/const5)
815: do jj = 2,26
816: d_eq(12,jj) = d_eq(12,1)
817: enddo
819: d_eq(12,3) = d_eq(12,3) &
820: & -k_eq(6)*(pt**1.5)*((an_r(14)+1.0d-50)/const3) &
821: & *(0.5/sqrt(an_r(3)+1.0d-50))
823: d_eq(12,5) = d_eq(12,5) + pt/an_t
824: d_eq(12,14) = d_eq(12,14) &
825: & -k_eq(6)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
827: d_eq(13,1) = pt*an_r(4)*(-1.0)/const2 &
828: & -k_eq(7)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
829: & *(an_r(13))*(-1.5/const5)
831: do jj = 2,26
832: d_eq(13,jj) = d_eq(13,1)
833: enddo
835: d_eq(13,3) = d_eq(13,3) &
836: & -k_eq(7)*(pt**1.5)*(an_r(13)/const3) &
837: & *(0.5/sqrt(an_r(3)+1.0d-50))
839: d_eq(13,4) = d_eq(13,4) + pt/an_t
840: d_eq(13,13) = d_eq(13,13) &
841: & -k_eq(7)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
843: d_eq(14,1) = pt*an_r(15)*(-1.0)/const2 &
844: & -k_eq(8)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
845: & *(an_r(9))*(-1.5/const5)
847: do jj = 2,26
848: d_eq(14,jj) = d_eq(14,1)
849: enddo
851: d_eq(14,3) = d_eq(14,3) &
852: & -k_eq(8)*(pt**1.5)*(an_r(9)/const3) &
853: & *(0.5/sqrt(an_r(3)+1.0d-50))
854: d_eq(14,9) = d_eq(14,9) &
855: & -k_eq(8)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
856: d_eq(14,15) = d_eq(14,15)+ pt/an_t
858: d_eq(15,1) = pt*an_r(16)*(-1.0)/const2 &
859: & -k_eq(9)*(pt**1.5)*sqrt(an_r(14)+1.0d-50) &
860: & *(an_r(3))*(-1.5/const5)
862: do jj = 2,26
863: d_eq(15,jj) = d_eq(15,1)
864: enddo
866: d_eq(15,3) = d_eq(15,3) &
867: & -k_eq(9)*(pt**1.5)*(sqrt(an_r(14)+1.0d-50)/const3)
868: d_eq(15,14) = d_eq(15,14) &
869: & -k_eq(9)*(pt**1.5)*(an_r(3)/const3) &
870: & *(0.5/sqrt(an_r(14)+1.0d-50))
871: d_eq(15,16) = d_eq(15,16) + pt/an_t
873: d_eq(16,1) = pt*an_r(12)*(-1.0)/const2 &
874: & -k_eq(10)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
875: & *(an_r(6))*(-1.5/const5)
877: do jj = 2,26
878: d_eq(16,jj) = d_eq(16,1)
879: enddo
881: d_eq(16,3) = d_eq(16,3) &
882: & -k_eq(10)*(pt**1.5)*(an_r(6)/const3) &
883: & *(0.5/sqrt(an_r(3)+1.0d-50))
885: d_eq(16,6) = d_eq(16,6) &
886: & -k_eq(10)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
887: d_eq(16,12) = d_eq(16,12) + pt/an_t
889: const_cube = an_t*an_t*an_t
890: const_four = const2*const2
892: d_eq(17,1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
893: & - k_eq(15) * an_r(17)*pt * (-1/const2)
894: do jj = 2,26
895: d_eq(17,jj) = d_eq(17,1)
896: enddo
897: d_eq(17,14) = d_eq(17,14) + an_r(18)*an_r(18)*(pt**3)/const_cube
898: d_eq(17,17) = d_eq(17,17) - k_eq(15)*pt/an_t
899: d_eq(17,18) = d_eq(17,18) + 2*an_r(18)*an_r(14) &
900: & *(pt**3)/const_cube
902: d_eq(18,1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube) &
903: & - k_eq(16) * an_r(3)*an_r(18)*an_r(18) &
904: & * (pt*pt*pt) * (-3/const_four)
905: do jj = 2,26
906: d_eq(18,jj) = d_eq(18,1)
907: enddo
908: d_eq(18,3) = d_eq(18,3) &
909: & - k_eq(16) *an_r(18)* an_r(18)*pt*pt*pt /const_cube
910: d_eq(18,13) = d_eq(18,13) &
911: & + 2* an_r(13)*pt*pt /const2
912: d_eq(18,18) = d_eq(18,18) -k_eq(16)*an_r(3) &
913: & * 2*an_r(18)*pt*pt*pt/const_cube
915: !====for eq 19
917: d_eq(19,1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube) &
918: & - k_eq(17)*an_r(13)*an_r(10)*pt*pt * (-2/const_cube)
919: do jj = 2,26
920: d_eq(19,jj) = d_eq(19,1)
921: enddo
922: d_eq(19,13) = d_eq(19,13) &
923: & - k_eq(17) *an_r(10)*pt*pt /const2
924: d_eq(19,10) = d_eq(19,10) &
925: & - k_eq(17) *an_r(13)*pt*pt /const2
926: d_eq(19,3) = d_eq(19,3) + an_r(19)*pt*pt/const2
927: d_eq(19,19) = d_eq(19,19) + an_r(3)*pt*pt/const2
928: !====for eq 20
930: d_eq(20,1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube) &
931: & - k_eq(18) * an_r(19)*an_r(8)*pt*pt * (-2/const_cube)
932: do jj = 2,26
933: d_eq(20,jj) = d_eq(20,1)
934: enddo
935: d_eq(20,8) = d_eq(20,8) &
936: & - k_eq(18) *an_r(19)*pt*pt /const2
937: d_eq(20,19) = d_eq(20,19) &
938: & - k_eq(18) *an_r(8)*pt*pt /const2
939: d_eq(20,20) = d_eq(20,20) + an_r(21)*pt*pt/const2
940: d_eq(20,21) = d_eq(20,21) + an_r(20)*pt*pt/const2
942: !========
943: !====for eq 21
945: d_eq(21,1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube) &
946: & - k_eq(19)*an_r(7)*an_r(8)*pt*pt * (-2/const_cube)
947: do jj = 2,26
948: d_eq(21,jj) = d_eq(21,1)
949: enddo
950: d_eq(21,7) = d_eq(21,7) &
951: & - k_eq(19) *an_r(8)*pt*pt /const2
952: d_eq(21,8) = d_eq(21,8) &
953: & - k_eq(19) *an_r(7)*pt*pt /const2
954: d_eq(21,21) = d_eq(21,21) + an_r(23)*pt*pt/const2
955: d_eq(21,23) = d_eq(21,23) + an_r(21)*pt*pt/const2
957: !========
958: ! for 22
959: d_eq(22,1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube) &
960: & -k_eq(20)*an_r(21)*an_r(22)*pt*pt * (-2/const_cube)
961: do jj = 2,26
962: d_eq(22,jj) = d_eq(22,1)
963: enddo
964: d_eq(22,21) = d_eq(22,21) &
965: & - k_eq(20) *an_r(22)*pt*pt /const2
966: d_eq(22,22) = d_eq(22,22) &
967: & - k_eq(20) *an_r(21)*pt*pt /const2
968: d_eq(22,11) = d_eq(22,11) + an_r(5)*pt*pt/(const2)
969: d_eq(22,5) = d_eq(22,5) + an_r(11)*pt*pt/(const2)
971: !========
972: ! for 23
974: d_eq(23,1) = an_r(24)*(pt)*(-1/const2) &
975: & - k_eq(21)*an_r(21)*an_r(3)*pt*pt * (-2/const_cube)
976: do jj = 2,26
977: d_eq(23,jj) = d_eq(23,1)
978: enddo
979: d_eq(23,3) = d_eq(23,3) &
980: & - k_eq(21) *an_r(21)*pt*pt /const2
981: d_eq(23,21) = d_eq(23,21) &
982: & - k_eq(21) *an_r(3)*pt*pt /const2
983: d_eq(23,24) = d_eq(23,24) + pt/(an_t)
985: !========
986: ! for 24
987: d_eq(24,1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube) &
988: & - k_eq(22)*an_r(24)*an_r(8)*pt*pt * (-2/const_cube)
989: do jj = 2,26
990: d_eq(24,jj) = d_eq(24,1)
991: enddo
992: d_eq(24,8) = d_eq(24,8) &
993: & - k_eq(22) *an_r(24)*pt*pt /const2
994: d_eq(24,24) = d_eq(24,24) &
995: & - k_eq(22) *an_r(8)*pt*pt /const2
996: d_eq(24,3) = d_eq(24,3) + an_r(25)*pt*pt/const2
997: d_eq(24,25) = d_eq(24,25) + an_r(3)*pt*pt/const2
999: !========
1000: !for 25
1002: d_eq(25,1) = an_r(26)*(pt)*(-1/const2) &
1003: & - k_eq(23)*an_r(21)*an_r(10)*pt*pt * (-2/const_cube)
1004: do jj = 2,26
1005: d_eq(25,jj) = d_eq(25,1)
1006: enddo
1007: d_eq(25,10) = d_eq(25,10) &
1008: & - k_eq(23) *an_r(21)*pt*pt /const2
1009: d_eq(25,21) = d_eq(25,21) &
1010: & - k_eq(23) *an_r(10)*pt*pt /const2
1011: d_eq(25,26) = d_eq(25,26) + pt/(an_t)
1013: !============
1014: ! for 26
1015: d_eq(26,20) = -1
1016: d_eq(26,22) = -1
1017: d_eq(26,23) = -1
1018: d_eq(26,21) = 1
1019: d_eq(26,24) = 1
1020: d_eq(26,25) = 1
1021: d_eq(26,26) = 1
1023: do j = 1,26
1024: do i = 1,26
1025: write(44,*)i,j,d_eq(i,j)
1026: enddo
1027: enddo
1029: return
1030: end
1032: subroutine ShashiPostCheck(ls,X,Y,W,c_Y,c_W,dummy)
1033: #include <petsc/finclude/petscsnes.h>
1034: use petscsnes
1035: implicit none
1036: SNESLineSearch ls
1037: PetscErrorCode ierr
1038: Vec X,Y,W
1039: PetscObject dummy
1040: PetscBool c_Y,c_W
1041: PetscScalar,pointer :: xx(:)
1042: PetscInt i
1043: PetscCall(VecGetArrayF90(W,xx,ierr))
1044: do i=1,26
1045: if (xx(i) < 0.0) then
1046: xx(i) = 0.0
1047: c_W = PETSC_TRUE
1048: endif
1049: if (xx(i) > 3.0) then
1050: xx(i) = 3.0
1051: endif
1052: enddo
1053: PetscCall(VecRestoreArrayF90(W,xx,ierr))
1054: return
1055: end