Actual source code: shashi.F90
1: #include <petsc/finclude/petsc.h>
2: module shashimodule
3: use petscsnes
4: implicit none
5: PetscScalar, parameter :: &
6: an_c_additive = 2.0, &
7: an_h_additive = 6.0, &
8: an_o_additive = 1.0
9: PetscInt, parameter :: an_h = 1, an_c = 4
10: PetscScalar, parameter :: pt = 0.1
11: PetscScalar, parameter :: k_eq(23) = &
12: [1.75149D-05, &
13: 4.01405D-06, &
14: 6.04663D-14, &
15: 2.73612D-01, &
16: 3.25592D-03, &
17: 5.33568D+05, &
18: 2.07479D+05, &
19: 1.11841D-02, &
20: 1.72684D-03, &
21: 1.98588D-07, &
22: 7.23600D+27, &
23: 5.73926D+49, &
24: 1.00000D+00, &
25: 1.64493D+16, &
26: 2.73837D-29, &
27: 3.27419D+50, &
28: 1.72447D-23, &
29: 4.24657D-06, &
30: 1.16065D-14, &
31: 3.28020D+25, &
32: 1.06291D+00, &
33: 9.11507D+02, &
34: 6.02837D+03]
36: contains
37: !
38: ! ------------------------------------------------------------------------
39: !
40: ! FormFunction - Evaluates nonlinear function, F(x).
41: !
42: ! Input Parameters:
43: ! snes - the SNES context
44: ! x - input vector
45: ! dummy - optional user-defined context (not used here)
46: !
47: ! Output Parameter:
48: ! f - function vector
49: !
50: subroutine FormFunction(snes, x, f, dummy, ierr)
51: SNES snes
52: Vec x, f
53: PetscErrorCode, intent(out) :: ierr
54: integer dummy(*)
56: ! Declarations for use with local arrays
58: PetscScalar, pointer ::lx_v(:), lf_v(:)
60: ! Get pointers to vector data.
61: ! - For default PETSc vectors, VecGetArray() returns a pointer to
62: ! the data array. Otherwise, the routine is implementation dependent.
63: ! - You MUST call VecRestoreArray() when you no longer need access to
64: ! the array.
65: ! - Note that the Fortran interface to VecGetArray() differs from the
66: ! C version. See the Fortran chapter of the users manual for details.
68: PetscCall(VecGetArrayRead(x, lx_v, ierr))
69: PetscCall(VecGetArray(f, lf_v, ierr))
70: PetscCall(ShashiFormFunction(lx_v, lf_v))
71: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
72: PetscCall(VecRestoreArray(f, lf_v, ierr))
73: end
75: ! ---------------------------------------------------------------------
76: !
77: ! FormJacobian - Evaluates Jacobian matrix.
78: !
79: ! Input Parameters:
80: ! snes - the SNES context
81: ! x - input vector
82: ! dummy - optional user-defined context (not used here)
83: !
84: ! Output Parameters:
85: ! A - Jacobian matrix
86: ! B - optionally different matrix used to construct the preconditioner
87: !
88: subroutine FormJacobian(snes, X, jac, B, dummy, ierr)
89: SNES snes
90: Vec X
91: Mat jac, B
92: PetscErrorCode, intent(out) :: ierr
93: integer dummy(*)
95: ! Declarations for use with local arrays
96: PetscScalar, pointer ::lx_v(:), lf_v(:, :)
98: ! Get pointer to vector data
100: PetscCall(VecGetArrayRead(x, lx_v, ierr))
101: PetscCall(MatDenseGetArray(B, lf_v, ierr))
102: PetscCall(ShashiFormJacobian(lx_v, lf_v))
103: PetscCall(MatDenseRestoreArray(B, lf_v, ierr))
104: PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
106: ! Assemble matrix
108: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
109: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
111: end
113: subroutine ShashiLowerBound(an_r)
114: PetscScalar an_r(26)
115: PetscInt i
117: an_r(2:26) = 1000.0/6.023D+23
118: end
120: subroutine ShashiInitialGuess(an_r)
121: PetscScalar an_r(26)
123: an_r = 0.0
124: an_r(3) = .5
125: an_r(6) = 1.88000
126: an_r(14) = 1.
128: #if defined(solution)
129: an_r(1) = 3.802208D-33
130: an_r(2) = 1.298287D-29
131: an_r(3) = 2.533067D-04
132: an_r(4) = 6.865078D-22
133: an_r(5) = 9.993125D-01
134: an_r(6) = 1.879964D+00
135: an_r(7) = 4.449489D-13
136: an_r(8) = 3.428687D-07
137: an_r(9) = 7.105138D-05
138: an_r(10) = 1.094368D-04
139: an_r(11) = 2.362305D-06
140: an_r(12) = 1.107145D-09
141: an_r(13) = 1.276162D-24
142: an_r(14) = 6.315538D-04
143: an_r(15) = 2.356540D-09
144: an_r(16) = 2.048248D-09
145: an_r(17) = 1.966187D-22
146: an_r(18) = 7.856497D-29
147: an_r(19) = 1.987840D-36
148: an_r(20) = 8.182441D-22
149: an_r(21) = 2.684880D-16
150: an_r(22) = 2.680473D-16
151: an_r(23) = 6.594967D-18
152: an_r(24) = 2.509714D-21
153: an_r(25) = 3.096459D-21
154: an_r(26) = 6.149551D-18
155: #endif
156: end
158: subroutine ShashiFormFunction(an_r, f_eq)
159: PetscScalar, intent(in) :: an_r(26)
160: PetscScalar, intent(out) :: f_eq(26)
161: PetscScalar, parameter :: &
162: atom_c_init = 6.7408177364816552D-022, &
163: atom_h_init = 2.0, &
164: atom_o_init = 1.0, &
165: atom_n_init = 3.76
166: PetscScalar an_r(26), f_eq(26)
167: PetscScalar part_p(26), idiff
168: PetscInt i_cc, i_hh, i_h2o
169: PetscScalar an_t
170: PetscScalar a_io2
171: PetscInt i, ip
173: an_t = sum(an_r)
175: f_eq(1) = atom_h_init &
176: - (an_h*an_r(1) + an_h_additive*an_r(2) &
177: + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) &
178: + an_r(16) + 2*an_r(17) + an_r(19) &
179: + an_r(20) + 3*an_r(22) + an_r(26))
181: f_eq(2) = atom_o_init &
182: - (an_o_additive*an_r(2) + 2*an_r(3) &
183: + 2*an_r(4) + an_r(5) &
184: + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
185: + 2*an_r(15) + 2*an_r(16) + an_r(20) + an_r(22) &
186: + an_r(23) + 2*an_r(24) + 1*an_r(25) + an_r(26))
188: f_eq(3) = an_r(2) - 1.0d-150
190: f_eq(4) = atom_c_init &
191: - (an_c*an_r(1) + an_c_additive*an_r(2) &
192: + an_r(4) + an_r(13) + 2*an_r(17) + an_r(18) &
193: + an_r(19) + an_r(20))
195: part_p(1:26) = (an_r(1:26)/an_t)*pt
197: i_cc = an_c
198: i_hh = an_h
199: a_io2 = i_cc + i_hh/4.0
200: i_h2o = i_hh/2
201: idiff = (i_cc + i_h2o) - (a_io2 + 1)
203: f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2 &
204: - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
205: ! write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II'
206: ! stop
207: f_eq(6) = atom_n_init &
208: - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) &
209: + an_r(15) &
210: + an_r(23))
212: f_eq(7) = part_p(11) &
213: - (k_eq(1)*sqrt(part_p(14) + 1d-23))
214: f_eq(8) = part_p(8) &
215: - (k_eq(2)*sqrt(part_p(3) + 1d-23))
217: f_eq(9) = part_p(7) &
218: - (k_eq(3)*sqrt(part_p(6) + 1d-23))
220: f_eq(10) = part_p(10) &
221: - (k_eq(4)*sqrt(part_p(3) + 1d-23)) &
222: *sqrt(part_p(14))
224: f_eq(11) = part_p(9) &
225: - (k_eq(5)*sqrt(part_p(3) + 1d-23)) &
226: *sqrt(part_p(6) + 1d-23)
227: f_eq(12) = part_p(5) &
228: - (k_eq(6)*sqrt(part_p(3) + 1d-23)) &
229: *(part_p(14))
231: f_eq(13) = part_p(4) &
232: - (k_eq(7)*sqrt(part_p(3) + 1.0d-23)) &
233: *(part_p(13))
235: f_eq(14) = part_p(15) &
236: - (k_eq(8)*sqrt(part_p(3) + 1.0d-50)) &
237: *(part_p(9))
238: f_eq(15) = part_p(16) &
239: - (k_eq(9)*part_p(3)) &
240: *sqrt(part_p(14) + 1d-23)
242: f_eq(16) = part_p(12) &
243: - (k_eq(10)*sqrt(part_p(3) + 1d-23)) &
244: *(part_p(6))
246: f_eq(17) = part_p(14)*part_p(18)**2 &
247: - (k_eq(15)*part_p(17))
249: f_eq(18) = (part_p(13)**2) &
250: - (k_eq(16)*part_p(3)*part_p(18)**2)
251: print *, f_eq(18), part_p(3), part_p(18), part_p(13), k_eq(16)
253: f_eq(19) = part_p(19)*part_p(3) - k_eq(17)*part_p(13)*part_p(10)
255: f_eq(20) = part_p(21)*part_p(20) - k_eq(18)*part_p(19)*part_p(8)
257: f_eq(21) = part_p(21)*part_p(23) - k_eq(19)*part_p(7)*part_p(8)
259: f_eq(22) = part_p(5)*part_p(11) - k_eq(20)*part_p(21)*part_p(22)
261: f_eq(23) = part_p(24) - k_eq(21)*part_p(21)*part_p(3)
263: f_eq(24) = part_p(3)*part_p(25) - k_eq(22)*part_p(24)*part_p(8)
265: f_eq(25) = part_p(26) - k_eq(23)*part_p(21)*part_p(10)
267: f_eq(26) = -(an_r(20) + an_r(22) + an_r(23)) &
268: + (an_r(21) + an_r(24) + an_r(25) + an_r(26))
270: do i = 1, 26
271: write (44, *) i, f_eq(i) ! MarDiehl: What is "44"?
272: end do
274: end
276: subroutine ShashiFormJacobian(an_r, d_eq)
277: PetscScalar an_t, ai_o2
278: PetscScalar const5, const3
279: PetscInt jj, jb, ii3, id, ib, i
280: PetscScalar an_r(26)
281: PetscScalar d_eq(26, 26)
282: PetscScalar ai_pwr1, idiff
283: PetscInt i_cc, i_hh
284: PetscInt i_h2o, i_pwr2, i_co2_h2o
285: PetscInt j
287: d_eq = 0.0d0
288: an_t = sum(an_r)
290: d_eq(1, 1) = -an_h
291: d_eq(1, 2) = -an_h_additive
292: d_eq(1, 5) = -2
293: d_eq(1, 10) = -1
294: d_eq(1, 11) = -1
295: d_eq(1, 14) = -2
296: d_eq(1, 16) = -1
297: d_eq(1, 17) = -2
298: d_eq(1, 19) = -1
299: d_eq(1, 20) = -1
300: d_eq(1, 22) = -3
301: d_eq(1, 26) = -1
303: d_eq(2, 2) = -1*an_o_additive
304: d_eq(2, 3) = -2
305: d_eq(2, 4) = -2
306: d_eq(2, 5) = -1
307: d_eq(2, 8) = -1
308: d_eq(2, 9) = -1
309: d_eq(2, 10) = -1
310: d_eq(2, 12) = -1
311: d_eq(2, 13) = -1
312: d_eq(2, 15) = -2
313: d_eq(2, 16) = -2
314: d_eq(2, 20) = -1
315: d_eq(2, 22) = -1
316: d_eq(2, 23) = -1
317: d_eq(2, 24) = -2
318: d_eq(2, 25) = -1
319: d_eq(2, 26) = -1
321: d_eq(6, 6) = -2
322: d_eq(6, 7) = -1
323: d_eq(6, 9) = -1
324: d_eq(6, 12) = -2
325: d_eq(6, 15) = -1
326: d_eq(6, 23) = -1
328: d_eq(4, 1) = -an_c
329: d_eq(4, 2) = -an_c_additive
330: d_eq(4, 4) = -1
331: d_eq(4, 13) = -1
332: d_eq(4, 17) = -2
333: d_eq(4, 18) = -1
334: d_eq(4, 19) = -1
335: d_eq(4, 20) = -1
337: !----------
338: const3 = (an_t**2)*sqrt(an_t)
339: const5 = an_t*const3
341: i_cc = an_c
342: i_hh = an_h
343: ai_o2 = i_cc + float(i_hh)/4.0
344: i_co2_h2o = i_cc + i_hh/2
345: i_h2o = i_hh/2
346: ai_pwr1 = 1 + i_cc + float(i_hh)/4.0
347: i_pwr2 = i_cc + i_hh/2
348: idiff = (i_cc + i_h2o) - (ai_o2 + 1)
350: d_eq(5, :) = &
351: -(an_r(4)**i_cc)*(an_r(5)**i_h2o) &
352: *((pt/an_t)**idiff)*(-idiff/an_t)
353: d_eq(5, 1) = d_eq(5, 1) + k_eq(11)*(an_r(3)**ai_o2)
354: d_eq(5, 3) = d_eq(5, 3) + k_eq(11)*(ai_o2*an_r(3)**(ai_o2 - 1))*an_r(1)
355: d_eq(5, 4) = d_eq(5, 4) - (i_cc*an_r(4)**(i_cc - 1))* &
356: (an_r(5)**i_h2o)*((pt/an_t)**idiff)
357: d_eq(5, 5) = d_eq(5, 5) &
358: - (i_h2o*(an_r(5)**(i_h2o - 1))) &
359: *(an_r(4)**i_cc)*((pt/an_t)**idiff)
361: d_eq(3, :) = 0.0d0
362: d_eq(3, 2) = 1.0d0
364: d_eq(7, :) = pt*an_r(11)*(-1.0)/an_t**2 &
365: - k_eq(1)*sqrt(pt)*sqrt(an_r(14) + 1d-50)*(-0.5/const3)
366: d_eq(7, 11) = d_eq(7, 11) + pt/an_t
367: d_eq(7, 14) = d_eq(7, 14) &
368: - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14) + 1d-50)*an_t)))
370: d_eq(8, :) = pt*an_r(8)*(-1.0)/an_t**2 &
371: - k_eq(2)*sqrt(pt)*sqrt(an_r(3) + 1.0d-50)*(-0.5/const3)
372: d_eq(8, 3) = d_eq(8, 3) &
373: - k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3) + 1.0d-50)*an_t)))
374: d_eq(8, 8) = d_eq(8, 8) + pt/an_t
376: d_eq(9, :) = pt*an_r(7)*(-1.0)/an_t**2 &
377: - k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3)
378: d_eq(9, 7) = d_eq(9, 7) + pt/an_t
379: d_eq(9, 6) = d_eq(9, 6) &
380: - k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))
382: d_eq(10, :) = pt*an_r(10)*(-1.0)/an_t**2 &
383: - k_eq(4)*(pt)*sqrt((an_r(3) + 1.0d-50) &
384: *an_r(14))*(-1.0/an_t**2)
385: d_eq(10, 3) = d_eq(10, 3) &
386: - k_eq(4)*(pt)*sqrt(an_r(14)) &
387: *(0.5/(sqrt(an_r(3) + 1.0d-50)*an_t))
388: d_eq(10, 10) = d_eq(10, 10) + pt/an_t
389: d_eq(10, 14) = d_eq(10, 14) &
390: - k_eq(4)*(pt)*sqrt(an_r(3) + 1.0d-50) &
391: *(0.5/(sqrt(an_r(14) + 1.0d-50)*an_t))
393: d_eq(11, :) = pt*an_r(9)*(-1.0)/an_t**2 &
394: - k_eq(5)*(pt)*sqrt((an_r(3) + 1.0d-50)*an_r(6)) &
395: *(-1.0/an_t**2)
396: d_eq(11, 3) = d_eq(11, 3) &
397: - k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ &
398: (sqrt(an_r(3) + 1.0d-50)*an_t))
399: d_eq(11, 6) = d_eq(11, 6) &
400: - k_eq(5)*(pt)*sqrt(an_r(3) + 1.0d-50) &
401: *(0.5/(sqrt(an_r(6))*an_t))
402: d_eq(11, 9) = d_eq(11, 9) + pt/an_t
404: d_eq(12, :) = pt*an_r(5)*(-1.0)/an_t**2 &
405: - k_eq(6)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
406: *(an_r(14))*(-1.5/const5)
407: d_eq(12, 3) = d_eq(12, 3) &
408: - k_eq(6)*(pt**1.5)*((an_r(14) + 1.0d-50)/const3) &
409: *(0.5/sqrt(an_r(3) + 1.0d-50))
410: d_eq(12, 5) = d_eq(12, 5) + pt/an_t
411: d_eq(12, 14) = d_eq(12, 14) &
412: - k_eq(6)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
414: d_eq(13, :) = pt*an_r(4)*(-1.0)/an_t**2 &
415: - k_eq(7)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
416: *(an_r(13))*(-1.5/const5)
417: d_eq(13, 3) = d_eq(13, 3) &
418: - k_eq(7)*(pt**1.5)*(an_r(13)/const3) &
419: *(0.5/sqrt(an_r(3) + 1.0d-50))
420: d_eq(13, 4) = d_eq(13, 4) + pt/an_t
421: d_eq(13, 13) = d_eq(13, 13) &
422: - k_eq(7)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
424: d_eq(14, :) = pt*an_r(15)*(-1.0)/an_t**2 &
425: - k_eq(8)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
426: *(an_r(9))*(-1.5/const5)
427: d_eq(14, 3) = d_eq(14, 3) &
428: - k_eq(8)*(pt**1.5)*(an_r(9)/const3) &
429: *(0.5/sqrt(an_r(3) + 1.0d-50))
430: d_eq(14, 9) = d_eq(14, 9) &
431: - k_eq(8)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
432: d_eq(14, 15) = d_eq(14, 15) + pt/an_t
434: d_eq(15, :) = pt*an_r(16)*(-1.0)/an_t**2 &
435: - k_eq(9)*(pt**1.5)*sqrt(an_r(14) + 1.0d-50) &
436: *(an_r(3))*(-1.5/const5)
437: d_eq(15, 3) = d_eq(15, 3) &
438: - k_eq(9)*(pt**1.5)*(sqrt(an_r(14) + 1.0d-50)/const3)
439: d_eq(15, 14) = d_eq(15, 14) &
440: - k_eq(9)*(pt**1.5)*(an_r(3)/const3) &
441: *(0.5/sqrt(an_r(14) + 1.0d-50))
442: d_eq(15, 16) = d_eq(15, 16) + pt/an_t
444: d_eq(16, :) = pt*an_r(12)*(-1.0)/an_t**2 &
445: - k_eq(10)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
446: *(an_r(6))*(-1.5/const5)
447: d_eq(16, 3) = d_eq(16, 3) &
448: - k_eq(10)*(pt**1.5)*(an_r(6)/const3) &
449: *(0.5/sqrt(an_r(3) + 1.0d-50))
450: d_eq(16, 6) = d_eq(16, 6) &
451: - k_eq(10)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
452: d_eq(16, 12) = d_eq(16, 12) + pt/an_t
454: d_eq(17, :) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/an_t**4) &
455: - k_eq(15)*an_r(17)*pt*(-1/an_t**2)
456: d_eq(17, 14) = d_eq(17, 14) + an_r(18)*an_r(18)*(pt**3)/an_t**3
457: d_eq(17, 17) = d_eq(17, 17) - k_eq(15)*pt/an_t
458: d_eq(17, 18) = d_eq(17, 18) + 2*an_r(18)*an_r(14)*(pt**3)/an_t**3
460: d_eq(18, :) = an_r(13)*an_r(13)*(pt**2)*(-2/an_t**3) &
461: - k_eq(16)*an_r(3)*an_r(18)*an_r(18)*(pt**3)*(-3/an_t**4)
462: d_eq(18, 3) = d_eq(18, 3) &
463: - k_eq(16)*an_r(18)*an_r(18)*pt**3/an_t**3
464: d_eq(18, 13) = d_eq(18, 13) &
465: + 2*an_r(13)*pt**2/an_t**2
466: d_eq(18, 18) = d_eq(18, 18) - k_eq(16)*an_r(3)*2*an_r(18)*pt**3/an_t**3
468: d_eq(19, :) = an_r(3)*an_r(19)*(pt**2)*(-2/an_t**3) &
469: - k_eq(17)*an_r(13)*an_r(10)*pt**2*(-2/an_t**3)
470: d_eq(19, 13) = d_eq(19, 13) &
471: - k_eq(17)*an_r(10)*pt**2/an_t**2
472: d_eq(19, 10) = d_eq(19, 10) &
473: - k_eq(17)*an_r(13)*pt**2/an_t**2
474: d_eq(19, 3) = d_eq(19, 3) + an_r(19)*pt**2/an_t**2
475: d_eq(19, 19) = d_eq(19, 19) + an_r(3)*pt**2/an_t**2
477: d_eq(20, :) = an_r(21)*an_r(20)*(pt**2)*(-2/an_t**3) &
478: - k_eq(18)*an_r(19)*an_r(8)*pt**2*(-2/an_t**3)
479: d_eq(20, 8) = d_eq(20, 8) &
480: - k_eq(18)*an_r(19)*pt**2/an_t**2
481: d_eq(20, 19) = d_eq(20, 19) &
482: - k_eq(18)*an_r(8)*pt**2/an_t**2
483: d_eq(20, 20) = d_eq(20, 20) + an_r(21)*pt**2/an_t**2
484: d_eq(20, 21) = d_eq(20, 21) + an_r(20)*pt**2/an_t**2
486: d_eq(21, :) = an_r(21)*an_r(23)*(pt**2)*(-2/an_t**3) &
487: - k_eq(19)*an_r(7)*an_r(8)*pt**2*(-2/an_t**3)
488: d_eq(21, 7) = d_eq(21, 7) &
489: - k_eq(19)*an_r(8)*pt**2/an_t**2
490: d_eq(21, 8) = d_eq(21, 8) &
491: - k_eq(19)*an_r(7)*pt**2/an_t**2
492: d_eq(21, 21) = d_eq(21, 21) + an_r(23)*pt**2/an_t**2
493: d_eq(21, 23) = d_eq(21, 23) + an_r(21)*pt**2/an_t**2
495: d_eq(22, :) = an_r(5)*an_r(11)*(pt**2)*(-2/an_t**3) &
496: - k_eq(20)*an_r(21)*an_r(22)*pt**2*(-2/an_t**3)
497: d_eq(22, 21) = d_eq(22, 21) &
498: - k_eq(20)*an_r(22)*pt**2/an_t**2
499: d_eq(22, 22) = d_eq(22, 22) &
500: - k_eq(20)*an_r(21)*pt**2/an_t**2
501: d_eq(22, 11) = d_eq(22, 11) + an_r(5)*pt**2/(an_t**2)
502: d_eq(22, 5) = d_eq(22, 5) + an_r(11)*pt**2/(an_t**2)
504: d_eq(23, :) = an_r(24)*(pt)*(-1/an_t**2) &
505: - k_eq(21)*an_r(21)*an_r(3)*pt**2*(-2/an_t**3)
506: d_eq(23, 3) = d_eq(23, 3) &
507: - k_eq(21)*an_r(21)*pt**2/an_t**2
508: d_eq(23, 21) = d_eq(23, 21) &
509: - k_eq(21)*an_r(3)*pt**2/an_t**2
510: d_eq(23, 24) = d_eq(23, 24) + pt/(an_t)
512: d_eq(24, :) = an_r(3)*an_r(25)*(pt**2)*(-2/an_t**3) &
513: - k_eq(22)*an_r(24)*an_r(8)*pt**2*(-2/an_t**3)
514: d_eq(24, 8) = d_eq(24, 8) &
515: - k_eq(22)*an_r(24)*pt**2/an_t**2
516: d_eq(24, 24) = d_eq(24, 24) &
517: - k_eq(22)*an_r(8)*pt**2/an_t**2
518: d_eq(24, 3) = d_eq(24, 3) + an_r(25)*pt**2/an_t**2
519: d_eq(24, 25) = d_eq(24, 25) + an_r(3)*pt**2/an_t**2
521: d_eq(25, :) = an_r(26)*(pt)*(-1/an_t**2) &
522: - k_eq(23)*an_r(21)*an_r(10)*pt**2*(-2/an_t**3)
523: d_eq(25, 10) = d_eq(25, 10) &
524: - k_eq(23)*an_r(21)*pt**2/an_t**2
525: d_eq(25, 21) = d_eq(25, 21) &
526: - k_eq(23)*an_r(10)*pt**2/an_t**2
527: d_eq(25, 26) = d_eq(25, 26) + pt/(an_t)
529: d_eq(26, 20) = -1
530: d_eq(26, 22) = -1
531: d_eq(26, 23) = -1
532: d_eq(26, 21) = 1
533: d_eq(26, 24) = 1
534: d_eq(26, 25) = 1
535: d_eq(26, 26) = 1
537: do j = 1, 26
538: do i = 1, 26
539: write (44, *) i, j, d_eq(i, j)
540: end do
541: end do
543: end
545: subroutine ShashiPostCheck(ls, X, Y, W, c_Y, c_W, dummy)
546: SNESLineSearch ls
547: PetscErrorCode ierr
548: Vec X, Y, W
549: PetscObject dummy
550: PetscBool c_Y, c_W
551: PetscScalar, pointer :: xx(:)
552: PetscInt i
553: PetscCall(VecGetArray(W, xx, ierr))
554: do i = 1, 26
555: if (xx(i) < 0.0) then
556: xx(i) = 0.0
557: c_W = PETSC_TRUE
558: end if
559: if (xx(i) > 3.0) xx(i) = 3.0
560: end do
561: PetscCall(VecRestoreArray(W, xx, ierr))
562: end
563: end module shashimodule
565: program main
566: use petsc
567: use shashimodule
568: implicit none
570: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
571: ! Variable declarations
572: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
573: !
574: ! Variables:
575: ! snes - nonlinear solver
576: ! x, r - solution, residual vectors
577: ! J - Jacobian matrix
578: !
579: SNES snes
580: Vec x, r, lb, ub
581: Mat J
582: PetscErrorCode ierr
583: PetscMPIInt size
584: PetscScalar, pointer :: xx(:)
585: PetscScalar, parameter :: zero = 0.0
586: PetscScalar big
587: SNESLineSearch ls
588: !
589: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
590: ! Beginning of program
591: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
593: PetscCallA(PetscInitialize(ierr))
594: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
595: PetscCheckA(size == 1, PETSC_COMM_WORLD, 1, 'requires one process')
597: big = PETSC_INFINITY
598: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
599: ! Create nonlinear solver context
600: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
602: PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
604: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
605: ! Create matrix and vector data structures; set corresponding routines
606: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
608: ! Create vectors for solution and nonlinear function
609: PetscCallA(VecCreateSeq(PETSC_COMM_SELF, 26_PETSC_INT_KND, x, ierr))
610: PetscCallA(VecDuplicate(x, r, ierr))
612: ! Create Jacobian matrix data structure
613: PetscCallA(MatCreateDense(PETSC_COMM_SELF, 26, 26, 26, 26, PETSC_NULL_SCALAR, J, ierr))
615: ! Set function evaluation routine and vector
616: PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr))
618: ! Set Jacobian matrix data structure and Jacobian evaluation routine
619: PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr))
621: PetscCallA(VecDuplicate(x, lb, ierr))
622: PetscCallA(VecDuplicate(x, ub, ierr))
623: PetscCallA(VecSet(lb, zero, ierr))
624: ! PetscCallA(VecGetArray(lb,xx,ierr))
625: ! PetscCallA(ShashiLowerBound(xx))
626: ! PetscCallA(VecRestoreArray(lb,xx,ierr))
627: PetscCallA(VecSet(ub, big, ierr))
628: ! PetscCallA(SNESVISetVariableBounds(snes,lb,ub,ierr))
630: PetscCallA(SNESGetLineSearch(snes, ls, ierr))
631: PetscCallA(SNESLineSearchSetPostCheck(ls, ShashiPostCheck, 0, ierr))
632: PetscCallA(SNESSetType(snes, SNESVINEWTONRSLS, ierr))
634: PetscCallA(SNESSetFromOptions(snes, ierr))
636: ! set initial guess
637: PetscCallA(VecGetArray(x, xx, ierr))
638: PetscCallA(ShashiInitialGuess(xx))
639: PetscCallA(VecRestoreArray(x, xx, ierr))
641: PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
643: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
644: ! Free work space. All PETSc objects should be destroyed when they
645: ! are no longer needed.
646: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
647: PetscCallA(VecDestroy(lb, ierr))
648: PetscCallA(VecDestroy(ub, ierr))
649: PetscCallA(VecDestroy(x, ierr))
650: PetscCallA(VecDestroy(r, ierr))
651: PetscCallA(MatDestroy(J, ierr))
652: PetscCallA(SNESDestroy(snes, ierr))
653: PetscCallA(PetscFinalize(ierr))
654: end