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