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