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