1c4762a1bSJed Brown#include <petsc/finclude/petsc.h> 2*01fa2b5aSMartin Diehlmodule shashimodule 3e7a95102SMartin Diehl use petscsnes 4c4762a1bSJed Brown implicit none 5c4762a1bSJed Brown 6e7a95102SMartin 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 & 280e7a95102SMartin Diehl - (an_h(1)*an_r(1) + an_h_additive*an_r(2) & 281e7a95102SMartin Diehl + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) & 282e7a95102SMartin Diehl + an_r(16) + 2*an_r(17) + an_r(19) & 283e7a95102SMartin Diehl + an_r(20) + 3*an_r(22) + an_r(26)) 284c4762a1bSJed Brown 285c4762a1bSJed Brown f_eq(2) = atom_o_init & 286e7a95102SMartin Diehl - (an_o_additive*an_r(2) + 2*an_r(3) & 287e7a95102SMartin Diehl + 2*an_r(4) + an_r(5) & 288e7a95102SMartin Diehl + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) & 289e7a95102SMartin Diehl + 2*an_r(15) + 2*an_r(16) + an_r(20) + an_r(22) & 290e7a95102SMartin 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 & 295e7a95102SMartin Diehl - (an_c(1)*an_r(1) + an_c_additive*an_r(2) & 296e7a95102SMartin Diehl + an_r(4) + an_r(13) + 2*an_r(17) + an_r(18) & 297e7a95102SMartin 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 & 310e7a95102SMartin 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 & 314e7a95102SMartin Diehl - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) & 315e7a95102SMartin Diehl + an_r(15) & 316e7a95102SMartin Diehl + an_r(23)) 317c4762a1bSJed Brown 318c4762a1bSJed Brown f_eq(7) = part_p(11) & 319e7a95102SMartin Diehl - (k_eq(1)*sqrt(part_p(14) + 1d-23)) 320c4762a1bSJed Brown f_eq(8) = part_p(8) & 321e7a95102SMartin Diehl - (k_eq(2)*sqrt(part_p(3) + 1d-23)) 322c4762a1bSJed Brown 323c4762a1bSJed Brown f_eq(9) = part_p(7) & 324e7a95102SMartin Diehl - (k_eq(3)*sqrt(part_p(6) + 1d-23)) 325c4762a1bSJed Brown 326c4762a1bSJed Brown f_eq(10) = part_p(10) & 327e7a95102SMartin Diehl - (k_eq(4)*sqrt(part_p(3) + 1d-23)) & 328e7a95102SMartin Diehl *sqrt(part_p(14)) 329c4762a1bSJed Brown 330c4762a1bSJed Brown f_eq(11) = part_p(9) & 331e7a95102SMartin Diehl - (k_eq(5)*sqrt(part_p(3) + 1d-23)) & 332e7a95102SMartin Diehl *sqrt(part_p(6) + 1d-23) 333c4762a1bSJed Brown f_eq(12) = part_p(5) & 334e7a95102SMartin Diehl - (k_eq(6)*sqrt(part_p(3) + 1d-23)) & 335e7a95102SMartin Diehl *(part_p(14)) 336c4762a1bSJed Brown 337c4762a1bSJed Brown f_eq(13) = part_p(4) & 338e7a95102SMartin Diehl - (k_eq(7)*sqrt(part_p(3) + 1.0d-23)) & 339e7a95102SMartin Diehl *(part_p(13)) 340c4762a1bSJed Brown 341c4762a1bSJed Brown f_eq(14) = part_p(15) & 342e7a95102SMartin Diehl - (k_eq(8)*sqrt(part_p(3) + 1.0d-50)) & 343e7a95102SMartin Diehl *(part_p(9)) 344c4762a1bSJed Brown f_eq(15) = part_p(16) & 345e7a95102SMartin Diehl - (k_eq(9)*part_p(3)) & 346e7a95102SMartin Diehl *sqrt(part_p(14) + 1d-23) 347c4762a1bSJed Brown 348c4762a1bSJed Brown f_eq(16) = part_p(12) & 349e7a95102SMartin Diehl - (k_eq(10)*sqrt(part_p(3) + 1d-23)) & 350e7a95102SMartin Diehl *(part_p(6)) 351c4762a1bSJed Brown 352c4762a1bSJed Brown f_eq(17) = part_p(14)*part_p(18)**2 & 353e7a95102SMartin Diehl - (k_eq(15)*part_p(17)) 354c4762a1bSJed Brown 355c4762a1bSJed Brown f_eq(18) = (part_p(13)**2) & 356e7a95102SMartin 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)) & 374e7a95102SMartin 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) = & 576e7a95102SMartin Diehl -(an_r(4)**i_cc)*(an_r(5)**i_h2o) & 577e7a95102SMartin 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))* & 589e7a95102SMartin Diehl (an_r(5)**i_h2o)*((pt/an_t)**idiff) 590c4762a1bSJed Brown d_eq(5, 5) = d_eq(5, 5) & 591e7a95102SMartin Diehl - (i_h2o*(an_r(5)**(i_h2o - 1))) & 592e7a95102SMartin 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 & 614e7a95102SMartin 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) & 622e7a95102SMartin 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 & 625e7a95102SMartin 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) & 632e7a95102SMartin 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 & 636e7a95102SMartin 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) & 644e7a95102SMartin 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 & 647e7a95102SMartin Diehl - k_eq(4)*(pt)*sqrt((an_r(3) + 1.0d-50) & 648e7a95102SMartin 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) & 654e7a95102SMartin Diehl - k_eq(4)*(pt)*sqrt(an_r(14)) & 655e7a95102SMartin 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) & 658e7a95102SMartin Diehl - k_eq(4)*(pt)*sqrt(an_r(3) + 1.0d-50) & 659e7a95102SMartin 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 & 662e7a95102SMartin Diehl - k_eq(5)*(pt)*sqrt((an_r(3) + 1.0d-50)*an_r(6)) & 663e7a95102SMartin 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) & 670e7a95102SMartin Diehl - k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ & 671e7a95102SMartin Diehl (sqrt(an_r(3) + 1.0d-50)*an_t)) 672c4762a1bSJed Brown d_eq(11, 6) = d_eq(11, 6) & 673e7a95102SMartin Diehl - k_eq(5)*(pt)*sqrt(an_r(3) + 1.0d-50) & 674e7a95102SMartin 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 & 678e7a95102SMartin Diehl - k_eq(6)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & 679e7a95102SMartin 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) & 686e7a95102SMartin Diehl - k_eq(6)*(pt**1.5)*((an_r(14) + 1.0d-50)/const3) & 687e7a95102SMartin 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) & 691e7a95102SMartin 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 & 694e7a95102SMartin Diehl - k_eq(7)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & 695e7a95102SMartin 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) & 702e7a95102SMartin Diehl - k_eq(7)*(pt**1.5)*(an_r(13)/const3) & 703e7a95102SMartin 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) & 707e7a95102SMartin 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 & 710e7a95102SMartin Diehl - k_eq(8)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & 711e7a95102SMartin 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) & 718e7a95102SMartin Diehl - k_eq(8)*(pt**1.5)*(an_r(9)/const3) & 719e7a95102SMartin Diehl *(0.5/sqrt(an_r(3) + 1.0d-50)) 720c4762a1bSJed Brown d_eq(14, 9) = d_eq(14, 9) & 721e7a95102SMartin 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 & 725e7a95102SMartin Diehl - k_eq(9)*(pt**1.5)*sqrt(an_r(14) + 1.0d-50) & 726e7a95102SMartin 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) & 733e7a95102SMartin 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) & 735e7a95102SMartin Diehl - k_eq(9)*(pt**1.5)*(an_r(3)/const3) & 736e7a95102SMartin 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 & 740e7a95102SMartin Diehl - k_eq(10)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & 741e7a95102SMartin 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) & 748e7a95102SMartin Diehl - k_eq(10)*(pt**1.5)*(an_r(6)/const3) & 749e7a95102SMartin Diehl *(0.5/sqrt(an_r(3) + 1.0d-50)) 750c4762a1bSJed Brown 751c4762a1bSJed Brown d_eq(16, 6) = d_eq(16, 6) & 752e7a95102SMartin 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) & 759e7a95102SMartin 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) & 769e7a95102SMartin Diehl - k_eq(16)*an_r(3)*an_r(18)*an_r(18) & 770e7a95102SMartin 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) & 775e7a95102SMartin 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) & 777e7a95102SMartin 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) & 784e7a95102SMartin 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) & 789e7a95102SMartin Diehl - k_eq(17)*an_r(10)*pt*pt/const2 790c4762a1bSJed Brown d_eq(19, 10) = d_eq(19, 10) & 791e7a95102SMartin 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) & 797e7a95102SMartin 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) & 802e7a95102SMartin Diehl - k_eq(18)*an_r(19)*pt*pt/const2 803c4762a1bSJed Brown d_eq(20, 19) = d_eq(20, 19) & 804e7a95102SMartin 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) & 812e7a95102SMartin 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) & 817e7a95102SMartin Diehl - k_eq(19)*an_r(8)*pt*pt/const2 818c4762a1bSJed Brown d_eq(21, 8) = d_eq(21, 8) & 819e7a95102SMartin 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) & 826e7a95102SMartin 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) & 831e7a95102SMartin Diehl - k_eq(20)*an_r(22)*pt*pt/const2 832c4762a1bSJed Brown d_eq(22, 22) = d_eq(22, 22) & 833e7a95102SMartin 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) & 841e7a95102SMartin 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) & 846e7a95102SMartin Diehl - k_eq(21)*an_r(21)*pt*pt/const2 847c4762a1bSJed Brown d_eq(23, 21) = d_eq(23, 21) & 848e7a95102SMartin 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) & 854e7a95102SMartin 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) & 859e7a95102SMartin Diehl - k_eq(22)*an_r(24)*pt*pt/const2 860c4762a1bSJed Brown d_eq(24, 24) = d_eq(24, 24) & 861e7a95102SMartin 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) & 869e7a95102SMartin 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) & 874e7a95102SMartin Diehl - k_eq(23)*an_r(21)*pt*pt/const2 875c4762a1bSJed Brown d_eq(25, 21) = d_eq(25, 21) & 876e7a95102SMartin 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*01fa2b5aSMartin Diehlend module shashimodule 918e7a95102SMartin Diehl 919e7a95102SMartin Diehlprogram main 920e7a95102SMartin Diehl use petsc 921*01fa2b5aSMartin Diehl use shashimodule 922e7a95102SMartin Diehl implicit none 923e7a95102SMartin Diehl 924e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 925e7a95102SMartin Diehl! Variable declarations 926e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 927e7a95102SMartin Diehl! 928e7a95102SMartin Diehl! Variables: 929e7a95102SMartin Diehl! snes - nonlinear solver 930e7a95102SMartin Diehl! x, r - solution, residual vectors 931e7a95102SMartin Diehl! J - Jacobian matrix 932e7a95102SMartin Diehl! 933e7a95102SMartin Diehl SNES snes 934e7a95102SMartin Diehl Vec x, r, lb, ub 935e7a95102SMartin Diehl Mat J 936e7a95102SMartin Diehl PetscErrorCode ierr 937e7a95102SMartin Diehl PetscInt i2 938e7a95102SMartin Diehl PetscMPIInt size 939e7a95102SMartin Diehl PetscScalar, pointer :: xx(:) 940e7a95102SMartin Diehl PetscScalar zero, big 941e7a95102SMartin Diehl SNESLineSearch ls 942e7a95102SMartin Diehl! 943e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 944e7a95102SMartin Diehl! Beginning of program 945e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 946e7a95102SMartin Diehl 947e7a95102SMartin Diehl PetscCallA(PetscInitialize(ierr)) 948e7a95102SMartin Diehl PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)) 949e7a95102SMartin Diehl PetscCheckA(size == 1, PETSC_COMM_WORLD, 1, 'requires one process') 950e7a95102SMartin Diehl 951e7a95102SMartin Diehl big = 2.88 952e7a95102SMartin Diehl big = PETSC_INFINITY 953e7a95102SMartin Diehl zero = 0.0 954e7a95102SMartin Diehl i2 = 26 955e7a95102SMartin Diehl! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - 956e7a95102SMartin Diehl! Create nonlinear solver context 957e7a95102SMartin Diehl! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - 958e7a95102SMartin Diehl 959e7a95102SMartin Diehl PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr)) 960e7a95102SMartin Diehl 961e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 962e7a95102SMartin Diehl! Create matrix and vector data structures; set corresponding routines 963e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 964e7a95102SMartin Diehl 965e7a95102SMartin Diehl! Create vectors for solution and nonlinear function 966e7a95102SMartin Diehl 967e7a95102SMartin Diehl PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i2, x, ierr)) 968e7a95102SMartin Diehl PetscCallA(VecDuplicate(x, r, ierr)) 969e7a95102SMartin Diehl 970e7a95102SMartin Diehl! Create Jacobian matrix data structure 971e7a95102SMartin Diehl 972e7a95102SMartin Diehl PetscCallA(MatCreateDense(PETSC_COMM_SELF, 26, 26, 26, 26, PETSC_NULL_SCALAR, J, ierr)) 973e7a95102SMartin Diehl 974e7a95102SMartin Diehl! Set function evaluation routine and vector 975e7a95102SMartin Diehl 976e7a95102SMartin Diehl PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr)) 977e7a95102SMartin Diehl 978e7a95102SMartin Diehl! Set Jacobian matrix data structure and Jacobian evaluation routine 979e7a95102SMartin Diehl 980e7a95102SMartin Diehl PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr)) 981e7a95102SMartin Diehl 982e7a95102SMartin Diehl PetscCallA(VecDuplicate(x, lb, ierr)) 983e7a95102SMartin Diehl PetscCallA(VecDuplicate(x, ub, ierr)) 984e7a95102SMartin Diehl PetscCallA(VecSet(lb, zero, ierr)) 985e7a95102SMartin Diehl! PetscCallA(VecGetArray(lb,xx,ierr)) 986e7a95102SMartin Diehl! PetscCallA(ShashiLowerBound(xx)) 987e7a95102SMartin Diehl! PetscCallA(VecRestoreArray(lb,xx,ierr)) 988e7a95102SMartin Diehl PetscCallA(VecSet(ub, big, ierr)) 989e7a95102SMartin Diehl! PetscCallA(SNESVISetVariableBounds(snes,lb,ub,ierr)) 990e7a95102SMartin Diehl 991e7a95102SMartin Diehl PetscCallA(SNESGetLineSearch(snes, ls, ierr)) 992e7a95102SMartin Diehl PetscCallA(SNESLineSearchSetPostCheck(ls, ShashiPostCheck, 0, ierr)) 993e7a95102SMartin Diehl PetscCallA(SNESSetType(snes, SNESVINEWTONRSLS, ierr)) 994e7a95102SMartin Diehl 995e7a95102SMartin Diehl PetscCallA(SNESSetFromOptions(snes, ierr)) 996e7a95102SMartin Diehl 997e7a95102SMartin Diehl! set initial guess 998e7a95102SMartin Diehl 999e7a95102SMartin Diehl PetscCallA(VecGetArray(x, xx, ierr)) 1000e7a95102SMartin Diehl PetscCallA(ShashiInitialGuess(xx)) 1001e7a95102SMartin Diehl PetscCallA(VecRestoreArray(x, xx, ierr)) 1002e7a95102SMartin Diehl 1003e7a95102SMartin Diehl PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr)) 1004e7a95102SMartin Diehl 1005e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1006e7a95102SMartin Diehl! Free work space. All PETSc objects should be destroyed when they 1007e7a95102SMartin Diehl! are no longer needed. 1008e7a95102SMartin Diehl! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1009e7a95102SMartin Diehl PetscCallA(VecDestroy(lb, ierr)) 1010e7a95102SMartin Diehl PetscCallA(VecDestroy(ub, ierr)) 1011e7a95102SMartin Diehl PetscCallA(VecDestroy(x, ierr)) 1012e7a95102SMartin Diehl PetscCallA(VecDestroy(r, ierr)) 1013e7a95102SMartin Diehl PetscCallA(MatDestroy(J, ierr)) 1014e7a95102SMartin Diehl PetscCallA(SNESDestroy(snes, ierr)) 1015e7a95102SMartin Diehl PetscCallA(PetscFinalize(ierr)) 1016e7a95102SMartin Diehlend 1017