c ... 11/20/2001:: Sophia Lefantzi c The input might be on a patch by patch basis, but c the implementation of the Jacobian is for a single cell. c ... 07/20/2003:: Sophia Lefantzi c ... Extended the code to return the Jacobian for all the species too. c ... Now the user has two options: c ... 1. When the component is compiled without the ALL_SPECIES flag the call to c ... jacobian returns (ns x ns) array corresponding to the solution vector c ... [T, Y_1, Y_2, ..., Y_{N-1}]. c ... 2. When the component is compiled with the ALL_SPECIES flag the call to c ... jacobian returns ( (ns+1) x (ns+1) ) array corresponding to the solution vector c ... [T, Y_1, Y_2, ..., Y_{N}]. subroutine jacobian(pressure,n_slack_species_index, 1 nspecies,ncells,n_size,data_in,data_out, 2 ntmx,temptab,temp1tab,odtab,cgsptab, 3 hsptab,cpsptab,cgsptb_t,hsptb_t,cpsptb_t) IMPLICIT NONE include 'parameter.par' include 'chem.com' include 'chem_ref.com' integer ntmx real*8 temptab(ntmx) real*8 cgsptab(nsmx,ntmx) real*8 hsptab(nsmx,ntmx) real*8 cpsptab(nsmx,ntmx) real*8 cgsptb_t(nsmx,ntmx) real*8 hsptb_t(nsmx,ntmx) real*8 cpsptb_t(nsmx,ntmx) double precision temp1tab,odtab integer n_size,nspecies,ncells,n,j,k integer n_slack_species_index,kj,ki real*8 s0(nvmx),s(nvmx),capr0(nvmx),capr(nvmx),ds(nvmx) real*8 dpd(nvmx,nvmx) real*8 ajac0(nvmx,nvmx) c ... The solution vector in. If we are solving for (ns-1) species it is: c ... [Temp, Y(ns-1)], otherwise if we are solving for all the species (ns) it is: c ... [Temp, Y(ns)] real*8 data_in(n_size) real*8 data_out(n_size,n_size) real*8 pressure(ncells),temp(ncells) real*8 ys(nspecies,ncells) double precision epss,sum,err,dY c Crap here real*8 cgspl_t(nsmx) real*8 hspl_t(nsmx) real*8 cpspl_t(nsmx) real*8 ttab(50000) double precision t1,t2,dt,tv,tc,vtemp integer ng,nc,kk character*1024 fmt if (ns.ne. nspecies)stop 1 'species number wrong' n = ncells c--------------------------------------------------------------------------- c assign the input data do j=1,n ! j=cell id Sum = 0.0 temp(j) = data_in( (j-1) * (nspecies) + 1 ) #ifdef NON_DIM tc = temp(j) #else tc = temp(j) / tempref #endif s0(1) = temp(j) #ifdef ALL_SPECIES do k=1,ns s0(k+1) = data_in( (j-1)*(nspecies) + (1+k) ) enddo #else c Fill in the mass fractions of the species do k=1,ns-1 s0(k+1) = data_in( (j-1)*(nspecies) + (1+k) ) enddo #endif #ifdef NON_DIM p0 = pressure(j) #else p0 = pressure(j)/pref #endif end do call gjac(s0,ntmx,tempref,temp1tab,odtab,temptab,cgsptab,hsptab, 1 cpsptab,cgsptb_t,hsptb_t,cpsptb_t, 2 n_slack_species_index,n_size,data_out) return end c******************************************************************************* subroutine gjac(s,ntmx,tempref,temp1tab,odtab,temptab,cgsptab,hsptab, 1 cpsptab,cgsptb_t,hsptb_t,cpsptb_t, 2 n_slack_species_index,n_size,pdnc) IMPLICIT NONE include 'parameter.par' include 'chem.com' integer k,n_slack_species_index,kj,ki,ntmx,n_size double precision ys_sum,sum,wb,tc,rc real*8 s(nsmx),capr(nvmx),cspl(nsmx) real*8 ysl(nsmx) real*8 drda(nsmx+3,idim),pdnc(n_size,n_size) real*8 pd(nsmx+1,nsmx+1) double precision tempref,temp1tab,odtab real*8 temptab(ntmx) real*8 cgsptab(nsmx,ntmx) real*8 hsptab(nsmx,ntmx) real*8 cpsptab(nsmx,ntmx) real*8 cgsptb_t(nsmx,ntmx) real*8 hsptb_t(nsmx,ntmx) real*8 cpsptb_t(nsmx,ntmx) c------------------------------------------------------------------------------- ! ysl local ys vector ys_sum = 0.0d0 Sum = 0.0d0 #ifdef ALL_SPECIES do k=1,ns ysl(k) = s(k+1) end do #else c Compute the species to take the blame from the rest do k=2,ns Sum = Sum + s(k) end do ysl(n_slack_species_index) = 1.0 - Sum c Make sure that no negative mass fraction appears c if ( ysl(n_slack_species_index).lt.0.0 ) then c ysl(n_slack_species_index) = 0.0 c end if c Fill in the mass fractions of the rest of the species if (n_slack_species_index.eq.1) then do k=2,ns ysl(k) = s(k) end do elseif (n_slack_species_index.eq.ns) then do k=1,(ns-1) ysl(k) = s(k+1) end do else do k=1,(n_slack_species_index-1) ysl(k) = s(k+1) end do do k=(n_slack_species_index+1),ns ysl(k) = s(k) end do endif #endif ! wbar corresp to ysl sum = 0.0 do k=1,ns sum = sum + ysl(k)*ospwt(k) end do wb = 1.0d0/sum ! temperature #ifdef NON_DIM tc = s(1) #else tc = s(1)/tempref #endif ! density rc = p0 * wb / tc ! concentrations do k=1,ns cspl(k) = rc*ysl(k)*ospwt(k) end do call dojac(n_slack_species_index,tc,rc,ntmx,temptab,temp1tab, 1 odtab,cgsptab,hsptab,cpsptab,cgsptb_t,hsptb_t, 2 cpsptb_t,cspl,0,0,drda,pd,n_size,pdnc) return end c******************************************************************************* c ... get h for each of the pure species at given tc function gwtl(tc,temp1tab,odtab,temptab,ntmx,hsptab,wcheml) IMPLICIT NONE include 'parameter.par' include 'chem.com' double precision tc,sum,gwtl integer k,ntmx real*8 wcheml(nsmx),hspl(nsmx) double precision temp1tab,odtab real*8 temptab(nsmx) real*8 hsptab(nsmx,ntmx) c------------------------------------------------------------------------------- call gethspl(tc,temp1tab,odtab,temptab,ntmx,hsptab,hspl) sum = 0.0 do k=1,ns sum = sum + hspl(k) * wcheml(k) end do gwtl = - sum c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get h for each of the pure species at given tc subroutine gethspl(tc,temp1tab,odtab,temptab,ntmx,hsptab,hspl) IMPLICIT NONE include 'parameter.par' include 'chem.com' integer ntmx real*8 temptab(ntmx) real*8 hsptab(nsmx,ntmx) double precision temp1tab,odtab double precision jtab,tc,rtab,jtp,hsv integer k real*8 hspl(nsmx) c------------------------------------------------------------------------------- c ... linear interpolation from table of h_k(T) jtab = int((tc-temp1tab)*odtab) + 1 rtab = (tc-temptab(jtab))*odtab jtp = jtab + 1 do k=1,ns hsv = hsptab(k,jtab) hspl(k) = hsv + rtab*(hsptab(k,jtp)-hsv) end do c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get cp for each of the pure species at given tc function gcpl(tc,temp1tab,odtab,temptab,ntmx,cpsptab,ysl) IMPLICIT NONE include 'parameter.par' include 'chem.com' real*8 ysl(nsmx),cpspl(nsmx) double precision gcpl,tc,sum integer k, ntmx double precision temp1tab,odtab real*8 temptab(ntmx) real*8 cpsptab(nsmx,ntmx) c------------------------------------------------------------------------------- call getcpspl(tc,temp1tab,odtab,temptab,ntmx,cpsptab,cpspl) sum = 0.0 do k = 1,ns sum = sum + ysl(k)*cpspl(k) end do gcpl = sum c------------------------------------------------------------------------------- return end c******************************************************************************* c ... get cp for each of the pure species at given tc subroutine getcpspl(tc,temp1tab,odtab,temptab,ntmx,cpsptab, 1 cpspl) IMPLICIT NONE include 'parameter.par' include 'chem.com' integer ntmx real*8 temptab(ntmx) real*8 cpsptab(nsmx,ntmx) double precision temp1tab,odtab double precision tc,sum,jtab,rtab,jtp,cpv integer k real*8 cpspl(nsmx) c------------------------------------------------------------------------------- c ... linear interpolation from table of cp_k(T) jtab = int((tc-temp1tab)*odtab) + 1 rtab = (tc-temptab(jtab))*odtab jtp = jtab + 1 do k=1,ns cpv = cpsptab(k,jtab) cpspl(k) = cpv + rtab*(cpsptab(k,jtp)-cpv) end do c------------------------------------------------------------------------------- return end