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