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