c----------------------------------------------------------------------- c----------------------------------------------------------------------- c c This file is part of the Test Set for IVP solvers c http://www.dm.uniba.it/~testset/ c c Charge pump, DAE format c index 2 DAE of dimension 9 c c DISCLAIMER: see c http://www.dm.uniba.it/~testset/disclaimer.php c c The most recent version of this source file can be found at c http://www.dm.uniba.it/~testset/src/problems/pump.f c c This is revision c $Id: pump.F,v 1.3 2006/10/06 12:12:27 testset Exp $ c c----------------------------------------------------------------------- C C Copyright by Georg Denk, 1995 C C Reference: Michael Guenther, Georg Denk, Uwe Feldmann: C How models for MOS transistors reflect charge distribution C effects. C C Submitted to the journal Mathematical Modelling of Systems, C also: Preprint 1745, May 1995, Technische Hochschule Darmstadt, C Fachbereich Mathematik. C c----------------------------------------------------------------------- integer function pidate() pidate = 20060828 return end c----------------------------------------------------------------------- subroutine prob(fullnm,problm,type, + neqn,ndisc,t, + numjac,mljac,mujac, + nummas,mlmas,mumas, + ind) character*(*) fullnm, problm, type integer neqn,ndisc,mljac,mujac,mlmas,mumas,ind(*) double precision t(0:*) logical numjac, nummas integer i fullnm = 'Charge pump' problm = 'pump' type = 'DAE' neqn = 9 ndisc = 39 t(0) = 0d0 t(1) = 50d-9 t(2) = 60d-9 t(3) = 110d-9 do 10 i=4,39 t(i) = t(mod(i,4)) + dble(int(i/4))*120d-9 10 continue t(40) = 1200d-9 numjac = .true. mujac = neqn mljac = neqn mumas = 2 mlmas = 0 do 20 i=1,8 ind(i) = 1 20 continue ind(9) = 2 return end c----------------------------------------------------------------------- subroutine init(neqn,t,y,yprime,consis) integer neqn double precision t,y(neqn),yprime(neqn) logical consis double precision qgate,qsrc,qdrain y(1) = qgate (0d0,0d0,0d0) y(2) = 0d0 y(3) = qsrc (0d0,0d0,0d0) y(4) = 0d0 y(5) = qdrain (0d0,0d0,0d0) y(6) = 0d0 y(7) = 0d0 y(8) = 0d0 y(9) = 0d0 yprime(1) = 0d0 yprime(2) = 0d0 yprime(3) = 0d0 yprime(4) = 0d0 yprime(5) = 0d0 yprime(6) = 0d0 yprime(7) = 0d0 yprime(8) = 0d0 yprime(9) = 0d0 consis = .true. return end c----------------------------------------------------------------------- subroutine settolerances(neqn,rtol,atol,tolvec) integer neqn logical tolvec double precision rtol(neqn), atol(neqn) integer i do 10 i=9,6,-1 atol(i) = atol(1) rtol(i) = rtol(1) 10 continue c c charges are much smaller: c do 20 i=5,1,-1 atol(i) = 1d-6*atol(1) rtol(i) = rtol(1) 20 continue tolvec = .true. return end c----------------------------------------------------------------------- subroutine setoutput(neqn,solref,printsolout, + nindsol,indsol) logical solref, printsolout integer neqn, nindsol integer indsol(neqn) c the reference solution is available solref = .true. c output file is required printsolout = .true. c default values if printsolout is .true. nindsol = 4 c only nindsol component of indsol are referenced indsol(1) = 6 indsol(2) = 7 indsol(3) = 8 indsol(4) = 9 return end c----------------------------------------------------------------------- subroutine feval(neqn,t,y,yprime,f,ierr,rpar,ipar) integer neqn,ierr,ipar(*) double precision t,y(neqn),yprime(neqn),f(neqn),rpar(*) external vin double precision qgate, qsrc, qdrain, vin double precision capd, caps parameter (capd = 0.40d-12, caps = 1.60d-12) f(1) = -y(9) f(2) = 0d0 f(3) = 0d0 f(4) = -y(6) + vin(t) f(5) = y(1) - qgate (y(6), y(6)-y(7), y(6)-y(8)) f(6) = y(2) - caps*y(7) f(7) = y(3) - qsrc (y(6), y(6)-y(7), y(6)-y(8)) f(8) = y(4) - capd*y(8) f(9) = y(5) - qdrain(y(6), y(6)-y(7), y(6)-y(8)) return end c----------------------------------------------------------------------- subroutine jeval(ldim,neqn,t,y,yprime,dfdy,ierr,rpar,ipar) integer ldim,neqn,ierr,ipar(*) double precision t,y(neqn),yprime(neqn),dfdy(ldim,neqn),rpar(*) C C ... Dummy subroutine for the Jacobian C return end c----------------------------------------------------------------------- subroutine meval(ldim,neqn,t,y,yprime,dfddy,ierr,rpar,ipar) integer ldim,neqn,ierr,ipar(*) double precision t,y(neqn),yprime(neqn),dfddy(ldim,neqn),rpar(*) integer i,j do 20 j=1,neqn do 10 i=1,3 dfddy(i,j) = 0d0 10 continue 20 continue dfddy(3,1) = 1d0 dfddy(3,2) = 1d0 dfddy(2,3) = 1d0 dfddy(2,4) = 1d0 dfddy(1,5) = 1d0 return end c----------------------------------------------------------------------- subroutine solut(neqn,t,y) integer neqn double precision t,y(neqn) integer i c old reference solution computed with PSIDE, tol=2d-8 c c y(1) = 0.1262800429876759d-12 c do 10 i=2,8 c y(i) = 0d0 c 10 continue c y(9) = 0.1522565920949845d-03 c c computed using quadruple precision GAMD on an c Alphaserver DS20E, with a 667 MHz EV67 processor. c c rtol = atol = 1q-18 c h0 = 1.0q-37 c c y(1)= 0.126280042987675933170893D-12 c as for the CWI version we force to zero the c components 2 to 8 do 10 i=2,8 y(i) = 0d0 10 continue c y(2)= 0.675964116033495173748909D-40 c y(3)= -0.675964116033495173748909D-40 c y(4)= 0.675964116033495173748909D-40 c y(5)= -0.675964116033495173748909D-40 c y(6)= 0.371380333028958567423113D-30 c y(7)= 0.422477572520934483593064D-28 c y(8)= 0.168991029008373793437243D-27 y(9)= 0.152255686815577679043511D-03 return end c----------------------------------------------------------------------- c auxiliary functions required by subroutine feval c----------------------------------------------------------------------- C Copyright by Georg Denk, 1995 C C Reference: Michael Guenther, Georg Denk, Uwe Feldmann: C How models for MOS transistors reflect charge distribution C effects. C C Submitted to the journal Mathematical Modelling of Systems, C also: Preprint 1745, May 1995, Technische Hochschule Darmstadt, C Fachbereich Mathematik. C C------------------ Begincharge.f double precision function qgate (vgb,vgs,vgd) C C ... computes gate charge C double precision vgb, vgs, vgd C double precision vt0, gamma, phi, cox parameter (vt0 = 0.20d+0, gamma = 0.350d-1, phi = 1.010d+0, , cox = 4.0d-12) C double precision ugs, ugd, ugb, ubs, vfb, vte, ugst, ugdt C if ((vgs - vgd) .le. 0.d0) then ugs = vgd ugd = vgs else ugs = vgs ugd = vgd endif C ugb = vgb ubs = ugs - ugb C vfb = vt0 - gamma * dsqrt(phi) - phi vte = vt0 + gamma *( dsqrt(phi - ubs) - dsqrt(phi) ) C if ( ugb .le. vfb) then qgate = cox * (ugb - vfb) else if ( ugb .gt. vfb .and. ugs .le. vte ) then qgate = cox * gamma * * ( dsqrt((gamma/2.d0)**2.d0 + ugb - vfb) - * gamma/2.d0 ) else if ( ugb .gt .vfb .and. ugs .gt. vte ) then ugst = ugs - vte if (ugd.gt.vte) then ugdt = ugd - vte else ugdt = 0.d0 endif qgate = cox * ( (2.d0/3.d0) * (ugdt + ugst - * ((ugdt * ugst)/(ugdt+ugst)) ) + * gamma * dsqrt(phi-ubs) ) endif C return end C double precision function qbulk (vgb,vgs,vgd) C C ... computes bulk charge C double precision vgb, vgs, vgd C double precision vt0, gamma, phi, cox parameter (vt0 = 0.20d+0, gamma = 0.350d-1, phi = 1.010d+0, , cox = 4.0d-12) C double precision ugs, ugd, ugb, ubs, vfb, vte C if ((vgs - vgd) .le. 0.d0) then ugs = vgd ugd = vgs else ugs = vgs ugd = vgd endif C ugb = vgb ubs = ugs - ugb C vfb = vt0 - gamma * dsqrt(phi) - phi vte = vt0 + gamma *( dsqrt(phi - ubs) - dsqrt(phi) ) C if (ugb.le.vfb) then qbulk = - cox * (ugb - vfb) else if ((ugb.gt.vfb).and.(ugs.le.vte)) then qbulk = - cox * gamma * * ( dsqrt((gamma/2.d0)**2.d0 + ugb - vfb) - * gamma/2.d0 ) else if ( ugb .gt .vfb .and. ugs .gt. vte ) then qbulk = - cox * gamma * dsqrt(phi-ubs) endif C return end C double precision function qsrc (vgb,vgs,vgd) C C ... computes gate charge C double precision vgb, vgs, vgd C double precision vt0, gamma, phi, cox parameter (vt0 = 0.20d+0, gamma = 0.350d-1, phi = 1.010d+0, , cox = 4.0d-12) C double precision ugs, ugd, ugb, ubs, vfb, vte, ugst, ugdt C if ((vgs - vgd) .le. 0.d0) then ugs = vgd ugd = vgs else ugs = vgs ugd = vgd endif C ugb = vgb ubs = ugs - ugb C vfb = vt0 - gamma * dsqrt(phi) - phi vte = vt0 + gamma *( dsqrt(phi - ubs) - dsqrt(phi) ) C if (ugb.le.vfb) then qsrc = 0.d0 else if ((ugb.gt.vfb).and.(ugs.le.vte)) then qsrc = 0.d0 else if ( ugb .gt .vfb .and. ugs .gt. vte ) then ugst = ugs - vte if (ugd.ge.vte) then ugdt = ugd - vte else ugdt = 0.d0 endif qsrc = - cox * (1.d0/3.d0) * (ugdt + ugst - * ((ugdt * ugst)/(ugdt+ugst)) ) endif C return end C double precision function qdrain (vgb,vgs,vgd) C C ... computes drain charge C double precision vgb, vgs, vgd C double precision vt0, gamma, phi, cox parameter (vt0 = 0.20d+0, gamma = 0.350d-1, phi = 1.010d+0, , cox = 4.0d-12) C double precision ugs, ugd, ugb, ubs, vfb, vte, ugst, ugdt C if ((vgs - vgd) .le. 0.d0) then ugs = vgd ugd = vgs else ugs = vgs ugd = vgd endif C ugb = vgb ubs = ugs - ugb C vfb = vt0 - gamma * dsqrt(phi) - phi vte = vt0 + gamma *( dsqrt(phi - ubs) - dsqrt(phi) ) C if (ugb.le.vfb) then qdrain = 0.d0 else if ((ugb.gt.vfb).and.(ugs.le.vte)) then qdrain = 0.d0 else if ( ugb .gt .vfb .and. ugs .gt. vte ) then ugst = ugs - vte if (ugd.ge.vte) then ugdt = ugd - vte else ugdt = 0.d0 endif qdrain = - cox * (1.d0/3.d0) * (ugdt + ugst - * ((ugdt * ugst)/(ugdt+ugst)) ) endif C return end C double precision function vin(t) C C ... computes pulsed input voltage C double precision t C double precision vhigh, deltat, t1, t2, t3, dummy parameter (vhigh = 20.0d+0) C deltat = 120.0d-9 t1 = 50.0d-9 t2 = 60.0d-9 t3 = 110.0d-9 C dummy = dmod(t, deltat) if (dummy .lt. t1) then C ... 0 <= t' < 50ns vin = 0.0d+0 else if (dummy .lt. t2) then C ... 50ns <= t' < 60ns vin = (dummy - t1) * 0.10d+9 * vhigh else if (dummy .lt. t3) then C ... 60ns <= t' < 110ns vin = vhigh else C ... 110ns <= t' < 120ns vin = (deltat - dummy) * 0.10d+9 * vhigh endif endif endif C return end