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