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