c c This file is part of the Test Set for BVP solvers c http://www.dm.uniba.it/~bvpsolvers/ c c generic TWPBVPLC.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/twpbvplcd.f c c c----------------------------------------------------------------------- c program twpbvplcd 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 PARAMETER(NUDIM=20,NXXDIM=10000) PARAMETER(LWRKFL=26*NUDIM+16*NUDIM*NUDIM+14*NUDIM*NXXDIM * +5*NUDIM*NUDIM*NXXDIM+4*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) INTEGER iseries(liseries) INTEGER nmguess INTEGER iset(7) DOUBLE PRECISION precis(3) integer n0, maxmesh LOGICAL LINEAR, GIVEU, GIVMSH logical useC logical full, ACCMESHFAIL LOGICAL use_c, comp_c, tolvec EXTERNAL TWPBVPLC EXTERNAL feval EXTERNAL jeval EXTERNAL bceval EXTERNAL dbceval EXTERNAL fnumjac EXTERNAL bcnumjac INTEGER lnblnk EXTERNAL lnblnk 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 character fullnm*25, problm*15, bvptype*5 character driver*9, solver*8 parameter (driver = 'twpbvplcd', solver='TWPBVPLC') 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 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)=' ' 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 *, 'TWPBVPLCD: ERROR: ', + 'TWPBVPLC cannot solve higher order problems' stop endif if ( bvptype(1:lnblnk(bvptype)).ne.'BVP' .and. + bvptype(1:lnblnk(bvptype)) .ne. 'SPBVP' )then print *, 'TWPBVPLCD: ERROR: ', + 'unknown Test Set problem bvptype ', bvptype 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 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) 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 TWPBVPLC c----------------------------------------------------------------------- timtim = gettim() timtim = gettim() - timtim cputim = gettim() if (bvptype(1:3).eq.'BVP' .or. + bvptype(1:5).eq.'SPBVP') then if (numjac .and. .not. numbcjac) then CALL TWPBVPLC(neqn,NLBC,ALEFT,ARIGHT,NFXPNT,FIXPNT,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,NMSH,NXXDIM, + XX,NUDIM,U,NMAX,LWRKFL,WRK,LWRKIN,IWRK,precis, + feval,fnumjac,bceval,dbceval, + ckappa1,gamma1,sigma,ckappa,ckappa2,rpar,ipar,IFLBVP, + liseries, iseries, indnms,full,useC,nmguess, + xguess,nudim,yguess,iset) elseif (numjac .and. numbcjac) then CALL TWPBVPLC(neqn,NLBC,ALEFT,ARIGHT,NFXPNT,FIXPNT,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,NMSH,NXXDIM, + XX,nudim,U,NMAX,LWRKFL,WRK, + LWRKIN,IWRK,precis, + feval,fnumjac,bceval,bcnumjac, + ckappa1,gamma1,sigma,ckappa,ckappa2,rpar,ipar,IFLBVP, + liseries, iseries, indnms,full,useC,nmguess, + xguess,nudim,yguess,iset) else CALL TWPBVPLC(neqn,NLBC,ALEFT,ARIGHT,NFXPNT,FIXPNT,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,NMSH,NXXDIM, + XX,NUDIM,U,NMAX,LWRKFL,WRK,LWRKIN,IWRK,precis, + feval,jeval,bceval,dbceval, + ckappa1,gamma1,sigma,ckappa,ckappa2,rpar,ipar,IFLBVP, + liseries, iseries, indnms,full,useC,nmguess, + xguess,nugdim,yguess,iset) end if else print *, 'TWPBVPLCD: ERROR: ', + 'unknown Test Set problem bvptype ', bvptype stop end if if (iflbvp.gt.0) then print *, 'TWPBVPLC: ERROR: ', + 'TWPBVPLC returned IFLBVP = ', IFLBVP stop endif cputim = gettim() - cputim - timtim c----------------------------------------------------------------------- c print numerical solution in endpoint and integration statistics c----------------------------------------------------------------------- c c printsolout = .true. maxmesh = 0 DO k=1,liseries maxmesh = max(maxmesh,iseries(k)) ENDDO iset(7) = maxmesh 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, 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 c---------------------------------------------------------------------- c Numerical Jacobian c SUBROUTINE FNUMJAC(R, X0, Y0, JF0, RPAR,IPAR) IMPLICIT NONE INTEGER, INTENT(IN) :: R, IPAR(*) DOUBLE PRECISION, INTENT(IN) :: X0, Y0(R), RPAR(*) DOUBLE PRECISION, INTENT(OUT) :: JF0(1:R,1:R) EXTERNAL feval DOUBLE PRECISION :: DN(R), Y0N(R),YSAFE, DELT, X,FP(R),EPSPREC c DOUBLE PRECISION :: SQRTEPS INTEGER :: I EPSPREC= epsilon(1.0d0) c SQRTEPS = sqrt(EPSPREC) CALL feval(R,X0,Y0,FP,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))) c IF ( ABS(YSAFE) > 0D0 ) THEN c DELT = SQRTEPS*ABS(YSAFE) c ELSE c DELT = SQRTEPS c END IF Y0N(I)=YSAFE+DELT X = X0 CALL feval(R,X,Y0N,DN,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,RPAR,IPAR) IMPLICIT NONE INTEGER, INTENT(IN) :: I, R DOUBLE PRECISION, INTENT(OUT)::dbc(R) DOUBLE PRECISION, INTENT(IN ):: y(R) INTEGER, INTENT(IN):: ipar(*) DOUBLE PRECISION, INTENT(IN) :: rpar(*) DOUBLE PRECISION :: DN, Y0N(R),YSAFE, DELT, BC, EPSPREC c DOUBLE PRECISION :: SQRTEPS INTEGER :: J EXTERNAL bceval EPSPREC= epsilon(1.0d0) c SQRTEPS = sqrt(EPSPREC) Y0N(1:R) = Y(1:R) call bceval(i,R,y,bc,rpar,ipar) DO J=1,R YSAFE=Y(J) DELT=SQRT(1.D-1*EPSPREC*MAX(1.D-5,ABS(YSAFE))) c IF ( ABS(YSAFE) > 0D0 ) THEN c DELT = SQRTEPS*ABS(YSAFE) c ELSE c DELT = SQRTEPS c END IF Y0N(J)=YSAFE+DELT call bceval(i,R,y0N,DN,rpar,ipar) DBC(J)=(DN-BC)/DELT Y0N(J)=YSAFE END DO RETURN END SUBROUTINE BCNUMJAC c c---------------------------------------------------------------------- c Computation of the reference solution c c SUBROUTINE APPREFSOL(neqn,nmsh,aleft,aright, + xx,u,uacc,tol,ntol,ltol, + nlbc,linear, wrk, lwrkfl,iwrk,lwrkin, + ckappa, stab_kappa, accmeshfail, + nudim,nxxdim,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) PARAMETER (LISERIES=100) dimension iseries_a(liseries) INTEGER indnms_a,iseries_a,nmguess, neqn INTEGER iset_a(6) DOUBLE PRECISION precis(3) LOGICAL LINEAR, GIVEU, GIVMSH logical useC EXTERNAL TWPBVPLC EXTERNAL feval EXTERNAL jeval EXTERNAL bceval EXTERNAL dbceval EXTERNAL fnumjac EXTERNAL bcnumjac logical full, ACCMESHFAIL logical stab_kappa LOGICAL use_c, comp_c, PRINTSOL 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 precis(3) = epsilon(1.0d0) precis(2) = huge(1.0d0) precis(1) = tiny(1.0d0) PRINTSOL = .FALSE. full=.false. 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) = XX(1) DO I=2,2*(NMSH-2)+2 Xguess(I) =FIXPNTACC(I-1) END DO Xguess(NMSHACC) = XX(NMSH) XXACC(1:NMSHACC)=XGUESS(1:NMSHACC) XXACC(1) = aleft XXACC(NMSHACC) = aright nmguess=nmshacc givmsh = .TRUE. giveu = .TRUE. Call Interp(Neqn, Nmshacc, Xguess, Nudim, yguess,nudim, * Nmsh, Xx, U) YACC(1:neqn,1:nmsh) = yguess(1:neqn,1:nmsh) 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 TWPBVPLC(neqn,NLBC,ALEFT,ARIGHT,NFXPNT,FIXPNTACC,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,NMSHACC,NXXDIM, + XXACC,NUDIM,YACC,NMAX,LWRKFL,WRK,LWRKIN,IWRK,precis, + feval,fnumjac,bceval,dbceval, + ckappa1_a,gamma1_a,sigma_a,ckappa_a,ckappa2_a, + rpar,ipar,IFLBVP_a, + liseries, iseries_a, indnms_a,full,useC,nmguess, + xguess,nudim,yguess,iset_a) elseif (numjac .and. numbcjac) then CALL TWPBVPLC(neqn,NLBC,ALEFT,ARIGHT,NFXPNT,FIXPNTACC,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,NMSHACC,NXXDIM, + XXACC,NUDIM,YACC,NMAX,LWRKFL,WRK,LWRKIN,IWRK,precis, + feval,fnumjac,bceval,bcnumjac, + ckappa1_a,gamma1_a,sigma_a,ckappa_a,ckappa2_a, + rpar,ipar,IFLBVP_a, + liseries, iseries_a, indnms_a,full,useC,nmguess, + xguess,nudim,yguess,iset_a) else CALL TWPBVPLC(neqn,NLBC,ALEFT,ARIGHT,NFXPNT,FIXPNTACC,NTOL, + LTOL,TOL,LINEAR,GIVMSH,GIVEU,NMSHACC,NXXDIM, + XXACC,NUDIM,YACC,NMAX,LWRKFL,WRK,LWRKIN,IWRK,precis, + feval,jeval,bceval,dbceval, + ckappa1_a,gamma1_a,sigma_a,ckappa_a,ckappa2_a, + rpar,ipar,IFLBVP_a, + liseries, iseries_a, indnms_a,full,useC,nmguess, + xguess,nudim,yguess,iset_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 .gt. 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