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