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