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