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