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