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