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