From 61ed78248552318ff7426e07879a1b28e411beff Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Sat, 31 Aug 2019 09:54:34 +0900
Subject: [PATCH] added convergence criterion

---
 code/m_parameters.f90 |  2 ++
 code/ysolve.f90       | 71 +++++++++++++++++++++++++++++--------------
 2 files changed, 50 insertions(+), 23 deletions(-)

diff --git a/code/m_parameters.f90 b/code/m_parameters.f90
index fa78b7c..e8ce0d0 100644
--- a/code/m_parameters.f90
+++ b/code/m_parameters.f90
@@ -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
diff --git a/code/ysolve.f90 b/code/ysolve.f90
index b123150..32ddcd9 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -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)