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 .eq. 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(VecGetArrayF90(lb,xx,ierr))
75: ! PetscCallA(ShashiLowerBound(xx))
76: ! PetscCallA(VecRestoreArrayF90(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(VecGetArrayF90(x,xx,ierr))
89: PetscCallA(ShashiInitialGuess(xx))
90: PetscCallA(VecRestoreArrayF90(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(VecGetArrayReadF90(x,lx_v,ierr))
141: PetscCall(VecGetArrayF90(f,lf_v,ierr))
142: PetscCall(ShashiFormFunction(lx_v,lf_v))
143: PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
144: PetscCall(VecRestoreArrayF90(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 preconditioning matrix
160: ! flag - flag indicating matrix structure
161: !
162: subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
163: #include <petsc/finclude/petscsnes.h>
164: use petscsnes
165: implicit none
166: SNES snes
167: Vec X
168: Mat jac,B
169: PetscErrorCode ierr
170: integer dummy(*)
172: ! Declarations for use with local arrays
173: PetscScalar,pointer ::lx_v(:),lf_v(:,:)
175: ! Get pointer to vector data
177: PetscCall(VecGetArrayReadF90(x,lx_v,ierr))
178: PetscCall(MatDenseGetArrayF90(B,lf_v,ierr))
179: PetscCall(ShashiFormJacobian(lx_v,lf_v))
180: PetscCall(MatDenseRestoreArrayF90(B,lf_v,ierr))
181: PetscCall(VecRestoreArrayReadF90(x,lx_v,ierr))
183: ! Assemble matrix
185: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
186: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
188: end
190: subroutine ShashiLowerBound(an_r)
191: ! implicit PetscScalar (a-h,o-z)
192: implicit none
193: PetscScalar an_r(26)
194: PetscInt i
196: do i=2,26
197: an_r(i) = 1000.0/6.023D+23
198: enddo
199: end
201: subroutine ShashiInitialGuess(an_r)
202: ! implicit PetscScalar (a-h,o-z)
203: implicit none
204: PetscScalar an_c_additive
205: PetscScalar an_h_additive
206: PetscScalar an_o_additive
207: PetscScalar atom_c_init
208: PetscScalar atom_h_init
209: PetscScalar atom_n_init
210: PetscScalar atom_o_init
211: PetscScalar h_init
212: PetscScalar p_init
213: PetscInt nfuel
214: PetscScalar temp,pt
215: PetscScalar an_r(26)
216: PetscInt an_h(1),an_c(1)
218: pt = 0.1
219: atom_c_init =6.7408177364816552D-022
220: atom_h_init = 2.0
221: atom_o_init = 1.0
222: atom_n_init = 3.76
223: nfuel = 1
224: an_c(1) = 1
225: an_h(1) = 4
226: an_c_additive = 2
227: an_h_additive = 6
228: an_o_additive = 1
229: h_init = 128799.7267952987
230: p_init = 0.1
231: temp = 1500
233: an_r( 1) = 1.66000D-24
234: an_r( 2) = 1.66030D-22
235: an_r( 3) = 5.00000D-01
236: an_r( 4) = 1.66030D-22
237: an_r( 5) = 1.66030D-22
238: an_r( 6) = 1.88000D+00
239: an_r( 7) = 1.66030D-22
240: an_r( 8) = 1.66030D-22
241: an_r( 9) = 1.66030D-22
242: an_r(10) = 1.66030D-22
243: an_r(11) = 1.66030D-22
244: an_r(12) = 1.66030D-22
245: an_r(13) = 1.66030D-22
246: an_r(14) = 1.00000D+00
247: an_r(15) = 1.66030D-22
248: an_r(16) = 1.66030D-22
249: an_r(17) = 1.66000D-24
250: an_r(18) = 1.66030D-24
251: an_r(19) = 1.66030D-24
252: an_r(20) = 1.66030D-24
253: an_r(21) = 1.66030D-24
254: an_r(22) = 1.66030D-24
255: an_r(23) = 1.66030D-24
256: an_r(24) = 1.66030D-24
257: an_r(25) = 1.66030D-24
258: an_r(26) = 1.66030D-24
260: an_r = 0
261: an_r( 3) = .5
262: an_r( 6) = 1.88000
263: an_r(14) = 1.
265: #if defined(solution)
266: an_r( 1) = 3.802208D-33
267: an_r( 2) = 1.298287D-29
268: an_r( 3) = 2.533067D-04
269: an_r( 4) = 6.865078D-22
270: an_r( 5) = 9.993125D-01
271: an_r( 6) = 1.879964D+00
272: an_r( 7) = 4.449489D-13
273: an_r( 8) = 3.428687D-07
274: an_r( 9) = 7.105138D-05
275: an_r(10) = 1.094368D-04
276: an_r(11) = 2.362305D-06
277: an_r(12) = 1.107145D-09
278: an_r(13) = 1.276162D-24
279: an_r(14) = 6.315538D-04
280: an_r(15) = 2.356540D-09
281: an_r(16) = 2.048248D-09
282: an_r(17) = 1.966187D-22
283: an_r(18) = 7.856497D-29
284: an_r(19) = 1.987840D-36
285: an_r(20) = 8.182441D-22
286: an_r(21) = 2.684880D-16
287: an_r(22) = 2.680473D-16
288: an_r(23) = 6.594967D-18
289: an_r(24) = 2.509714D-21
290: an_r(25) = 3.096459D-21
291: an_r(26) = 6.149551D-18
292: #endif
294: end
296: subroutine ShashiFormFunction(an_r,f_eq)
297: ! implicit PetscScalar (a-h,o-z)
298: implicit none
299: PetscScalar an_c_additive
300: PetscScalar an_h_additive
301: PetscScalar an_o_additive
302: PetscScalar atom_c_init
303: PetscScalar atom_h_init
304: PetscScalar atom_n_init
305: PetscScalar atom_o_init
306: PetscScalar h_init
307: PetscScalar p_init
308: PetscInt nfuel
309: PetscScalar temp,pt
310: PetscScalar an_r(26),k_eq(23),f_eq(26)
311: PetscScalar H_molar(26)
312: PetscInt an_h(1),an_c(1)
313: PetscScalar part_p(26),idiff
314: PetscInt i_cc,i_hh,i_h2o
315: PetscScalar an_t,sum_h
316: PetscScalar a_io2
317: PetscInt i,ip
318: pt = 0.1
319: atom_c_init =6.7408177364816552D-022
320: atom_h_init = 2.0
321: atom_o_init = 1.0
322: atom_n_init = 3.76
323: nfuel = 1
324: an_c(1) = 1
325: an_h(1) = 4
326: an_c_additive = 2
327: an_h_additive = 6
328: an_o_additive = 1
329: h_init = 128799.7267952987
330: p_init = 0.1
331: temp = 1500
333: k_eq( 1) = 1.75149D-05
334: k_eq( 2) = 4.01405D-06
335: k_eq( 3) = 6.04663D-14
336: k_eq( 4) = 2.73612D-01
337: k_eq( 5) = 3.25592D-03
338: k_eq( 6) = 5.33568D+05
339: k_eq( 7) = 2.07479D+05
340: k_eq( 8) = 1.11841D-02
341: k_eq( 9) = 1.72684D-03
342: k_eq(10) = 1.98588D-07
343: k_eq(11) = 7.23600D+27
344: k_eq(12) = 5.73926D+49
345: k_eq(13) = 1.00000D+00
346: k_eq(14) = 1.64493D+16
347: k_eq(15) = 2.73837D-29
348: k_eq(16) = 3.27419D+50
349: k_eq(17) = 1.72447D-23
350: k_eq(18) = 4.24657D-06
351: k_eq(19) = 1.16065D-14
352: k_eq(20) = 3.28020D+25
353: k_eq(21) = 1.06291D+00
354: k_eq(22) = 9.11507D+02
355: k_eq(23) = 6.02837D+03
357: H_molar( 1) = 3.26044D+03
358: H_molar( 2) = -8.00407D+04
359: H_molar( 3) = 4.05873D+04
360: H_molar( 4) = -3.31849D+05
361: H_molar( 5) = -1.93654D+05
362: H_molar( 6) = 3.84035D+04
363: H_molar( 7) = 4.97589D+05
364: H_molar( 8) = 2.74483D+05
365: H_molar( 9) = 1.30022D+05
366: H_molar(10) = 7.58429D+04
367: H_molar(11) = 2.42948D+05
368: H_molar(12) = 1.44588D+05
369: H_molar(13) = -7.16891D+04
370: H_molar(14) = 3.63075D+04
371: H_molar(15) = 9.23880D+04
372: H_molar(16) = 6.50477D+04
373: H_molar(17) = 3.04310D+05
374: H_molar(18) = 7.41707D+05
375: H_molar(19) = 6.32767D+05
376: H_molar(20) = 8.90624D+05
377: H_molar(21) = 2.49805D+04
378: H_molar(22) = 6.43473D+05
379: H_molar(23) = 1.02861D+06
380: H_molar(24) = -6.07503D+03
381: H_molar(25) = 1.27020D+05
382: H_molar(26) = -1.07011D+05
383: !=============
384: an_t = 0
385: sum_h = 0
387: do i = 1,26
388: an_t = an_t + an_r(i)
389: enddo
391: f_eq(1) = atom_h_init &
392: & - (an_h(1)*an_r(1) + an_h_additive*an_r(2) &
393: & + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) &
394: & + an_r(16)+ 2*an_r(17) + an_r(19) &
395: & +an_r(20) + 3*an_r(22)+an_r(26))
397: f_eq(2) = atom_o_init &
398: & - (an_o_additive*an_r(2) + 2*an_r(3) &
399: & + 2*an_r(4) + an_r(5) &
400: & + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
401: & + 2*an_r(15) + 2*an_r(16)+ an_r(20) + an_r(22) &
402: & + an_r(23) + 2*an_r(24) + 1*an_r(25)+an_r(26))
404: f_eq(3) = an_r(2)-1.0d-150
406: f_eq(4) = atom_c_init &
407: & - (an_c(1)*an_r(1) + an_c_additive * an_r(2) &
408: & + an_r(4) + an_r(13)+ 2*an_r(17) + an_r(18) &
409: & + an_r(19) + an_r(20))
411: do ip = 1,26
412: part_p(ip) = (an_r(ip)/an_t) * pt
413: enddo
415: i_cc = an_c(1)
416: i_hh = an_h(1)
417: a_io2 = i_cc + i_hh/4.0
418: i_h2o = i_hh/2
419: idiff = (i_cc + i_h2o) - (a_io2 + 1)
421: f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2 &
422: & - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
423: ! write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
424: ! stop
425: f_eq(6) = atom_n_init &
426: & - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) &
427: & + an_r(15) &
428: & + an_r(23))
430: f_eq( 7) = part_p(11) &
431: & - (k_eq( 1) * sqrt(part_p(14)+1d-23))
432: f_eq( 8) = part_p( 8) &
433: & - (k_eq( 2) * sqrt(part_p( 3)+1d-23))
435: f_eq( 9) = part_p( 7) &
436: & - (k_eq( 3) * sqrt(part_p( 6)+1d-23))
438: f_eq(10) = part_p(10) &
439: & - (k_eq( 4) * sqrt(part_p( 3)+1d-23)) &
440: & * sqrt(part_p(14))
442: f_eq(11) = part_p( 9) &
443: & - (k_eq( 5) * sqrt(part_p( 3)+1d-23)) &
444: & * sqrt(part_p( 6)+1d-23)
445: f_eq(12) = part_p( 5) &
446: & - (k_eq( 6) * sqrt(part_p( 3)+1d-23)) &
447: & * (part_p(14))
449: f_eq(13) = part_p( 4) &
450: & - (k_eq( 7) * sqrt(part_p(3)+1.0d-23)) &
451: & * (part_p(13))
453: f_eq(14) = part_p(15) &
454: & - (k_eq( 8) * sqrt(part_p(3)+1.0d-50)) &
455: & * (part_p( 9))
456: f_eq(15) = part_p(16) &
457: & - (k_eq( 9) * part_p( 3)) &
458: & * sqrt(part_p(14)+1d-23)
460: f_eq(16) = part_p(12) &
461: & - (k_eq(10) * sqrt(part_p( 3)+1d-23)) &
462: & * (part_p( 6))
464: f_eq(17) = part_p(14)*part_p(18)**2 &
465: & - (k_eq(15) * part_p(17))
467: f_eq(18) = (part_p(13)**2) &
468: & - (k_eq(16) * part_p(3)*part_p(18)**2)
469: print*,f_eq(18),part_p(3),part_p(18),part_p(13),k_eq(16)
471: f_eq(19) = part_p(19)*part_p(3) - k_eq(17)*part_p(13)*part_p(10)
473: f_eq(20) = part_p(21)*part_p(20) - k_eq(18)*part_p(19)*part_p(8)
475: f_eq(21) = part_p(21)*part_p(23) - k_eq(19)*part_p(7)*part_p(8)
477: f_eq(22) = part_p(5)*part_p(11) - k_eq(20)*part_p(21)*part_p(22)
479: f_eq(23) = part_p(24) - k_eq(21)*part_p(21)*part_p(3)
481: f_eq(24) = part_p(3)*part_p(25) - k_eq(22)*part_p(24)*part_p(8)
483: f_eq(25) = part_p(26) - k_eq(23)*part_p(21)*part_p(10)
485: f_eq(26) = -(an_r(20) + an_r(22) + an_r(23)) &
486: & +(an_r(21) + an_r(24) + an_r(25) + an_r(26))
488: do i = 1,26
489: write(44,*)i,f_eq(i)
490: enddo
492: end
494: subroutine ShashiFormJacobian(an_r,d_eq)
495: ! implicit PetscScalar (a-h,o-z)
496: implicit none
497: PetscScalar an_c_additive
498: PetscScalar an_h_additive
499: PetscScalar an_o_additive
500: PetscScalar atom_c_init
501: PetscScalar atom_h_init
502: PetscScalar atom_n_init
503: PetscScalar atom_o_init
504: PetscScalar h_init
505: PetscScalar p_init
506: PetscInt nfuel
507: PetscScalar temp,pt
508: PetscScalar an_t,ai_o2
509: PetscScalar an_tot1_d,an_tot1
510: PetscScalar an_tot2_d,an_tot2
511: PetscScalar const5,const3,const2
512: PetscScalar const_cube
513: PetscScalar const_five
514: PetscScalar const_four
515: PetscScalar const_six
516: PetscInt jj,jb,ii3,id,ib,i
517: PetscScalar pt2,pt1
518: PetscScalar an_r(26),k_eq(23)
519: PetscScalar d_eq(26,26),H_molar(26)
520: PetscInt an_h(1),an_c(1)
521: PetscScalar ai_pwr1,idiff
522: PetscInt i_cc,i_hh
523: PetscInt i_h2o,i_pwr2,i_co2_h2o
524: PetscScalar pt_cube,pt_five
525: PetscScalar pt_four
526: PetscScalar pt_val1,pt_val2
527: PetscInt j
529: pt = 0.1
530: atom_c_init =6.7408177364816552D-022
531: atom_h_init = 2.0
532: atom_o_init = 1.0
533: atom_n_init = 3.76
534: nfuel = 1
535: an_c(1) = 1
536: an_h(1) = 4
537: an_c_additive = 2
538: an_h_additive = 6
539: an_o_additive = 1
540: h_init = 128799.7267952987
541: p_init = 0.1
542: temp = 1500
544: k_eq( 1) = 1.75149D-05
545: k_eq( 2) = 4.01405D-06
546: k_eq( 3) = 6.04663D-14
547: k_eq( 4) = 2.73612D-01
548: k_eq( 5) = 3.25592D-03
549: k_eq( 6) = 5.33568D+05
550: k_eq( 7) = 2.07479D+05
551: k_eq( 8) = 1.11841D-02
552: k_eq( 9) = 1.72684D-03
553: k_eq(10) = 1.98588D-07
554: k_eq(11) = 7.23600D+27
555: k_eq(12) = 5.73926D+49
556: k_eq(13) = 1.00000D+00
557: k_eq(14) = 1.64493D+16
558: k_eq(15) = 2.73837D-29
559: k_eq(16) = 3.27419D+50
560: k_eq(17) = 1.72447D-23
561: k_eq(18) = 4.24657D-06
562: k_eq(19) = 1.16065D-14
563: k_eq(20) = 3.28020D+25
564: k_eq(21) = 1.06291D+00
565: k_eq(22) = 9.11507D+02
566: k_eq(23) = 6.02837D+03
568: H_molar( 1) = 3.26044D+03
569: H_molar( 2) = -8.00407D+04
570: H_molar( 3) = 4.05873D+04
571: H_molar( 4) = -3.31849D+05
572: H_molar( 5) = -1.93654D+05
573: H_molar( 6) = 3.84035D+04
574: H_molar( 7) = 4.97589D+05
575: H_molar( 8) = 2.74483D+05
576: H_molar( 9) = 1.30022D+05
577: H_molar(10) = 7.58429D+04
578: H_molar(11) = 2.42948D+05
579: H_molar(12) = 1.44588D+05
580: H_molar(13) = -7.16891D+04
581: H_molar(14) = 3.63075D+04
582: H_molar(15) = 9.23880D+04
583: H_molar(16) = 6.50477D+04
584: H_molar(17) = 3.04310D+05
585: H_molar(18) = 7.41707D+05
586: H_molar(19) = 6.32767D+05
587: H_molar(20) = 8.90624D+05
588: H_molar(21) = 2.49805D+04
589: H_molar(22) = 6.43473D+05
590: H_molar(23) = 1.02861D+06
591: H_molar(24) = -6.07503D+03
592: H_molar(25) = 1.27020D+05
593: H_molar(26) = -1.07011D+05
595: !
596: !=======
597: do jb = 1,26
598: do ib = 1,26
599: d_eq(ib,jb) = 0.0d0
600: end do
601: end do
603: an_t = 0.0
604: do id = 1,26
605: an_t = an_t + an_r(id)
606: enddo
608: !====
609: !====
610: d_eq(1,1) = -an_h(1)
611: d_eq(1,2) = -an_h_additive
612: d_eq(1,5) = -2
613: d_eq(1,10) = -1
614: d_eq(1,11) = -1
615: d_eq(1,14) = -2
616: d_eq(1,16) = -1
617: d_eq(1,17) = -2
618: d_eq(1,19) = -1
619: d_eq(1,20) = -1
620: d_eq(1,22) = -3
621: d_eq(1,26) = -1
623: d_eq(2,2) = -1*an_o_additive
624: d_eq(2,3) = -2
625: d_eq(2,4) = -2
626: d_eq(2,5) = -1
627: d_eq(2,8) = -1
628: d_eq(2,9) = -1
629: d_eq(2,10) = -1
630: d_eq(2,12) = -1
631: d_eq(2,13) = -1
632: d_eq(2,15) = -2
633: d_eq(2,16) = -2
634: d_eq(2,20) = -1
635: d_eq(2,22) = -1
636: d_eq(2,23) = -1
637: d_eq(2,24) = -2
638: d_eq(2,25) = -1
639: d_eq(2,26) = -1
641: d_eq(6,6) = -2
642: d_eq(6,7) = -1
643: d_eq(6,9) = -1
644: d_eq(6,12) = -2
645: d_eq(6,15) = -1
646: d_eq(6,23) = -1
648: d_eq(4,1) = -an_c(1)
649: d_eq(4,2) = -an_c_additive
650: d_eq(4,4) = -1
651: d_eq(4,13) = -1
652: d_eq(4,17) = -2
653: d_eq(4,18) = -1
654: d_eq(4,19) = -1
655: d_eq(4,20) = -1
657: !----------
658: const2 = an_t*an_t
659: const3 = (an_t)*sqrt(an_t)
660: const5 = an_t*const3
662: const_cube = an_t*an_t*an_t
663: const_four = const2*const2
664: const_five = const2*const_cube
665: const_six = const_cube*const_cube
666: pt_cube = pt*pt*pt
667: pt_four = pt_cube*pt
668: pt_five = pt_cube*pt*pt
670: i_cc = an_c(1)
671: i_hh = an_h(1)
672: ai_o2 = i_cc + float(i_hh)/4.0
673: i_co2_h2o = i_cc + i_hh/2
674: i_h2o = i_hh/2
675: ai_pwr1 = 1 + i_cc + float(i_hh)/4.0
676: i_pwr2 = i_cc + i_hh/2
677: idiff = (i_cc + i_h2o) - (ai_o2 + 1)
679: pt1 = pt**(ai_pwr1)
680: an_tot1 = an_t**(ai_pwr1)
681: pt_val1 = (pt/an_t)**(ai_pwr1)
682: an_tot1_d = an_tot1*an_t
684: pt2 = pt**(i_pwr2)
685: an_tot2 = an_t**(i_pwr2)
686: pt_val2 = (pt/an_t)**(i_pwr2)
687: an_tot2_d = an_tot2*an_t
689: d_eq(5,1) = &
690: & -(an_r(4)**i_cc)*(an_r(5)**i_h2o) &
691: & *((pt/an_t)**idiff) *(-idiff/an_t)
693: do jj = 2,26
694: d_eq(5,jj) = d_eq(5,1)
695: enddo
697: d_eq(5,1) = d_eq(5,1) + k_eq(11)* (an_r(3) **ai_o2)
699: d_eq(5,3) = d_eq(5,3) + k_eq(11)* (ai_o2*an_r(3)**(ai_o2-1)) &
700: & * an_r(1)
702: d_eq(5,4) = d_eq(5,4) - (i_cc*an_r(4)**(i_cc-1))* &
703: & (an_r(5)**i_h2o)* ((pt/an_t)**idiff)
704: d_eq(5,5) = d_eq(5,5) &
705: & - (i_h2o*(an_r(5)**(i_h2o-1))) &
706: & * (an_r(4)**i_cc)* ((pt/an_t)**idiff)
708: d_eq(3,1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
709: do jj = 2,26
710: d_eq(3,jj) = d_eq(3,1)
711: enddo
713: d_eq(3,2) = d_eq(3,2) + k_eq(12)* (an_r(3)**3)
715: d_eq(3,3) = d_eq(3,3) + k_eq(12)* (3*an_r(3)**2)*an_r(2)
717: d_eq(3,4) = d_eq(3,4) - 2*an_r(4)*(an_r(5)**3)*(pt/an_t)
719: d_eq(3,5) = d_eq(3,5) - 3*(an_r(5)**2)*(an_r(4)**2)*(pt/an_t)
720: ! & *(pt_five/const_five)
722: do ii3 = 1,26
723: d_eq(3,ii3) = 0.0d0
724: enddo
725: d_eq(3,2) = 1.0d0
727: d_eq(7,1) = pt*an_r(11)*(-1.0)/const2 &
728: & -k_eq(1)*sqrt(pt)*sqrt(an_r(14)+1d-50)*(-0.5/const3)
730: do jj = 2,26
731: d_eq(7,jj) = d_eq(7,1)
732: enddo
734: d_eq(7,11) = d_eq(7,11) + pt/an_t
735: d_eq(7,14) = d_eq(7,14) &
736: & - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14)+1d-50)*an_t)))
738: d_eq(8,1) = pt*an_r(8)*(-1.0)/const2 &
739: & -k_eq(2)*sqrt(pt)*sqrt(an_r(3)+1.0d-50)*(-0.5/const3)
741: do jj = 2,26
742: d_eq(8,jj) = d_eq(8,1)
743: enddo
745: d_eq(8,3) = d_eq(8,3) &
746: & -k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3)+1.0d-50)*an_t)))
747: d_eq(8,8) = d_eq(8,8) + pt/an_t
749: d_eq(9,1) = pt*an_r(7)*(-1.0)/const2 &
750: & -k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3)
752: do jj = 2,26
753: d_eq(9,jj) = d_eq(9,1)
754: enddo
756: d_eq(9,7) = d_eq(9,7) + pt/an_t
757: d_eq(9,6) = d_eq(9,6) &
758: & -k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))
760: d_eq(10,1) = pt*an_r(10)*(-1.0)/const2 &
761: & -k_eq(4)*(pt)*sqrt((an_r(3)+1.0d-50) &
762: & *an_r(14))*(-1.0/const2)
763: do jj = 2,26
764: d_eq(10,jj) = d_eq(10,1)
765: enddo
767: d_eq(10,3) = d_eq(10,3) &
768: & -k_eq(4)*(pt)*sqrt(an_r(14)) &
769: & *(0.5/(sqrt(an_r(3)+1.0d-50)*an_t))
770: d_eq(10,10) = d_eq(10,10) + pt/an_t
771: d_eq(10,14) = d_eq(10,14) &
772: & -k_eq(4)*(pt)*sqrt(an_r(3)+1.0d-50) &
773: & *(0.5/(sqrt(an_r(14)+1.0d-50)*an_t))
775: d_eq(11,1) = pt*an_r(9)*(-1.0)/const2 &
776: & -k_eq(5)*(pt)*sqrt((an_r(3)+1.0d-50)*an_r(6)) &
777: & *(-1.0/const2)
779: do jj = 2,26
780: d_eq(11,jj) = d_eq(11,1)
781: enddo
783: d_eq(11,3) = d_eq(11,3) &
784: & -k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ &
785: & (sqrt(an_r(3)+1.0d-50)*an_t))
786: d_eq(11,6) = d_eq(11,6) &
787: & -k_eq(5)*(pt)*sqrt(an_r(3)+1.0d-50) &
788: & *(0.5/(sqrt(an_r(6))*an_t))
789: d_eq(11,9) = d_eq(11,9) + pt/an_t
791: d_eq(12,1) = pt*an_r(5)*(-1.0)/const2 &
792: & -k_eq(6)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
793: & *(an_r(14))*(-1.5/const5)
795: do jj = 2,26
796: d_eq(12,jj) = d_eq(12,1)
797: enddo
799: d_eq(12,3) = d_eq(12,3) &
800: & -k_eq(6)*(pt**1.5)*((an_r(14)+1.0d-50)/const3) &
801: & *(0.5/sqrt(an_r(3)+1.0d-50))
803: d_eq(12,5) = d_eq(12,5) + pt/an_t
804: d_eq(12,14) = d_eq(12,14) &
805: & -k_eq(6)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
807: d_eq(13,1) = pt*an_r(4)*(-1.0)/const2 &
808: & -k_eq(7)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
809: & *(an_r(13))*(-1.5/const5)
811: do jj = 2,26
812: d_eq(13,jj) = d_eq(13,1)
813: enddo
815: d_eq(13,3) = d_eq(13,3) &
816: & -k_eq(7)*(pt**1.5)*(an_r(13)/const3) &
817: & *(0.5/sqrt(an_r(3)+1.0d-50))
819: d_eq(13,4) = d_eq(13,4) + pt/an_t
820: d_eq(13,13) = d_eq(13,13) &
821: & -k_eq(7)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
823: d_eq(14,1) = pt*an_r(15)*(-1.0)/const2 &
824: & -k_eq(8)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
825: & *(an_r(9))*(-1.5/const5)
827: do jj = 2,26
828: d_eq(14,jj) = d_eq(14,1)
829: enddo
831: d_eq(14,3) = d_eq(14,3) &
832: & -k_eq(8)*(pt**1.5)*(an_r(9)/const3) &
833: & *(0.5/sqrt(an_r(3)+1.0d-50))
834: d_eq(14,9) = d_eq(14,9) &
835: & -k_eq(8)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
836: d_eq(14,15) = d_eq(14,15)+ pt/an_t
838: d_eq(15,1) = pt*an_r(16)*(-1.0)/const2 &
839: & -k_eq(9)*(pt**1.5)*sqrt(an_r(14)+1.0d-50) &
840: & *(an_r(3))*(-1.5/const5)
842: do jj = 2,26
843: d_eq(15,jj) = d_eq(15,1)
844: enddo
846: d_eq(15,3) = d_eq(15,3) &
847: & -k_eq(9)*(pt**1.5)*(sqrt(an_r(14)+1.0d-50)/const3)
848: d_eq(15,14) = d_eq(15,14) &
849: & -k_eq(9)*(pt**1.5)*(an_r(3)/const3) &
850: & *(0.5/sqrt(an_r(14)+1.0d-50))
851: d_eq(15,16) = d_eq(15,16) + pt/an_t
853: d_eq(16,1) = pt*an_r(12)*(-1.0)/const2 &
854: & -k_eq(10)*(pt**1.5)*sqrt(an_r(3)+1.0d-50) &
855: & *(an_r(6))*(-1.5/const5)
857: do jj = 2,26
858: d_eq(16,jj) = d_eq(16,1)
859: enddo
861: d_eq(16,3) = d_eq(16,3) &
862: & -k_eq(10)*(pt**1.5)*(an_r(6)/const3) &
863: & *(0.5/sqrt(an_r(3)+1.0d-50))
865: d_eq(16,6) = d_eq(16,6) &
866: & -k_eq(10)*(pt**1.5)*(sqrt(an_r(3)+1.0d-50)/const3)
867: d_eq(16,12) = d_eq(16,12) + pt/an_t
869: const_cube = an_t*an_t*an_t
870: const_four = const2*const2
872: d_eq(17,1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
873: & - k_eq(15) * an_r(17)*pt * (-1/const2)
874: do jj = 2,26
875: d_eq(17,jj) = d_eq(17,1)
876: enddo
877: d_eq(17,14) = d_eq(17,14) + an_r(18)*an_r(18)*(pt**3)/const_cube
878: d_eq(17,17) = d_eq(17,17) - k_eq(15)*pt/an_t
879: d_eq(17,18) = d_eq(17,18) + 2*an_r(18)*an_r(14) &
880: & *(pt**3)/const_cube
882: d_eq(18,1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube) &
883: & - k_eq(16) * an_r(3)*an_r(18)*an_r(18) &
884: & * (pt*pt*pt) * (-3/const_four)
885: do jj = 2,26
886: d_eq(18,jj) = d_eq(18,1)
887: enddo
888: d_eq(18,3) = d_eq(18,3) &
889: & - k_eq(16) *an_r(18)* an_r(18)*pt*pt*pt /const_cube
890: d_eq(18,13) = d_eq(18,13) &
891: & + 2* an_r(13)*pt*pt /const2
892: d_eq(18,18) = d_eq(18,18) -k_eq(16)*an_r(3) &
893: & * 2*an_r(18)*pt*pt*pt/const_cube
895: !====for eq 19
897: d_eq(19,1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube) &
898: & - k_eq(17)*an_r(13)*an_r(10)*pt*pt * (-2/const_cube)
899: do jj = 2,26
900: d_eq(19,jj) = d_eq(19,1)
901: enddo
902: d_eq(19,13) = d_eq(19,13) &
903: & - k_eq(17) *an_r(10)*pt*pt /const2
904: d_eq(19,10) = d_eq(19,10) &
905: & - k_eq(17) *an_r(13)*pt*pt /const2
906: d_eq(19,3) = d_eq(19,3) + an_r(19)*pt*pt/const2
907: d_eq(19,19) = d_eq(19,19) + an_r(3)*pt*pt/const2
908: !====for eq 20
910: d_eq(20,1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube) &
911: & - k_eq(18) * an_r(19)*an_r(8)*pt*pt * (-2/const_cube)
912: do jj = 2,26
913: d_eq(20,jj) = d_eq(20,1)
914: enddo
915: d_eq(20,8) = d_eq(20,8) &
916: & - k_eq(18) *an_r(19)*pt*pt /const2
917: d_eq(20,19) = d_eq(20,19) &
918: & - k_eq(18) *an_r(8)*pt*pt /const2
919: d_eq(20,20) = d_eq(20,20) + an_r(21)*pt*pt/const2
920: d_eq(20,21) = d_eq(20,21) + an_r(20)*pt*pt/const2
922: !========
923: !====for eq 21
925: d_eq(21,1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube) &
926: & - k_eq(19)*an_r(7)*an_r(8)*pt*pt * (-2/const_cube)
927: do jj = 2,26
928: d_eq(21,jj) = d_eq(21,1)
929: enddo
930: d_eq(21,7) = d_eq(21,7) &
931: & - k_eq(19) *an_r(8)*pt*pt /const2
932: d_eq(21,8) = d_eq(21,8) &
933: & - k_eq(19) *an_r(7)*pt*pt /const2
934: d_eq(21,21) = d_eq(21,21) + an_r(23)*pt*pt/const2
935: d_eq(21,23) = d_eq(21,23) + an_r(21)*pt*pt/const2
937: !========
938: ! for 22
939: d_eq(22,1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube) &
940: & -k_eq(20)*an_r(21)*an_r(22)*pt*pt * (-2/const_cube)
941: do jj = 2,26
942: d_eq(22,jj) = d_eq(22,1)
943: enddo
944: d_eq(22,21) = d_eq(22,21) &
945: & - k_eq(20) *an_r(22)*pt*pt /const2
946: d_eq(22,22) = d_eq(22,22) &
947: & - k_eq(20) *an_r(21)*pt*pt /const2
948: d_eq(22,11) = d_eq(22,11) + an_r(5)*pt*pt/(const2)
949: d_eq(22,5) = d_eq(22,5) + an_r(11)*pt*pt/(const2)
951: !========
952: ! for 23
954: d_eq(23,1) = an_r(24)*(pt)*(-1/const2) &
955: & - k_eq(21)*an_r(21)*an_r(3)*pt*pt * (-2/const_cube)
956: do jj = 2,26
957: d_eq(23,jj) = d_eq(23,1)
958: enddo
959: d_eq(23,3) = d_eq(23,3) &
960: & - k_eq(21) *an_r(21)*pt*pt /const2
961: d_eq(23,21) = d_eq(23,21) &
962: & - k_eq(21) *an_r(3)*pt*pt /const2
963: d_eq(23,24) = d_eq(23,24) + pt/(an_t)
965: !========
966: ! for 24
967: d_eq(24,1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube) &
968: & - k_eq(22)*an_r(24)*an_r(8)*pt*pt * (-2/const_cube)
969: do jj = 2,26
970: d_eq(24,jj) = d_eq(24,1)
971: enddo
972: d_eq(24,8) = d_eq(24,8) &
973: & - k_eq(22) *an_r(24)*pt*pt /const2
974: d_eq(24,24) = d_eq(24,24) &
975: & - k_eq(22) *an_r(8)*pt*pt /const2
976: d_eq(24,3) = d_eq(24,3) + an_r(25)*pt*pt/const2
977: d_eq(24,25) = d_eq(24,25) + an_r(3)*pt*pt/const2
979: !========
980: !for 25
982: d_eq(25,1) = an_r(26)*(pt)*(-1/const2) &
983: & - k_eq(23)*an_r(21)*an_r(10)*pt*pt * (-2/const_cube)
984: do jj = 2,26
985: d_eq(25,jj) = d_eq(25,1)
986: enddo
987: d_eq(25,10) = d_eq(25,10) &
988: & - k_eq(23) *an_r(21)*pt*pt /const2
989: d_eq(25,21) = d_eq(25,21) &
990: & - k_eq(23) *an_r(10)*pt*pt /const2
991: d_eq(25,26) = d_eq(25,26) + pt/(an_t)
993: !============
994: ! for 26
995: d_eq(26,20) = -1
996: d_eq(26,22) = -1
997: d_eq(26,23) = -1
998: d_eq(26,21) = 1
999: d_eq(26,24) = 1
1000: d_eq(26,25) = 1
1001: d_eq(26,26) = 1
1003: do j = 1,26
1004: do i = 1,26
1005: write(44,*)i,j,d_eq(i,j)
1006: enddo
1007: enddo
1009: end
1011: subroutine ShashiPostCheck(ls,X,Y,W,c_Y,c_W,dummy)
1012: #include <petsc/finclude/petscsnes.h>
1013: use petscsnes
1014: implicit none
1015: SNESLineSearch ls
1016: PetscErrorCode ierr
1017: Vec X,Y,W
1018: PetscObject dummy
1019: PetscBool c_Y,c_W
1020: PetscScalar,pointer :: xx(:)
1021: PetscInt i
1022: PetscCall(VecGetArrayF90(W,xx,ierr))
1023: do i=1,26
1024: if (xx(i) < 0.0) then
1025: xx(i) = 0.0
1026: c_W = PETSC_TRUE
1027: endif
1028: if (xx(i) > 3.0) then
1029: xx(i) = 3.0
1030: endif
1031: enddo
1032: PetscCall(VecRestoreArrayF90(W,xx,ierr))
1033: end