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