Actual source code: shashi.F90
2: program main
3: #include <petsc/finclude/petsc.h>
4: use petsc
5: implicit none
7: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
8: ! Variable declarations
9: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
10: !
11: ! Variables:
12: ! snes - nonlinear solver
13: ! ksp - linear solver
14: ! pc - preconditioner context
15: ! ksp - Krylov subspace method context
16: ! x, r - solution, residual vectors
17: ! J - Jacobian matrix
18: ! its - iterations for convergence
19: !
20: SNES snes
21: PC pc
22: KSP ksp
23: Vec x,r,lb,ub
24: Mat J
25: SNESLineSearch linesearch
26: PetscErrorCode ierr
27: PetscInt its,i2,i20
28: PetscMPIInt size
29: PetscScalar pfive
30: PetscReal tol
31: PetscBool setls
32: PetscScalar,pointer :: xx(:)
33: PetscScalar zero,big
34: SNESLineSearch ls
36: ! Note: Any user-defined Fortran routines (such as FormJacobian)
37: ! MUST be declared as external.
39: external FormFunction, FormJacobian
40: external ShashiPostCheck
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Macro definitions
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: !
46: ! Macros to make clearer the process of setting values in vectors and
47: ! getting values from vectors. These vectors are used in the routines
48: ! FormFunction() and FormJacobian().
49: ! - The element lx_a(ib) is element ib in the vector x
50: !
51: #define lx_a(ib) lx_v(lx_i + (ib))
52: #define lf_a(ib) lf_v(lf_i + (ib))
53: !
54: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: ! Beginning of program
56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
59: if (ierr .ne. 0) then
60: print*,'Unable to initialize PETSc'
61: stop
62: endif
63: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
64: if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,1,'requires one process'); endif
66: big = 2.88
67: big = PETSC_INFINITY
68: zero = 0.0
69: i2 = 26
70: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
71: ! Create nonlinear solver context
72: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
74: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: ! Create matrix and vector data structures; set corresponding routines
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: ! Create vectors for solution and nonlinear function
82: call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
83: call VecDuplicate(x,r,ierr)
85: ! Create Jacobian matrix data structure
87: call MatCreateDense(PETSC_COMM_SELF,26,26,26,26, &
88: & PETSC_NULL_SCALAR,J,ierr)
90: ! Set function evaluation routine and vector
92: call SNESSetFunction(snes,r,FormFunction,0,ierr)
94: ! Set Jacobian matrix data structure and Jacobian evaluation routine
96: call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)
98: call VecDuplicate(x,lb,ierr)
99: call VecDuplicate(x,ub,ierr)
100: call VecSet(lb,zero,ierr)
101: ! call VecGetArrayF90(lb,xx,ierr)
102: ! call ShashiLowerBound(xx)
103: ! call VecRestoreArrayF90(lb,xx,ierr)
104: call VecSet(ub,big,ierr)
105: ! call SNESVISetVariableBounds(snes,lb,ub,ierr)
107: call SNESGetLineSearch(snes,ls,ierr)
108: call SNESLineSearchSetPostCheck(ls,ShashiPostCheck, &
109: & 0,ierr)
110: call SNESSetType(snes,SNESVINEWTONRSLS,ierr)
112: call SNESSetFromOptions(snes,ierr)
114: ! set initial guess
116: call VecGetArrayF90(x,xx,ierr)
117: call ShashiInitialGuess(xx)
118: call VecRestoreArrayF90(x,xx,ierr)
120: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: ! Free work space. All PETSc objects should be destroyed when they
124: ! are no longer needed.
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: call VecDestroy(lb,ierr)
127: call VecDestroy(ub,ierr)
128: call VecDestroy(x,ierr)
129: call VecDestroy(r,ierr)
130: call MatDestroy(J,ierr)
131: call SNESDestroy(snes,ierr)
132: call PetscFinalize(ierr)
133: end
134: !
135: ! ------------------------------------------------------------------------
136: !
137: ! FormFunction - Evaluates nonlinear function, F(x).
138: !
139: ! Input Parameters:
140: ! snes - the SNES context
141: ! x - input vector
142: ! dummy - optional user-defined context (not used here)
143: !
144: ! Output Parameter:
145: ! f - function vector
146: !
147: subroutine FormFunction(snes,x,f,dummy,ierr)
148: #include <petsc/finclude/petscsnes.h>
149: use petscsnes
150: implicit none
151: SNES snes
152: Vec x,f
153: PetscErrorCode ierr
154: integer dummy(*)
156: ! Declarations for use with local arrays
158: PetscScalar lx_v(2),lf_v(2)
159: PetscOffset lx_i,lf_i
161: ! Get pointers to vector data.
162: ! - For default PETSc vectors, VecGetArray() returns a pointer to
163: ! the data array. Otherwise, the routine is implementation dependent.
164: ! - You MUST call VecRestoreArray() when you no longer need access to
165: ! the array.
166: ! - Note that the Fortran interface to VecGetArray() differs from the
167: ! C version. See the Fortran chapter of the users manual for details.
169: call VecGetArrayRead(x,lx_v,lx_i,ierr)
170: call VecGetArray(f,lf_v,lf_i,ierr)
171: call ShashiFormFunction(lx_a(1),lf_a(1))
172: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
173: call VecRestoreArray(f,lf_v,lf_i,ierr)
175: return
176: end
178: ! ---------------------------------------------------------------------
179: !
180: ! FormJacobian - Evaluates Jacobian matrix.
181: !
182: ! Input Parameters:
183: ! snes - the SNES context
184: ! x - input vector
185: ! dummy - optional user-defined context (not used here)
186: !
187: ! Output Parameters:
188: ! A - Jacobian matrix
189: ! B - optionally different preconditioning matrix
190: ! flag - flag indicating matrix structure
191: !
192: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
193: #include <petsc/finclude/petscsnes.h>
194: use petscsnes
195: implicit none
196: SNES snes
197: Vec X
198: Mat jac,B
199: PetscScalar A(4)
200: PetscErrorCode ierr
201: PetscInt idx(2),i2
202: integer dummy(*)
204: ! Declarations for use with local arrays
206: PetscScalar lx_v(1),lf_v(1)
207: PetscOffset lx_i,lf_i
209: ! Get pointer to vector data
211: call VecGetArrayRead(x,lx_v,lx_i,ierr)
212: call MatDenseGetArray(B,lf_v,lf_i,ierr)
213: call ShashiFormJacobian(lx_a(1),lf_a(1))
214: call MatDenseRestoreArray(B,lf_v,lf_i,ierr)
215: call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
217: ! Assemble matrix
219: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
220: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
222: return
223: end
225: subroutine ShashiLowerBound(an_r)
226: ! implicit PetscScalar (a-h,o-z)
227: implicit none
228: PetscScalar an_r(26)
229: PetscInt i
231: do i=2,26
232: an_r(i) = 1000.0/6.023D+23
233: enddo
234: return
235: end
237: subroutine ShashiInitialGuess(an_r)
238: ! implicit PetscScalar (a-h,o-z)
239: implicit none
240: PetscScalar an_c_additive
241: PetscScalar an_h_additive
242: PetscScalar an_o_additive
243: PetscScalar atom_c_init
244: PetscScalar atom_h_init
245: PetscScalar atom_n_init
246: PetscScalar atom_o_init
247: PetscScalar h_init
248: PetscScalar p_init
249: PetscInt nfuel
250: PetscScalar temp,pt
251: PetscScalar an_r(26),k_eq(23),f_eq(26)
252: PetscScalar d_eq(26,26),H_molar(26)
253: PetscInt an_h(1),an_c(1)
254: PetscScalar part_p(26)
255: PetscInt i_cc,i_hh,i_h2o
256: PetscInt i_pwr2,i_co2_h2o
258: pt = 0.1
259: atom_c_init =6.7408177364816552D-022
260: atom_h_init = 2.0
261: atom_o_init = 1.0
262: atom_n_init = 3.76
263: nfuel = 1
264: an_c(1) = 1
265: an_h(1) = 4
266: an_c_additive = 2
267: an_h_additive = 6
268: an_o_additive = 1
269: h_init = 128799.7267952987
270: p_init = 0.1
271: temp = 1500
273: an_r( 1) = 1.66000D-24
274: an_r( 2) = 1.66030D-22
275: an_r( 3) = 5.00000D-01
276: an_r( 4) = 1.66030D-22
277: an_r( 5) = 1.66030D-22
278: an_r( 6) = 1.88000D+00
279: an_r( 7) = 1.66030D-22
280: an_r( 8) = 1.66030D-22
281: an_r( 9) = 1.66030D-22
282: an_r(10) = 1.66030D-22
283: an_r(11) = 1.66030D-22
284: an_r(12) = 1.66030D-22
285: an_r(13) = 1.66030D-22
286: an_r(14) = 1.00000D+00
287: an_r(15) = 1.66030D-22
288: an_r(16) = 1.66030D-22
289: an_r(17) = 1.66000D-24
290: an_r(18) = 1.66030D-24
291: an_r(19) = 1.66030D-24
292: an_r(20) = 1.66030D-24
293: an_r(21) = 1.66030D-24
294: an_r(22) = 1.66030D-24
295: an_r(23) = 1.66030D-24
296: an_r(24) = 1.66030D-24
297: an_r(25) = 1.66030D-24
298: an_r(26) = 1.66030D-24
300: an_r = 0
301: an_r( 3) = .5
302: an_r( 6) = 1.88000
303: an_r(14) = 1.
305: #if defined(solution)
306: an_r( 1) = 3.802208D-33
307: an_r( 2) = 1.298287D-29
308: an_r( 3) = 2.533067D-04
309: an_r( 4) = 6.865078D-22
310: an_r( 5) = 9.993125D-01
311: an_r( 6) = 1.879964D+00
312: an_r( 7) = 4.449489D-13
313: an_r( 8) = 3.428687D-07
314: an_r( 9) = 7.105138D-05
315: an_r(10) = 1.094368D-04
316: an_r(11) = 2.362305D-06
317: an_r(12) = 1.107145D-09
318: an_r(13) = 1.276162D-24
319: an_r(14) = 6.315538D-04
320: an_r(15) = 2.356540D-09
321: an_r(16) = 2.048248D-09
322: an_r(17) = 1.966187D-22
323: an_r(18) = 7.856497D-29
324: an_r(19) = 1.987840D-36
325: an_r(20) = 8.182441D-22
326: an_r(21) = 2.684880D-16
327: an_r(22) = 2.680473D-16
328: an_r(23) = 6.594967D-18
329: an_r(24) = 2.509714D-21
330: an_r(25) = 3.096459D-21
331: an_r(26) = 6.149551D-18
332: #endif
334: return
335: end
337: subroutine ShashiFormFunction(an_r,f_eq)
338: ! implicit PetscScalar (a-h,o-z)
339: implicit none
340: PetscScalar an_c_additive
341: PetscScalar an_h_additive
342: PetscScalar an_o_additive
343: PetscScalar atom_c_init
344: PetscScalar atom_h_init
345: PetscScalar atom_n_init
346: PetscScalar atom_o_init
347: PetscScalar h_init
348: PetscScalar p_init
349: PetscInt nfuel
350: PetscScalar temp,pt
351: PetscScalar an_r(26),k_eq(23),f_eq(26)
352: PetscScalar d_eq(26,26),H_molar(26)
353: PetscInt an_h(1),an_c(1)
354: PetscScalar part_p(26),idiff
355: PetscInt i_cc,i_hh,i_h2o
356: PetscInt i_pwr2,i_co2_h2o
357: PetscScalar an_t,sum_h,pt_cubed,pt_five
358: PetscScalar pt_four,pt_val1,pt_val2
359: PetscScalar a_io2
360: PetscInt i,ip
361: pt = 0.1
362: atom_c_init =6.7408177364816552D-022
363: atom_h_init = 2.0
364: atom_o_init = 1.0
365: atom_n_init = 3.76
366: nfuel = 1
367: an_c(1) = 1
368: an_h(1) = 4
369: an_c_additive = 2
370: an_h_additive = 6
371: an_o_additive = 1
372: h_init = 128799.7267952987
373: p_init = 0.1
374: temp = 1500
376: k_eq( 1) = 1.75149D-05
377: k_eq( 2) = 4.01405D-06
378: k_eq( 3) = 6.04663D-14
379: k_eq( 4) = 2.73612D-01
380: k_eq( 5) = 3.25592D-03
381: k_eq( 6) = 5.33568D+05
382: k_eq( 7) = 2.07479D+05
383: k_eq( 8) = 1.11841D-02
384: k_eq( 9) = 1.72684D-03
385: k_eq(10) = 1.98588D-07
386: k_eq(11) = 7.23600D+27
387: k_eq(12) = 5.73926D+49
388: k_eq(13) = 1.00000D+00
389: k_eq(14) = 1.64493D+16
390: k_eq(15) = 2.73837D-29
391: k_eq(16) = 3.27419D+50
392: k_eq(17) = 1.72447D-23
393: k_eq(18) = 4.24657D-06
394: k_eq(19) = 1.16065D-14
395: k_eq(20) = 3.28020D+25
396: k_eq(21) = 1.06291D+00
397: k_eq(22) = 9.11507D+02
398: k_eq(23) = 6.02837D+03
400: H_molar( 1) = 3.26044D+03
401: H_molar( 2) = -8.00407D+04
402: H_molar( 3) = 4.05873D+04
403: H_molar( 4) = -3.31849D+05
404: H_molar( 5) = -1.93654D+05
405: H_molar( 6) = 3.84035D+04
406: H_molar( 7) = 4.97589D+05
407: H_molar( 8) = 2.74483D+05
408: H_molar( 9) = 1.30022D+05
409: H_molar(10) = 7.58429D+04
410: H_molar(11) = 2.42948D+05
411: H_molar(12) = 1.44588D+05
412: H_molar(13) = -7.16891D+04
413: H_molar(14) = 3.63075D+04
414: H_molar(15) = 9.23880D+04
415: H_molar(16) = 6.50477D+04
416: H_molar(17) = 3.04310D+05
417: H_molar(18) = 7.41707D+05
418: H_molar(19) = 6.32767D+05
419: H_molar(20) = 8.90624D+05
420: H_molar(21) = 2.49805D+04
421: H_molar(22) = 6.43473D+05
422: H_molar(23) = 1.02861D+06
423: H_molar(24) = -6.07503D+03
424: H_molar(25) = 1.27020D+05
425: H_molar(26) = -1.07011D+05
426: !=============
427: an_t = 0
428: sum_h = 0
430: do i = 1,26
431: an_t = an_t + an_r(i)
432: enddo
434: f_eq(1) = atom_h_init &
435: & - (an_h(1)*an_r(1) + an_h_additive*an_r(2) &
436: & + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) &
437: & + an_r(16)+ 2*an_r(17) + an_r(19) &
438: & +an_r(20) + 3*an_r(22)+an_r(26))
440: f_eq(2) = atom_o_init &
441: & - (an_o_additive*an_r(2) + 2*an_r(3) &
442: & + 2*an_r(4) + an_r(5) &
443: & + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
444: & + 2*an_r(15) + 2*an_r(16)+ an_r(20) + an_r(22) &
445: & + an_r(23) + 2*an_r(24) + 1*an_r(25)+an_r(26))
447: f_eq(3) = an_r(2)-1.0d-150
449: f_eq(4) = atom_c_init &
450: & - (an_c(1)*an_r(1) + an_c_additive * an_r(2) &
451: & + an_r(4) + an_r(13)+ 2*an_r(17) + an_r(18) &
452: & + an_r(19) + an_r(20))
454: do ip = 1,26
455: part_p(ip) = (an_r(ip)/an_t) * pt
456: enddo
458: i_cc = an_c(1)
459: i_hh = an_h(1)
460: a_io2 = i_cc + i_hh/4.0
461: i_h2o = i_hh/2
462: idiff = (i_cc + i_h2o) - (a_io2 + 1)
464: f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2 &
465: & - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
466: ! write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
467: ! stop
468: f_eq(6) = atom_n_init &
469: & - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) &
470: & + an_r(15) &
471: & + an_r(23))
473: f_eq( 7) = part_p(11) &
474: & - (k_eq( 1) * sqrt(part_p(14)+1d-23))
475: f_eq( 8) = part_p( 8) &
476: & - (k_eq( 2) * sqrt(part_p( 3)+1d-23))
478: f_eq( 9) = part_p( 7) &
479: & - (k_eq( 3) * sqrt(part_p( 6)+1d-23))
481: f_eq(10) = part_p(10) &
482: & - (k_eq( 4) * sqrt(part_p( 3)+1d-23)) &
483: & * sqrt(part_p(14))
485: f_eq(11) = part_p( 9) &
486: & - (k_eq( 5) * sqrt(part_p( 3)+1d-23)) &
487: & * sqrt(part_p( 6)+1d-23)
488: f_eq(12) = part_p( 5) &
489: & - (k_eq( 6) * sqrt(part_p( 3)+1d-23)) &
490: & * (part_p(14))
492: f_eq(13) = part_p( 4) &
493: & - (k_eq( 7) * sqrt(part_p(3)+1.0d-23)) &
494: & * (part_p(13))
496: f_eq(14) = part_p(15) &
497: & - (k_eq( 8) * sqrt(part_p(3)+1.0d-50)) &
498: & * (part_p( 9))
499: f_eq(15) = part_p(16) &
500: & - (k_eq( 9) * part_p( 3)) &
501: & * sqrt(part_p(14)+1d-23)
503: f_eq(16) = part_p(12) &
504: & - (k_eq(10) * sqrt(part_p( 3)+1d-23)) &
505: & * (part_p( 6))
507: f_eq(17) = part_p(14)*part_p(18)**2 &
508: & - (k_eq(15) * part_p(17))
510: f_eq(18) = (part_p(13)**2) &
511: & - (k_eq(16) * part_p(3)*part_p(18)**2)
512: print*,f_eq(18),part_p(3),part_p(18),part_p(13),k_eq(16)
514: f_eq(19) = part_p(19)*part_p(3) - k_eq(17)*part_p(13)*part_p(10)
516: f_eq(20) = part_p(21)*part_p(20) - k_eq(18)*part_p(19)*part_p(8)
518: f_eq(21) = part_p(21)*part_p(23) - k_eq(19)*part_p(7)*part_p(8)
520: f_eq(22) = part_p(5)*part_p(11) - k_eq(20)*part_p(21)*part_p(22)
522: f_eq(23) = part_p(24) - k_eq(21)*part_p(21)*part_p(3)
524: f_eq(24) = part_p(3)*part_p(25) - k_eq(22)*part_p(24)*part_p(8)
526: f_eq(25) = part_p(26) - k_eq(23)*part_p(21)*part_p(10)
528: f_eq(26) = -(an_r(20) + an_r(22) + an_r(23)) &
529: & +(an_r(21) + an_r(24) + an_r(25) + an_r(26))
531: do i = 1,26
532: write(44,*)i,f_eq(i)
533: enddo
535: return
536: end
538: subroutine ShashiFormJacobian(an_r,d_eq)
539: ! implicit PetscScalar (a-h,o-z)
540: implicit none
541: PetscScalar an_c_additive
542: PetscScalar an_h_additive
543: PetscScalar an_o_additive
544: PetscScalar atom_c_init
545: PetscScalar atom_h_init
546: PetscScalar atom_n_init
547: PetscScalar atom_o_init
548: PetscScalar h_init
549: PetscScalar p_init
550: PetscInt nfuel
551: PetscScalar temp,pt
552: PetscScalar an_t,ai_o2,sum_h
553: PetscScalar an_tot1_d,an_tot1
554: PetscScalar an_tot2_d,an_tot2
555: PetscScalar const5,const3,const2
556: PetscScalar const_cube
557: PetscScalar const_five
558: PetscScalar const_four
559: PetscScalar const_six
560: PetscInt jj,jb,ii3,id,ib,ip,i
561: PetscScalar pt2,pt1
562: PetscScalar an_r(26),k_eq(23),f_eq(26)
563: PetscScalar d_eq(26,26),H_molar(26)
564: PetscInt an_h(1),an_c(1)
565: PetscScalar ai_pwr1,part_p(26),idiff
566: PetscInt i_cc,i_hh
567: PetscInt i_h2o,i_pwr2,i_co2_h2o
568: PetscScalar pt_cube,pt_five
569: PetscScalar pt_four
570: PetscScalar pt_val1,pt_val2,a_io2
571: PetscInt j
573: pt = 0.1
574: atom_c_init =6.7408177364816552D-022
575: atom_h_init = 2.0
576: atom_o_init = 1.0
577: atom_n_init = 3.76
578: nfuel = 1
579: an_c(1) = 1
580: an_h(1) = 4
581: an_c_additive = 2
582: an_h_additive = 6
583: an_o_additive = 1
584: h_init = 128799.7267952987
585: p_init = 0.1
586: temp = 1500
588: k_eq( 1) = 1.75149D-05
589: k_eq( 2) = 4.01405D-06
590: k_eq( 3) = 6.04663D-14
591: k_eq( 4) = 2.73612D-01
592: k_eq( 5) = 3.25592D-03
593: k_eq( 6) = 5.33568D+05
594: k_eq( 7) = 2.07479D+05
595: k_eq( 8) = 1.11841D-02
596: k_eq( 9) = 1.72684D-03
597: k_eq(10) = 1.98588D-07
598: k_eq(11) = 7.23600D+27
599: k_eq(12) = 5.73926D+49
600: k_eq(13) = 1.00000D+00
601: k_eq(14) = 1.64493D+16
602: k_eq(15) = 2.73837D-29
603: k_eq(16) = 3.27419D+50
604: k_eq(17) = 1.72447D-23
605: k_eq(18) = 4.24657D-06
606: k_eq(19) = 1.16065D-14
607: k_eq(20) = 3.28020D+25
608: k_eq(21) = 1.06291D+00
609: k_eq(22) = 9.11507D+02
610: k_eq(23) = 6.02837D+03
612: H_molar( 1) = 3.26044D+03
613: H_molar( 2) = -8.00407D+04
614: H_molar( 3) = 4.05873D+04
615: H_molar( 4) = -3.31849D+05
616: H_molar( 5) = -1.93654D+05
617: H_molar( 6) = 3.84035D+04
618: H_molar( 7) = 4.97589D+05
619: H_molar( 8) = 2.74483D+05
620: H_molar( 9) = 1.30022D+05
621: H_molar(10) = 7.58429D+04
622: H_molar(11) = 2.42948D+05
623: H_molar(12) = 1.44588D+05
624: H_molar(13) = -7.16891D+04
625: H_molar(14) = 3.63075D+04
626: H_molar(15) = 9.23880D+04
627: H_molar(16) = 6.50477D+04
628: H_molar(17) = 3.04310D+05
629: H_molar(18) = 7.41707D+05
630: H_molar(19) = 6.32767D+05
631: H_molar(20) = 8.90624D+05
632: H_molar(21) = 2.49805D+04
633: H_molar(22) = 6.43473D+05
634: H_molar(23) = 1.02861D+06
635: H_molar(24) = -6.07503D+03
636: H_molar(25) = 1.27020D+05
637: H_molar(26) = -1.07011D+05
639: !
640: !=======
641: do jb = 1,26
642: do ib = 1,26
643: d_eq(ib,jb) = 0.0d0
644: end do
645: end do
647: an_t = 0.0
648: do id = 1,26
649: an_t = an_t + an_r(id)
650: enddo
652: !====
653: !====
654: d_eq(1,1) = -an_h(1)
655: d_eq(1,2) = -an_h_additive
656: d_eq(1,5) = -2
657: d_eq(1,10) = -1
658: d_eq(1,11) = -1
659: d_eq(1,14) = -2
660: d_eq(1,16) = -1
661: d_eq(1,17) = -2
662: d_eq(1,19) = -1
663: d_eq(1,20) = -1
664: d_eq(1,22) = -3
665: d_eq(1,26) = -1
667: d_eq(2,2) = -1*an_o_additive
668: d_eq(2,3) = -2
669: d_eq(2,4) = -2
670: d_eq(2,5) = -1
671: d_eq(2,8) = -1
672: d_eq(2,9) = -1
673: d_eq(2,10) = -1
674: d_eq(2,12) = -1
675: d_eq(2,13) = -1
676: d_eq(2,15) = -2
677: d_eq(2,16) = -2
678: d_eq(2,20) = -1
679: d_eq(2,22) = -1
680: d_eq(2,23) = -1
681: d_eq(2,24) = -2
682: d_eq(2,25) = -1
683: d_eq(2,26) = -1
685: d_eq(6,6) = -2
686: d_eq(6,7) = -1
687: d_eq(6,9) = -1
688: d_eq(6,12) = -2
689: d_eq(6,15) = -1
690: d_eq(6,23) = -1
692: d_eq(4,1) = -an_c(1)
693: d_eq(4,2) = -an_c_additive
694: d_eq(4,4) = -1
695: d_eq(4,13) = -1
696: d_eq(4,17) = -2
697: d_eq(4,18) = -1
698: d_eq(4,19) = -1
699: d_eq(4,20) = -1
701: !----------
702: const2 = an_t*an_t
703: const3 = (an_t)*sqrt(an_t)
704: const5 = an_t*const3
706: const_cube = an_t*an_t*an_t
707: const_four = const2*const2
708: const_five = const2*const_cube
709: const_six = const_cube*const_cube
710: pt_cube = pt*pt*pt
711: pt_four = pt_cube*pt
712: pt_five = pt_cube*pt*pt
714: i_cc = an_c(1)
715: i_hh = an_h(1)
716: ai_o2 = i_cc + float(i_hh)/4.0
717: i_co2_h2o = i_cc + i_hh/2
718: i_h2o = i_hh/2
719: ai_pwr1 = 1 + i_cc + float(i_hh)/4.0
720: i_pwr2 = i_cc + i_hh/2
721: idiff = (i_cc + i_h2o) - (ai_o2 + 1)
723: pt1 = pt**(ai_pwr1)
724: an_tot1 = an_t**(ai_pwr1)
725: pt_val1 = (pt/an_t)**(ai_pwr1)
726: an_tot1_d = an_tot1*an_t
728: pt2 = pt**(i_pwr2)
729: an_tot2 = an_t**(i_pwr2)
730: pt_val2 = (pt/an_t)**(i_pwr2)
731: an_tot2_d = an_tot2*an_t
733: d_eq(5,1) = &
734: & -(an_r(4)**i_cc)*(an_r(5)**i_h2o) &
735: & *((pt/an_t)**idiff) *(-idiff/an_t)
737: do jj = 2,26
738: d_eq(5,jj) = d_eq(5,1)
739: enddo
741: d_eq(5,1) = d_eq(5,1) + k_eq(11)* (an_r(3) **ai_o2)
743: d_eq(5,3) = d_eq(5,3) + k_eq(11)* (ai_o2*an_r(3)**(ai_o2-1)) &
744: & * an_r(1)
746: d_eq(5,4) = d_eq(5,4) - (i_cc*an_r(4)**(i_cc-1))* &
747: & (an_r(5)**i_h2o)* ((pt/an_t)**idiff)
748: d_eq(5,5) = d_eq(5,5) &
749: & - (i_h2o*(an_r(5)**(i_h2o-1))) &
750: & * (an_r(4)**i_cc)* ((pt/an_t)**idiff)
752: d_eq(3,1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
753: do jj = 2,26
754: d_eq(3,jj) = d_eq(3,1)
755: enddo
757: d_eq(3,2) = d_eq(3,2) + k_eq(12)* (an_r(3)**3)
759: d_eq(3,3) = d_eq(3,3) + k_eq(12)* (3*an_r(3)**2)*an_r(2)
761: d_eq(3,4) = d_eq(3,4) - 2*an_r(4)*(an_r(5)**3)*(pt/an_t)
763: d_eq(3,5) = d_eq(3,5) - 3*(an_r(5)**2)*(an_r(4)**2)*(pt/an_t)
764: ! & *(pt_five/const_five)
766: do ii3 = 1,26
767: d_eq(3,ii3) = 0.0d0
768: enddo
769: d_eq(3,2) = 1.0d0
771: d_eq(7,1) = pt*an_r(11)*(-1.0)/const2 &
772: & -k_eq(1)*sqrt(pt)*sqrt(an_r(14)+1d-50)*(-0.5/const3)
774: do jj = 2,26
775: d_eq(7,jj) = d_eq(7,1)
776: enddo
778: d_eq(7,11) = d_eq(7,11) + pt/an_t
779: d_eq(7,14) = d_eq(7,14) &
780: & - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14)+1d-50)*an_t)))
782: d_eq(8,1) = pt*an_r(8)*(-1.0)/const2 &
783: & -k_eq(2)*sqrt(pt)*sqrt(an_r(3)+1.0d-50)*(-0.5/const3)
785: do jj = 2,26
786: d_eq(8,jj) = d_eq(8,1)
787: enddo
789: d_eq(8,3) = d_eq(8,3) &
790: & -k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3)+1.0d-50)*an_t)))
791: d_eq(8,8) = d_eq(8,8) + pt/an_t
793: d_eq(9,1) = pt*an_r(7)*(-1.0)/const2 &
794: & -k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3)
796: do jj = 2,26
797: d_eq(9,jj) = d_eq(9,1)
798: enddo
800: d_eq(9,7) = d_eq(9,7) + pt/an_t
801: d_eq(9,6) = d_eq(9,6) &
802: & -k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))
804: d_eq(10,1) = pt*an_r(10)*(-1.0)/const2 &
805: & -k_eq(4)*(pt)*sqrt((an_r(3)+1.0d-50) &
806: & *an_r(14))*(-1.0/const2)
807: do jj = 2,26
808: d_eq(10,jj) = d_eq(10,1)
809: enddo
811: d_eq(10,3) = d_eq(10,3) &
812: & -k_eq(4)*(pt)*sqrt(an_r(14)) &
813: & *(0.5/(sqrt(an_r(3)+1.0d-50)*an_t))
814: d_eq(10,10) = d_eq(10,10) + pt/an_t
815: d_eq(10,14) = d_eq(10,14) &
816: & -k_eq(4)*(pt)*sqrt(an_r(3)+1.0d-50) &
817: & *(0.5/(sqrt(an_r(14)+1.0d-50)*an_t))
819: d_eq(11,1) = pt*an_r(9)*(-1.0)/const2 &
820: & -k_eq(5)*(pt)*sqrt((an_r(3)+1.0d-50)*an_r(6)) &
821: & *(-1.0/const2)
823: do jj = 2,26
824: d_eq(11,jj) = d_eq(11,1)
825: enddo
827: d_eq(11,3) = d_eq(11,3) &
828: & -k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ &
829: & (sqrt(an_r(3)+1.0d-50)*an_t))
830: d_eq(11,6) = d_eq(11,6) &
831: & -k_eq(5)*(pt)*sqrt(an_r(3)+1.0d-50) &
832: & *(0.5/(sqrt(an_r(6))*an_t))
833: d_eq(11,9) = d_eq(11,9) + pt/an_t
835: d_eq(12,1) = pt*an_r(5)*(-1.0)/const2 &
836: & -k_eq(6)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
837: & *(an_r(14))*(-1.5/const5)
839: do jj = 2,26
840: d_eq(12,jj) = d_eq(12,1)
841: enddo
843: d_eq(12,3) = d_eq(12,3) &
844: & -k_eq(6)*(pt**1.5)*((an_r(14)+1.0d-50)/const3) &
845: & *(0.5/sqrt(an_r(3)+1.0d-50))
847: d_eq(12,5) = d_eq(12,5) + pt/an_t
848: d_eq(12,14) = d_eq(12,14) &
849: & -k_eq(6)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
851: d_eq(13,1) = pt*an_r(4)*(-1.0)/const2 &
852: & -k_eq(7)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
853: & *(an_r(13))*(-1.5/const5)
855: do jj = 2,26
856: d_eq(13,jj) = d_eq(13,1)
857: enddo
859: d_eq(13,3) = d_eq(13,3) &
860: & -k_eq(7)*(pt**1.5)*(an_r(13)/const3) &
861: & *(0.5/sqrt(an_r(3)+1.0d-50))
863: d_eq(13,4) = d_eq(13,4) + pt/an_t
864: d_eq(13,13) = d_eq(13,13) &
865: & -k_eq(7)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
867: d_eq(14,1) = pt*an_r(15)*(-1.0)/const2 &
868: & -k_eq(8)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
869: & *(an_r(9))*(-1.5/const5)
871: do jj = 2,26
872: d_eq(14,jj) = d_eq(14,1)
873: enddo
875: d_eq(14,3) = d_eq(14,3) &
876: & -k_eq(8)*(pt**1.5)*(an_r(9)/const3) &
877: & *(0.5/sqrt(an_r(3)+1.0d-50))
878: d_eq(14,9) = d_eq(14,9) &
879: & -k_eq(8)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
880: d_eq(14,15) = d_eq(14,15)+ pt/an_t
882: d_eq(15,1) = pt*an_r(16)*(-1.0)/const2 &
883: & -k_eq(9)*(pt**1.5)*sqrt(an_r(14)+1.0d-50) &
884: & *(an_r(3))*(-1.5/const5)
886: do jj = 2,26
887: d_eq(15,jj) = d_eq(15,1)
888: enddo
890: d_eq(15,3) = d_eq(15,3) &
891: & -k_eq(9)*(pt**1.5)*(sqrt(an_r(14)+1.0d-50)/const3)
892: d_eq(15,14) = d_eq(15,14) &
893: & -k_eq(9)*(pt**1.5)*(an_r(3)/const3) &
894: & *(0.5/sqrt(an_r(14)+1.0d-50))
895: d_eq(15,16) = d_eq(15,16) + pt/an_t
897: d_eq(16,1) = pt*an_r(12)*(-1.0)/const2 &
898: & -k_eq(10)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
899: & *(an_r(6))*(-1.5/const5)
901: do jj = 2,26
902: d_eq(16,jj) = d_eq(16,1)
903: enddo
905: d_eq(16,3) = d_eq(16,3) &
906: & -k_eq(10)*(pt**1.5)*(an_r(6)/const3) &
907: & *(0.5/sqrt(an_r(3)+1.0d-50))
909: d_eq(16,6) = d_eq(16,6) &
910: & -k_eq(10)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
911: d_eq(16,12) = d_eq(16,12) + pt/an_t
913: const_cube = an_t*an_t*an_t
914: const_four = const2*const2
916: d_eq(17,1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
917: & - k_eq(15) * an_r(17)*pt * (-1/const2)
918: do jj = 2,26
919: d_eq(17,jj) = d_eq(17,1)
920: enddo
921: d_eq(17,14) = d_eq(17,14) + an_r(18)*an_r(18)*(pt**3)/const_cube
922: d_eq(17,17) = d_eq(17,17) - k_eq(15)*pt/an_t
923: d_eq(17,18) = d_eq(17,18) + 2*an_r(18)*an_r(14) &
924: & *(pt**3)/const_cube
926: d_eq(18,1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube) &
927: & - k_eq(16) * an_r(3)*an_r(18)*an_r(18) &
928: & * (pt*pt*pt) * (-3/const_four)
929: do jj = 2,26
930: d_eq(18,jj) = d_eq(18,1)
931: enddo
932: d_eq(18,3) = d_eq(18,3) &
933: & - k_eq(16) *an_r(18)* an_r(18)*pt*pt*pt /const_cube
934: d_eq(18,13) = d_eq(18,13) &
935: & + 2* an_r(13)*pt*pt /const2
936: d_eq(18,18) = d_eq(18,18) -k_eq(16)*an_r(3) &
937: & * 2*an_r(18)*pt*pt*pt/const_cube
939: !====for eq 19
941: d_eq(19,1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube) &
942: & - k_eq(17)*an_r(13)*an_r(10)*pt*pt * (-2/const_cube)
943: do jj = 2,26
944: d_eq(19,jj) = d_eq(19,1)
945: enddo
946: d_eq(19,13) = d_eq(19,13) &
947: & - k_eq(17) *an_r(10)*pt*pt /const2
948: d_eq(19,10) = d_eq(19,10) &
949: & - k_eq(17) *an_r(13)*pt*pt /const2
950: d_eq(19,3) = d_eq(19,3) + an_r(19)*pt*pt/const2
951: d_eq(19,19) = d_eq(19,19) + an_r(3)*pt*pt/const2
952: !====for eq 20
954: d_eq(20,1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube) &
955: & - k_eq(18) * an_r(19)*an_r(8)*pt*pt * (-2/const_cube)
956: do jj = 2,26
957: d_eq(20,jj) = d_eq(20,1)
958: enddo
959: d_eq(20,8) = d_eq(20,8) &
960: & - k_eq(18) *an_r(19)*pt*pt /const2
961: d_eq(20,19) = d_eq(20,19) &
962: & - k_eq(18) *an_r(8)*pt*pt /const2
963: d_eq(20,20) = d_eq(20,20) + an_r(21)*pt*pt/const2
964: d_eq(20,21) = d_eq(20,21) + an_r(20)*pt*pt/const2
966: !========
967: !====for eq 21
969: d_eq(21,1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube) &
970: & - k_eq(19)*an_r(7)*an_r(8)*pt*pt * (-2/const_cube)
971: do jj = 2,26
972: d_eq(21,jj) = d_eq(21,1)
973: enddo
974: d_eq(21,7) = d_eq(21,7) &
975: & - k_eq(19) *an_r(8)*pt*pt /const2
976: d_eq(21,8) = d_eq(21,8) &
977: & - k_eq(19) *an_r(7)*pt*pt /const2
978: d_eq(21,21) = d_eq(21,21) + an_r(23)*pt*pt/const2
979: d_eq(21,23) = d_eq(21,23) + an_r(21)*pt*pt/const2
981: !========
982: ! for 22
983: d_eq(22,1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube) &
984: & -k_eq(20)*an_r(21)*an_r(22)*pt*pt * (-2/const_cube)
985: do jj = 2,26
986: d_eq(22,jj) = d_eq(22,1)
987: enddo
988: d_eq(22,21) = d_eq(22,21) &
989: & - k_eq(20) *an_r(22)*pt*pt /const2
990: d_eq(22,22) = d_eq(22,22) &
991: & - k_eq(20) *an_r(21)*pt*pt /const2
992: d_eq(22,11) = d_eq(22,11) + an_r(5)*pt*pt/(const2)
993: d_eq(22,5) = d_eq(22,5) + an_r(11)*pt*pt/(const2)
995: !========
996: ! for 23
998: d_eq(23,1) = an_r(24)*(pt)*(-1/const2) &
999: & - k_eq(21)*an_r(21)*an_r(3)*pt*pt * (-2/const_cube)
1000: do jj = 2,26
1001: d_eq(23,jj) = d_eq(23,1)
1002: enddo
1003: d_eq(23,3) = d_eq(23,3) &
1004: & - k_eq(21) *an_r(21)*pt*pt /const2
1005: d_eq(23,21) = d_eq(23,21) &
1006: & - k_eq(21) *an_r(3)*pt*pt /const2
1007: d_eq(23,24) = d_eq(23,24) + pt/(an_t)
1009: !========
1010: ! for 24
1011: d_eq(24,1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube) &
1012: & - k_eq(22)*an_r(24)*an_r(8)*pt*pt * (-2/const_cube)
1013: do jj = 2,26
1014: d_eq(24,jj) = d_eq(24,1)
1015: enddo
1016: d_eq(24,8) = d_eq(24,8) &
1017: & - k_eq(22) *an_r(24)*pt*pt /const2
1018: d_eq(24,24) = d_eq(24,24) &
1019: & - k_eq(22) *an_r(8)*pt*pt /const2
1020: d_eq(24,3) = d_eq(24,3) + an_r(25)*pt*pt/const2
1021: d_eq(24,25) = d_eq(24,25) + an_r(3)*pt*pt/const2
1023: !========
1024: !for 25
1026: d_eq(25,1) = an_r(26)*(pt)*(-1/const2) &
1027: & - k_eq(23)*an_r(21)*an_r(10)*pt*pt * (-2/const_cube)
1028: do jj = 2,26
1029: d_eq(25,jj) = d_eq(25,1)
1030: enddo
1031: d_eq(25,10) = d_eq(25,10) &
1032: & - k_eq(23) *an_r(21)*pt*pt /const2
1033: d_eq(25,21) = d_eq(25,21) &
1034: & - k_eq(23) *an_r(10)*pt*pt /const2
1035: d_eq(25,26) = d_eq(25,26) + pt/(an_t)
1037: !============
1038: ! for 26
1039: d_eq(26,20) = -1
1040: d_eq(26,22) = -1
1041: d_eq(26,23) = -1
1042: d_eq(26,21) = 1
1043: d_eq(26,24) = 1
1044: d_eq(26,25) = 1
1045: d_eq(26,26) = 1
1047: do j = 1,26
1048: do i = 1,26
1049: write(44,*)i,j,d_eq(i,j)
1050: enddo
1051: enddo
1053: return
1054: end
1056: subroutine ShashiPostCheck(ls,X,Y,W,c_Y,c_W,dummy)
1057: #include <petsc/finclude/petscsnes.h>
1058: use petscsnes
1059: implicit none
1060: SNESLineSearch ls
1061: PetscErrorCode ierr
1062: Vec X,Y,W
1063: PetscObject dummy
1064: PetscBool c_Y,c_W
1065: PetscScalar,pointer :: xx(:)
1066: PetscInt i
1067: call VecGetArrayF90(W,xx,ierr)
1068: do i=1,26
1069: if (xx(i) < 0.0) then
1070: xx(i) = 0.0
1071: c_W = PETSC_TRUE
1072: endif
1073: if (xx(i) > 3.0) then
1074: xx(i) = 3.0
1075: endif
1076: enddo
1077: call VecRestoreArrayF90(W,xx,ierr)
1078: return
1079: end