c c This file is part of the Test Set for BVP solvers c http://www.dm.uniba.it/~bvpsolvers/ c c Problem SWIRLING FLOW III, bvpT33 c SPBVP of dimesion 6 c c DISCLAIMER: see c http://www.dm.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=34 c c The most recent version of this source file can be found at c http://www.dm.uniba.it/~bvpsolvers/bvpTestSet/fortransrc/problems/bvpT33.f c c c----------------------------------------------------------------------- c integer function pidate() pidate = 20130531 return end c----------------------------------------------------------------------- subroutine prob(fullnm,problm,type, + neqn,aleft,aright,nlbc,ms, + numjac,numbcjac,linear,rpar,ipar) implicit none character*(25), intent(out) :: fullnm character*(15), intent(out) ::problm character*(5), intent(out) :: type integer, intent(out) :: neqn,ipar(1), nlbc, ms(6) double precision, intent(out) :: aleft, aright, rpar(2) logical, intent(out) :: numjac, numbcjac,linear double precision :: EPS common /pars/ EPS c c For Singularly perturbed BVP (SPBVP) eps should be the c first elements of parms and the first element of rpar c fullnm = 'Problem Swiling Flow III' problm = 'bvpT33' type = 'SPBVP' ms(1:6) = 1 neqn = 6 nlbc = 3 aleft = 0.0d0 aright = 1.0d0 numjac = .false. numbcjac = .false. linear = .false. c c EPS is an input parameter in report.f c c rpar(1) may contains epsin, the starting value for the ontinuation codes c rpar(1) = 0 means that epsin is the default value rpar(1) = 0.5d0 c for this problem ipar is a dummy variable ipar(1) = 0 return end c c----------------------------------------------------------------- c subroutine settolerances(neqn,rtol,atol,tolvec,ntol,ltol) integer, intent(in) :: neqn integer, intent(out) :: ntol, ltol(*) logical, intent(out) :: tolvec double precision, intent(in out) :: rtol(neqn), atol(neqn) tolvec = .false. ntol = neqn ltol(1)=1 DO i=2,ntol ltol(i) = i rtol(i) = rtol(1) atol(i) = atol(1) ENDDO return end c--------------------------------------------------------------------- subroutine setoutput(neqn,solref,printsolout, + nindsol,indsol) logical, intent(out) :: solref, printsolout integer, intent(in) :: neqn integer, intent(out) :: nindsol integer, intent(out) :: indsol(*) c the reference solution is available solref = .false. c output file is required printsolout = .true. c default values if printsolout is .true. nindsol = neqn c only nindsol component of indsol are referenced do i=1,nindsol indsol(i) = i end do return end c----------------------------------------------------- c Initialiser for parameter common block c used only for the R interface c c SUBROUTINE initrparms(bvpparms) EXTERNAL bvpparms c c For Singularly Perturbed BVP (SPBVP) eps should be the c first elements of parms and the first element of rpar c DOUBLE PRECISION parms(1) COMMON / pars / parms CALL bvpparms(1, parms) RETURN END c----------------------------------------------------------------------- subroutine init(neqn,aleft,aright,ms,yval0, nmsh, givmsh, + givey,xguess,yguess,dmguess,nudim) integer, intent(in) :: neqn,nudim,ms(neqn) double precision, intent(in) :: aleft,aright integer, intent(out) :: nmsh double precision, intent(out) :: xguess(*),yguess(nudim,*) double precision, intent(out) :: dmguess(nudim,*) double precision, intent(out) :: yval0(1:neqn) logical, intent(out) :: givmsh, givey double precision h givmsh = .true. givey = .true. yval0(1:neqn) = 0.0d0 nmsh = 3*5+1 c---- dummy values, not used if givmsh=FALSE, givey=false h =(aright-aleft)/(nmsh-1) xguess(1) = aleft do i=2,nmsh xguess(i)=xguess(i-1)+h end do xguess(nmsh)= aright ij = 1 DO i=1,neqn do j=1,ms(i) yguess(ij,1:nmsh) = 0.0d0 ij = ij+1 end do end do DO i=1,neqn dmguess(i,1:nmsh) = 0.0d0 end do return end c----------------------------------------------------------------------- SUBROUTINE feval(NCOMP,X,Z,F,RPAR,IPAR) IMPLICIT NONE INTEGER NCOMP, IPAR , I DOUBLE PRECISION F, Z, RPAR, X DIMENSION Z(*),F(*) DIMENSION RPAR(*), IPAR(*) DOUBLE PRECISION EPS COMMON / pars / EPS F(1)=Z(2) F(2)=(Z(1)*Z(4) - Z(3)*Z(2))/EPS F(3)=Z(4) F(4)=Z(5) F(5)=Z(6) F(6)=(-Z(3)*Z(6) - Z(1)*Z(2))/EPS RETURN END subroutine bceval(i,neqn,z,G,rpar,ipar) integer, INTENT(IN) :: i, neqn, ipar(*) double precision, INTENT(IN) :: z(neqn), rpar(*) double precision, INTENT(OUT) :: G double precision :: EPS common /pars/ EPS IF (I.EQ.1) G=Z(1)+1.D0 IF (I.EQ.2) G=Z(3) IF (I.EQ.3) G=z(4) IF (I.EQ.4) G=Z(1)-1.D0 IF (I.EQ.5) G=Z(3) IF (I.EQ.6) G=Z(4) return end c c---------------------------------------------------------------------- c c SUBROUTINE jeval(NCOMP,X,Z,DF,RPAR,IPAR) IMPLICIT NONE INTEGER NCOMP, IPAR, I, J DOUBLE PRECISION X, Z, DF, RPAR DIMENSION Z(*),DF(NCOMP,*) DIMENSION RPAR(*), IPAR(*) DOUBLE PRECISION EPS COMMON / pars / EPS DF(1:ncomp,1:ncomp)=0.0d0 DF(1,2)=1.0D0 DF(2,1)=Z(4)/EPS DF(2,2)=-Z(3)/EPS DF(2,3)=-Z(2)/EPS DF(2,4)=Z(1)/EPS DF(3,4)=1.0D0 DF(4,5)=1.0D0 DF(5,6)=1.0D0 DF(6,1)=-Z(2)/EPS DF(6,2)=-Z(1)/EPS DF(6,3)=-Z(6)/EPS DF(6,6)=-Z(3)/EPS RETURN END c---------------------------------------------------------------------- c computation of the numerical jacobian for the boundary conditions c c SUBROUTINE dbceval(I,NCOMP,Z,DG,RPAR,IPAR) IMPLICIT NONE INTEGER I, NCOMP, IPAR DOUBLE PRECISION Z, DG, RPAR DIMENSION Z(*),DG(*) DIMENSION RPAR(*), IPAR(*) DOUBLE PRECISION EPS COMMON / pars / EPS DG(1)=0.D0 DG(2)=0.D0 DG(3)=0.D0 DG(4)=0.D0 DG(5)=0.D0 DG(6)=0.D0 IF (I.EQ.1.OR.I.EQ.4) DG(1)=1.D0 IF (I.EQ.2.OR.I.EQ.5) DG(3)=1.D0 IF (I.EQ.3.OR.I.EQ.6) DG(4)=1.D0 RETURN END c----------------------------------------------------------------------- subroutine solut(neqn,x,y,nmsh,rpar,ipar) integer, intent(in) :: neqn, ipar(*) double precision, intent(in) :: x(nmsh),rpar(*) double precision, intent(out) :: y(1:neqn,1:nmsh) double precision :: lambda, lambda2 common /pars/ lambda return end