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