research/incomp-flame-1d
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

added convergence criterion

Browse source
This commit is contained in:
ignis 2019-08-31 09:54:34 +09:00
parent ea467d5e27
commit 61ed782485
2 changed files with 50 additions and 23 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
code/m_parameters.f90
View file

@ -24,6 +24,8 @@ module m_parameters
real, parameter :: me=1.00e-20
real :: absolute_tolerence=1e-8
REAL :: hx,dt,vis,sc,diff,tf,t_now,t_uf,dt_uf
REAL :: lewis

71
code/ysolve.f90
View file

@ -6,7 +6,7 @@
IMPLICIT NONE
REAL, DIMENSION(:), ALLOCATABLE :: u
REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf
REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf,yold
REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
REAL, DIMENSION(:), ALLOCATABLE :: dm
@ -59,7 +59,21 @@
END IF
END SUBROUTINE parse
!------------------------------------------------------------------------
LOGICAL FUNCTION converged (x0, x1)
REAL, DIMENSION(:,:) :: x0, x1
REAL :: residual = 0.
residual = sum(abs(x0-x1))
converged = residual < absolute_tolerence
END FUNCTION converged
!------------------------------------------------------------------------
SUBROUTINE solve
INTEGER :: i,j,k,savenum
@ -73,35 +87,44 @@
t_uf=0.
DO
CALL SET_BC
CALL RK4(fns)
IF(t_now.ge.tf) EXIT
IF(t_uf.ge.dt_uf) THEN
pflold=pflame
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1)*hx
uxt(i)=y1(i,1)
ENDDO
CALL dfnonp(nx,hx,uxt,duxt,1,1)
delf=1./MAXVAL(ABS(duxt))
oldu=u(1)
u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame)
t_uf=0.
WRITE(*,'(a3,f8.3,a10,f6.3,a10,f6.3,a10,f7.4,a10,f7.4,a8,f7.4)') &
' T:',t_now,' // Tar_L:',l_0*tar_lo,' // cur_L:',pflame/hx/REAL(nx)*l_0, &
' // Old_U:',oldu,' // New_U:',u(1),' // L_f:',delf
! WRITE(*,'(a7,f7.4,a7,f7.4,a10,f7.4)') ' cur_U:',oldu,' // dU:',u(1,1,1)-oldu,' // new_U:',u(1,1,1)
! WRITE(*,*)
ENDIF
IF (t_now.ge.tf) EXIT
IF (converged(yold, y1)) EXIT
ncyc=ncyc+1
t_uf=t_uf+dt
t_now=t_now+dt
yold = y1
CALL SET_BC
CALL RK4(fns)
IF (t_uf.ge.dt_uf) THEN
pflold=pflame
pflame=0.
DO i=1,nx
pflame=pflame+y1(i,1)*hx
uxt(i)=y1(i,1)
ENDDO
CALL dfnonp(nx,hx,uxt,duxt,1,1)
delf=1./MAXVAL(ABS(duxt))
oldu=u(1)
u=u+0.5*(hx*REAL(nx)*tar_lo-pflame)+0.5*(pflold-pflame)
t_uf=0.
WRITE(*,'(a3,f8.3,a10,f6.3,a10,f6.3,a10,f7.4,a10,f7.4,a8,f7.4)') &
' T:',t_now,' // Tar_L:',l_0*tar_lo,' // cur_L:',pflame/hx/REAL(nx)*l_0, &
' // Old_U:',oldu,' // New_U:',u(1),' // L_f:',delf
ENDIF
IF (MOD(ncyc,int_pr).eq.0) THEN
! WRITE(*,'(a2,f8.3,a9,f10.7,a11,i6,a7,f9.5)') &
! 'T:',t_now,' // dT:',dt,' // NCYC:',ncyc,' // U:',u(1,1,1)
ENDIF
ENDIF
ENDDO
CALL write_sd
@ -301,6 +324,7 @@
ALLOCATE(y1(nx,2),STAT=ierr) ; y1=0.
ALLOCATE(y2(nx,2),STAT=ierr) ; y2=0.
ALLOCATE(yf(nx,2),STAT=ierr) ; yf=0.
ALLOCATE(yold(nx,2),STAT=ierr) ; yold=0.
ALLOCATE(uxt(nx),STAT=ierr) ; uxt=0.
ALLOCATE(duxt(nx),STAT=ierr) ; duxt=0.
ALLOCATE(dm(nx),STAT=ierr) ; dm=diff
@ -313,6 +337,7 @@
DEALLOCATE(y1)
DEALLOCATE(y2)
DEALLOCATE(yf)
DEALLOCATE(yold)
DEALLOCATE(uxt)
DEALLOCATE(duxt)
DEALLOCATE(dm)