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