From d5bb7d388a84475f65e5be512fca0940e7647b1f Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Sun, 1 Sep 2019 03:34:37 +0900
Subject: [PATCH] added an equation

---
 code/ysolve.f90 | 24 ++++++++++++++++--------
 1 file changed, 16 insertions(+), 8 deletions(-)

diff --git a/code/ysolve.f90 b/code/ysolve.f90
index 25ea98f..3933233 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -223,7 +223,7 @@
 
       OPEN (305,FILE='sfield.dat')
       DO i=1,nx
-       WRITE (305,'(e30.20)') y1(i,1)
+       WRITE (305,'(e30.20)') y1(i,2)
       ENDDO
       CLOSE (305)
 
@@ -237,7 +237,7 @@
       
       DO i=1,nx
        yr=y1(i,1)
-       theta=y1(i,2)
+       theta=y1(i,3)
        c=1.-yr
        ux(1,i)=yr
 
@@ -274,7 +274,8 @@
 
      SUBROUTINE SET_BC
       y1(1,1)=1.-ctmp
-      y1(1,2)=0.
+      y1(1,2)=1.
+      y1(1,3)=0.
      END SUBROUTINE SET_BC
 
      SUBROUTINE SET_IC
@@ -304,8 +305,11 @@
        y1(i,1)=(1.-xi)   ! reactant mass fraction
       ENDDO
 
+      ! initialize intermediate species field
+      y1(:,2) = y1(:,1)
+
       ! initialize temperature field
-      y1(:,2) = 1.0 - y1(:,1)
+      y1(:,3) = 1.0 - y1(:,1)
 
       ! calculate flame position
       pflame=0.
@@ -409,7 +413,7 @@
 ! x-direction
       DO i=1,nx
          ux(1,i)=r1(i,1)                              ! Y
-         ux(2,i)=r1(i,2)                              ! T
+         ux(2,i)=r1(i,2)                              ! Y
          ux(3,i)=r1(i,3)                              ! T
       ENDDO
 
@@ -418,18 +422,22 @@
       CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),3,1)
 
       DO i=1,nx
-        wrate=rate_1step(ux(1,i), ux(2,i))
+        wrate=rate_1step(ux(1,i), ux(3,i))
 
         ! - u*dY/dx + D*d2Y/d2x
         f(i,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate
 
+        ! - u*dY/dx + D*d2Y/d2x
+        f(i,2) = - ( u(i)*dux(2,i) ) + diff * d2ux(2,i) - wrate
+
         ! -  u*dT/dx + D*d2T/d2x
-        f(i,2) = - ( u(i)*dux(2,i) ) + (lewis * diff) * d2ux(2,i) + wrate
+        f(i,3) = - ( u(i)*dux(3,i) ) + (lewis * diff) * d2ux(3,i) + wrate
 
         ! Boundary conditions
         IF (i.eq.nx) THEN
             f(nx,1) = -wrate - u(nx)*dux(1,nx)
-            f(nx,2) = wrate - u(nx)*dux(2,nx)
+            f(nx,2) = -wrate - u(nx)*dux(2,nx)
+            f(nx,3) = wrate - u(nx)*dux(3,nx)
         ENDIF
       ENDDO