c c This file is part of the Test Set for BVP solvers c http://www.dm.uniba.it/~bvpsolvers/ c c generic ACDCC.f driver c 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/drivers/acdccd.f c c c----------------------------------------------------------------------- c program acdccd IMPLICIT NONE INTEGER LWRKFL, LWRKIN, NUDIM, NXXDIM, LISERIES, IPRINT, IND, + NTOL, NLBC, NMSH, NFXPNT, NMAX, + IFLBVP, INDNMS, K, NMINIT, IDUM, + NUGDIM DOUBLE PRECISION ALEFT, ARIGHT, CKAPPA1, GAMMA1, + CKAPPA,CKAPPA2,SIGMA, EPS, EPSIN, EPSMIN PARAMETER(NUDIM=20,NXXDIM=10000) PARAMETER(LWRKFL=21*NUDIM+16*NUDIM*NUDIM+14*NUDIM*NXXDIM * +5*NUDIM*NUDIM*NXXDIM+6*NXXDIM) PARAMETER(LWRKIN=2*NUDIM*NXXDIM+2*NXXDIM+3*NUDIM) PARAMETER(NUGDIM=NUDIM) DOUBLE PRECISION WRK(LWRKFL) INTEGER IWRK(LWRKIN) DOUBLE PRECISION U(NUDIM,NXXDIM), XX(NXXDIM) DOUBLE PRECISION YREF(NUDIM,NXXDIM) DOUBLE PRECISION XGUESS(NXXDIM),YGUESS(NUDIM,NXXDIM) DOUBLE PRECISION DMGUESS(NUDIM,NXXDIM) DOUBLE PRECISION UACC(NUDIM,NXXDIM) DOUBLE PRECISION TOL(NUDIM), FIXPNT(1) INTEGER MS(NUDIM) DOUBLE PRECISION RPAR(10) DOUBLE PRECISION ATOL(NUDIM) double precision yval0(1:nudim) INTEGER LTOL(NUDIM),IPAR(10) PARAMETER (LISERIES=100) DOUBLE PRECISION iseries(liseries) INTEGER nmguess INTEGER iset(8) DOUBLE PRECISION precis(3) integer n0 LOGICAL LINEAR, GIVEU, GIVMSH,giv_u integer ureset logical useC logical full, ACCMESHFAIL LOGICAL use_c, comp_c, tolvec EXTERNAL ACDCC EXTERNAL odef EXTERNAL odejac EXTERNAL gsub EXTERNAL dgsub EXTERNAL fnumjac EXTERNAL bcnumjac c this common need to be defined in order to run TWPBVPC successfully c give information about some common parameters common/algprs/ nminit, iprint, idum, use_c, comp_c c Common/acAlgprs/Maxcon,Itsaim,Uval0 common/acgu/ giv_u, ureset character fullnm*25, problm*15, bvptype*5 character driver*9, solver*8 parameter (driver = 'acdccd', solver='ACDCC') real gettim, timtim, cputim external gettim double precision scd, mescd character fileout*140, filepath*100 character namefile*100 logical printsolout, solref, numjac, numbcjac integer nindsol, indsol(nudim),neqn, mstar, i logical stab_kappa, GIVEPS common /pars/ EPS c----------------------------------------------------------------------- c check the problem definition interface date c----------------------------------------------------------------------- call chkdat(driver,solver,20130531) c----------------------------------------------------------------------- c get the problem dependent parameters c----------------------------------------------------------------------- fullnm(1:20)=' ' fileout(1:140)=' ' filepath(1:100)=' ' namefile(1:100)=' ' problm(1:12)=' ' bvptype(1:5)=' ' rpar(1:10) = 0.0d0 ipar(1:10) = 0 call prob(fullnm,problm,bvptype, + neqn,aleft,aright,nlbc,ms, + numjac,numbcjac,linear,rpar,ipar) mstar=0 do i=1,neqn mstar=mstar+ms(i) end do if ( mstar .ne. neqn )then print *, 'ACDCCD: ERROR: ', + 'ACDCC cannot solve higher order problems' stop endif if ( bvptype(1:lnblnk(bvptype)) .ne. 'BVP' .and. + bvptype(1:lnblnk(bvptype)) .ne. 'SPBVP' )then print *, 'ACDCCD: ERROR: ', + 'unknown Test Set problem bvptype ', bvptype stop endif if ( bvptype(1:lnblnk(bvptype)) .ne. 'SPBVP' ) then print *, 'ACDCCD: ERROR: ', + 'ACDCC cannot solve ', bvptype, ' problems ' stop endif c----------------------------------------------------------------------- c read the tolerances and initial stepsize c----------------------------------------------------------------------- call getinp(driver,problm,solver,fullnm, + tolvec,tol,atol,bvptype) c----------------------------------------------------------------------- c get the initial values c----------------------------------------------------------------------- call init(neqn,aleft,aright,ms,yval0,nmsh,givmsh,giveu,xguess, + yguess,dmguess,nudim) n0 = nmsh giv_u=giveu nmguess = nmsh XX(1:NMSH)=XGUESS(1:NMSH) XX(1) = aleft XX(NMSH) = aright U(1:neqn,1:nmsh) = yguess(1:neqn,1:nmsh) call settolerances(neqn,tol,atol,tolvec,ntol,ltol) call setoutput(neqn,solref,printsolout, + nindsol,indsol) call getout(printsolout) EPSMIN = EPS EPS = RPAR(1) IF (EPS .EQ. 0.0d0) THEN GIVEPS = .FALSE. ELSE GIVEPS = .TRUE. END IF USE_C = .true. COMP_C = .true. precis(3) = epsilon(1.0d0) precis(2) = huge(1.0d0) precis(1) = tiny(1.0d0) IPRINT = -1 DO ind=1,LWRKFL WRK(ind) = 0d0 ENDDO DO ind=1,LWRKIN IWRK(ind) = 0 ENDDO NFXPNT= 0 indnms=0 full=.false. useC=use_c if (printsolout) then call getarg(1,filepath) call getarg(2,namefile) if (lnblnk(namefile) .gt. 0) then write(fileout,'(a,a,a,a,a)') filepath(1:lnblnk(filepath)), + namefile(1:lnblnk(namefile)),'_', + solver(1:lnblnk(solver)),'.txt' open(UNIT=90,FILE=fileout) call mfileout(namefile,solver,filepath,nindsol,indsol) else write(fileout,'(a,a,a,a,a)') filepath(1:lnblnk(filepath)), + problm(1:lnblnk(problm)),'_', + solver(1:lnblnk(solver)),'.txt' open(UNIT=90,FILE=fileout) call mfileout(problm,solver,filepath,nindsol,indsol) end if end if c----------------------------------------------------------------------- c call of the subroutine ACDCC c----------------------------------------------------------------------- timtim = gettim() timtim = gettim() - timtim cputim = gettim() if ( bvptype .eq. 'SPBVP') then if (numjac .and. .not. numbcjac) then CALL ACDCC(neqn,NLBC,NXXDIM,ALEFT,ARIGHT,NFXPNT,FIXPNT,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,full,nmguess, + xguess,nudim,yguess,NMSH, + XX,NUDIM,U,NMAX,LWRKFL,WRK,LWRKIN,IWRK,giveps, + eps,epsmin, + odef,fnumjac,gsub,dgsub, + ckappa1,gamma1,sigma,ckappa,ckappa2,rpar,ipar,iset, + precis, use_c,IFLBVP) elseif (numjac .and. numbcjac) then CALL ACDCC(neqn,NLBC,NXXDIM,ALEFT,ARIGHT,NFXPNT,FIXPNT,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,full,nmguess, + xguess,nudim,yguess,NMSH, + XX,NUDIM,U,NMAX,LWRKFL,WRK,LWRKIN,IWRK,giveps, + eps,epsmin, + odef,fnumjac,gsub,bcnumjac, + ckappa1,gamma1,sigma,ckappa,ckappa2,rpar,ipar,iset, + precis, use_c,IFLBVP) else CALL ACDCC(neqn,NLBC,NXXDIM,ALEFT,ARIGHT,NFXPNT,FIXPNT,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,full,nmguess, + xguess,nudim,yguess,NMSH, + XX,NUDIM,U,NMAX,LWRKFL,WRK,LWRKIN,IWRK,giveps, + eps,epsmin, + odef,odejac,gsub,dgsub, + ckappa1,gamma1,sigma,ckappa,ckappa2,rpar,ipar,iset, + precis, use_c,IFLBVP) end if elseif (bvptype .eq. 'BVP') then print *, 'ACDCCD: ERROR: ', + 'ACDCC cannot solve ', bvptype, ' problems ' stop end if if (iflbvp.gt.0) then print *, 'ACDCC: ERROR: ', + 'ACDCC returned IFLBVP = ', IFLBVP stop endif if (iflbvp.ne.0) then print *, 'ACDCC: WARNING: ', + 'ACDCC RETURNED A SOLUTION FOR EPSMIN = ', EPS, + 'INSTEAD OF', EPSMIN stop endif cputim = gettim() - cputim - timtim c----------------------------------------------------------------------- c print numerical solution in endpoint and integration statistics c----------------------------------------------------------------------- c c printsolout = .true. if (solref) then call solut(neqn,xx(1:nmsh),yref(1:neqn,1:nmsh),nmsh,rpar,ipar) call getscd(mescd,scd,neqn,xx(1:nmsh),yref(1:neqn,1:nmsh), + u(1:neqn,1:nmsh), + problm, tolvec,tol,atol, nmsh,printsolout, + fileout,nindsol,indsol) call report( + driver,problm,solver, + tol,tol, + iset,cputim,scd,mescd,nmsh, + ckappa,ckappa1,ckappa2,gamma1,sigma) else c add the computation of a reference solution in the mesh point CALL APPREFSOL(neqn,nmsh,aleft,aright, + xx,u,uacc,tol,ntol,ltol, + nlbc,linear, wrk, lwrkfl,iwrk,lwrkin, + ckappa, stab_kappa, accmeshfail, + nudim,nxxdim,eps,numjac,numbcjac, rpar, ipar) if (.not. accmeshfail) THEN call getscd(mescd,scd,neqn,xx(1:nmsh),uacc(1:neqn,1:nmsh), + u(1:neqn,1:nmsh), + problm, tolvec,tol,atol, nmsh,printsolout, + fileout,nindsol,indsol) call report( + driver,problm,solver, + tol,tol, + iset,cputim,scd,mescd,nmsh, + ckappa,ckappa1,ckappa2,gamma1,sigma) ELSE call printsol(neqn,nmsh,xx(1:nmsh),u(1:neqn,1:nmsh), + problm,fileout,nindsol,indsol,printsolout) end if end if if (printsolout) then close(90) end if end C======================================================================= C `Test Set for BVP Solvers' solout function TWPBVPLC C======================================================================= C======================================================================= C `Test Set for BVP Solvers' ODE wrappers for TWPBVPLC C======================================================================= c c since in TWPBVPLC the format of the subroutines for the c function f and the Jacobian differ from the format c in the testset, we transform them c c----------------------------------------------------------------------- subroutine odef(n,x,y,f,eps,rpar,ipar) integer n,ipar(*) double precision x,y(*),f(*),rpar(*) double precision EPS double precision :: EPSVAL common /parms/ EPSVAL EPSVAL = EPS call feval(n,x,y(1:n),f(1:n),rpar,ipar) return end subroutine odejac(n,x,y,dfdy,eps,rpar,ipar) implicit none integer n,ipar(*) double precision x,y(*),dfdy(n,*),rpar(*) double precision EPS double precision :: EPSVAL common /parms/ EPSVAL EPSVAL = EPS call jeval(n,x,y,dfdy,rpar,ipar) return end SUBROUTINE GSUB(I,NEQN,Z,G,EPS,RPAR,IPAR) implicit none integer, INTENT(IN) :: i, neqn integer, INTENT(IN) :: ipar(*) double precision, INTENT(IN) :: Z(neqn) double precision, INTENT(IN) :: rpar(*) double precision, INTENT(OUT) :: G double precision EPS double precision :: EPSVAL common /parms/ EPSVAL EPSVAL = EPS call bceval(I,NEQN,Z,G,rpar,ipar) return end SUBROUTINE DGSUB(I,NCOMP,Z,DG,EPS, RPAR,IPAR) implicit none integer i,ncomp, ipar(*) double precision Z(*), RPAR(*), DG(ncomp) double precision EPS double precision :: EPSVAL common /parms/ EPSVAL EPSVAL = EPS call dbceval(I,NCOMP,Z,DG, RPAR,IPAR) return end c c---------------------------------------------------------------------- c Numerical Jacobian c SUBROUTINE FNUMJAC(R, X0, Y0, JF0, EPS, RPAR,IPAR) IMPLICIT NONE INTEGER, INTENT(IN) :: R, IPAR(*) DOUBLE PRECISION, INTENT(IN) :: X0, Y0(R), RPAR(*),EPS DOUBLE PRECISION, INTENT(OUT) :: JF0(1:R,1:R) EXTERNAL odef DOUBLE PRECISION :: DN(R), Y0N(R),YSAFE, DELT, X,FP(R),EPSPREC INTEGER :: I EPSPREC= epsilon(1.0d0) CALL odef(R,X0,Y0,FP,EPS,RPAR,IPAR) Y0N(1:R) = Y0(1:R) DO I=1,R YSAFE=Y0(I) DELT=SQRT(1.D-1*EPSPREC*MAX(1.D-5,ABS(YSAFE))) Y0N(I)=YSAFE+DELT X = X0 CALL odef(R,X,Y0N,DN,EPS,RPAR,IPAR) JF0(1:R,I)=(DN(1:R)-FP(1:R))/DELT Y0N(I)=YSAFE END DO RETURN END SUBROUTINE FNUMJAC c c---------------------------------------------------------------------- c computation of the numerical jacobian for the boundary conditions SUBROUTINE BCNUMJAC(I,R,y,dbc,EPS,RPAR,IPAR) IMPLICIT NONE INTEGER, INTENT(IN) :: I, R DOUBLE PRECISION, INTENT(OUT)::dbc(R) DOUBLE PRECISION, INTENT(IN ):: y(R),EPS INTEGER, INTENT(IN):: ipar(*) DOUBLE PRECISION, INTENT(IN) :: rpar(*) DOUBLE PRECISION :: DN, Y0N(R),YSAFE, DELT, BC, EPSPREC INTEGER :: J EXTERNAL gsub EPSPREC= epsilon(1.0d0) Y0N(1:R) = Y(1:R) call gsub(i,R,y,bc,eps,rpar,ipar) DO J=1,R YSAFE=Y(J) DELT=SQRT(1.D-1*EPSPREC*MAX(1.D-5,ABS(YSAFE))) Y0N(J)=YSAFE+DELT call gsub(i,R,y0N,DN,eps,rpar,ipar) DBC(J)=(DN-BC)/DELT Y0N(J)=YSAFE END DO RETURN END SUBROUTINE BCNUMJAC c c---------------------------------------------------------------------- c computation of the numerical jacobian for the boundary conditions SUBROUTINE APPREFSOL(neqn,nmsh,aleft,aright, + xx,u,uacc,tol,ntol,ltol, + nlbc,linear, wrk, lwrkfl,iwrk,lwrkin, + ckappa, stab_kappa, accmeshfail, + nudim,nxxdim,eps,numjac,numbcjac, rpar, ipar) IMPLICIT NONE INTEGER LWRKFL, LWRKIN, NUDIM, NXXDIM, LISERIES, IPRINT, IND,IWRK, + NTOL, IPAR, LTOL, NLBC, NMSH, NFXPNT, NMAX, + I, NMSHACC, K, NMINIT, IDUM, IFLBVP_A DOUBLE PRECISION TOL, WRK, RPAR, ALEFT, ARIGHT, XX, U, + CKAPPA, FIXPNTACC, XXACC, UACC,YACC DOUBLE PRECISION CKAPPA_A,CKAPPA1_A,GAMMA1_A,SIGMA_A, CKAPPA2_A LOGICAL numjac,numbcjac DIMENSION WRK(*), IWRK(*), rpar(*),ipar(*) DIMENSION U(NUDIM,NXXDIM), XX(NXXDIM) DIMENSION UACC(NUDIM,NXXDIM), XXACC(NXXDIM) DIMENSION YACC(NUDIM,NXXDIM) DOUBLE PRECISION XGUESS(NXXDIM),YGUESS(NUDIM,NXXDIM) DIMENSION TOL(NUDIM), LTOL(NUDIM), FIXPNTACC(NXXDIM) INTEGER indnms_a,iseries_a,nmguess, neqn INTEGER iset_a(8) DOUBLE PRECISION precis(3) LOGICAL LINEAR, GIVEU, GIVMSH logical useC EXTERNAL ACDCC EXTERNAL odef EXTERNAL odejac EXTERNAL gsub EXTERNAL dgsub EXTERNAL fnumjac EXTERNAL bcnumjac logical full, ACCMESHFAIL logical stab_kappa LOGICAL use_c, comp_c, PRINTSOL, GIVEPS, giv_u integer ureset DOUBLE PRECISION EPS, EPSMIN, EPSIN c this common need to be defined in order to run TWPBVPC successfully c give information about some parameters common/algprs/ nminit, iprint, idum, use_c, comp_c common/acgu/ giv_u, ureset precis(3) = epsilon(1.0d0) precis(2) = huge(1.0d0) precis(1) = tiny(1.0d0) PRINTSOL = .FALSE. full=.false. use_c = .true. useC=use_c NMSHACC= 2*NMSH-1 IF (NMSHACC .gt. NXXDIM ) THEN ACCMESHFAIL = .true. RETURN ENDIF DO I=2,NMSH-1 FIXPNTACC(2*(I-1)-1)=(XX(I-1)+XX(I))/2.0D0 FIXPNTACC(2*(I-1))=XX(I) END DO FIXPNTACC(2*(NMSH-2)+1)=(XX(NMSH-1)+XX(NMSH))/2.0 NFXPNT= 2*(NMSH-2)+1 Xguess(1) = aleft DO I=2,2*(NMSH-2)+2 Xguess(I) =FIXPNTACC(I-1) END DO Xguess(NMSHACC) = aright XXACC(1:NMSHACC)=XGUESS(1:NMSHACC) XXACC(1) = aleft XXACC(NMSHACC) = aright nmguess=nmshacc givmsh = .TRUE. giveu = .TRUE. giveps = .true. giv_u=giveu epsmin= eps Call acInterp(Neqn, Nmshacc, Xguess, Nudim, yguess,nudim, * Nmsh, Xx, U) YACC(1:neqn,1:nmshacc) = yguess(1:neqn,1:nmshacc) DO I=1,NTOL TOL(I) = TOL(I)/1.0D2 END DO DO ind=1,LWRKFL WRK(ind) = 0d0 ENDDO DO ind=1,LWRKIN IWRK(ind) = 0 ENDDO accmeshfail = .false. if (numjac .and. .not. numbcjac) then CALL ACDCC(neqn,NLBC,NXXDIM,ALEFT,ARIGHT,NFXPNT,FIXPNTACC,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,full,nmguess, + xguess,nudim,yguess,NMSHACC, + XXACC,NUDIM,YACC,NMAX,LWRKFL,WRK,LWRKIN,IWRK,giveps, + eps,epsmin, + odef,fnumjac,gsub,dgsub, + ckappa1_a,gamma1_a,sigma_a,ckappa_a,ckappa2_a, + rpar,ipar,iset_a, + precis, use_c,IFLBVP_a) elseif (numjac .and. numbcjac) then CALL ACDCC(neqn,NLBC,NXXDIM,ALEFT,ARIGHT,NFXPNT,FIXPNTACC,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,full,nmguess, + xguess,nudim,yguess,NMSHACC, + XXACC,NUDIM,YACC,NMAX,LWRKFL,WRK,LWRKIN,IWRK,giveps, + eps,epsmin, + odef,fnumjac,gsub,bcnumjac, + ckappa1_a,gamma1_a,sigma_a,ckappa_a,ckappa2_a, + rpar,ipar,iset_a, + precis, use_c,IFLBVP_a) else CALL ACDCC(neqn,NLBC,NXXDIM,ALEFT,ARIGHT,NFXPNT,FIXPNTACC,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,full,nmguess, + xguess,nudim,yguess,NMSHACC, + XXACC,NUDIM,YACC,NMAX,LWRKFL,WRK,LWRKIN,IWRK,giveps, + eps,epsmin, + odef,odejac,gsub,dgsub, + ckappa1_a,gamma1_a,sigma_a,ckappa_a,ckappa2_a, + rpar,ipar,iset_a, + precis, use_c,IFLBVP_a) endif IF (IFLBVP_A .gt. 0 ) THEN ACCMESHFAIL = .true. RETURN ENDIF IF (NMSHACC .gt. NXXDIM ) THEN ACCMESHFAIL = .true. RETURN ENDIF IF (NMSHACC .eq. NMSH) THEN ACCMESHFAIL = .true. RETURN ENDIF IF (NMSHACC .eq. 0 .or. IFLBVP_A .GT. 0) THEN ACCMESHFAIL = .true. RETURN ENDIF IF (comp_c) THEN STAB_KAPPA = abs(CKAPPA-CKAPPA_A)/CKAPPA_A < 5d-2 ELSE STAB_KAPPA=.FALSE. END IF UACC(1:neqn,1)=YACC(1:neqn,1) K = 1 DO I=1,NMSH-1 DO WHILE(DABS(XX(I)-XXACC(K)) .gt. 1.0D-10) K = K + 1 IF (K .ge. NMSHACC) THEN accmeshfail = .true. RETURN END IF END DO UACC(1:neqn,I)=YACC(1:neqn,K) END DO UACC(1:neqn,nmsh)=YACC(1:neqn,nmshacc) RETURN END