solve T

added debug flag to makefile
added solution variable T
2019-03-28 15:58:09 +09:00 · 2019-03-28 14:43:35 +09:00 · 2019-03-27 13:35:36 +09:00 · 2019-03-27 13:34:32 +09:00 · 2019-03-27 10:39:43 +09:00 · 2019-03-24 13:56:57 +09:00
7 changed files with 220 additions and 130 deletions
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@ -16,4 +16,6 @@ test:
  script:
    - cd ${CI_PROJECT_DIR}
    - cd code
-    - make test
+    - make clean && make
+    - cd sample/4pi-IC1
+    - ../../ex
--- a/code/makefile
+++ b/code/makefile
@ -1,17 +1,25 @@
 flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8  -march=native 

+ifdef CHECK
+flags += -fcheck=all
+endif
+
+ifdef DEBUG
+flags += -g
+endif
+
 compiler = gfortran

-ex  : test.o ysolve.mod compact.mod
+ex  : test.o ysolve.o Compact.o
 	${compiler} -o ex test.o ysolve.o Compact.o 

-test.o  : test.f90 ysolve.mod
+test.o  : test.f90 ysolve.mod compact.mod
 	${compiler} -c ${flags} test.f90

-ysolve.mod: ysolve.f90 compact.mod
+ysolve.o ysolve.mod : ysolve.f90 compact.mod
 	${compiler} -c ${flags} ysolve.f90

-compact.mod: Compact.f90
+Compact.o compact.mod : Compact.f90
 	${compiler} -c ${flags} Compact.f90

 test: ex
--- a/code/sample/4pi-IC1/itape
+++ b/code/sample/4pi-IC1/itape
@ -16,4 +16,5 @@
     c_ref     0.01  !0.003     ! 
    min_wr     0.        ! 5.0e-14
   prof_wr     1.0       
+     lewis     1.0

--- a/code/sample/itape.py
+++ b/code/sample/itape.py
@ -0,0 +1,83 @@
+#!/bin/python
+
+import os, sys
+
+import argparse
+
+parser = argparse.ArgumentParser()
+file_parser = argparse.ArgumentParser()
+
+param_dict = {
+"nx"		:512,
+"l_0"		:4.,
+"int_pr"	:400,
+"tar_lo"	:0.60,
+"dt"		:0.001,
+"sc"		:0.75,
+"vis"		:0.020,
+"pre"		:2.10E+4,
+"ac"		:26.7,
+"bc"		:3.,
+"u0"		:0.21,
+"tf"		:100.,
+"dt_uf"		:0.4,
+"ctmp"		:0.,
+"c_cut"		:0.001,
+"c_ref"		:0.01,
+"min_wr"	:0.,
+"prof_wr"	:1.0,
+"rvis"		:1.0,
+"cstar"		:1.0,
+"sigw"		:100.0,
+}
+
+file_parser.add_argument("--default", help="file containing override default values")
+
+file_args, others = file_parser.parse_known_args()
+
+if file_args.default is not None:
+    with open(file_args.default) as df:
+        try:
+            for line in df:
+                k, v = line.split()[:2]
+                if k in param_dict:
+                    param_dict[k] = v
+        except ValueError:
+          pass
+
+
+for k, v in param_dict.iteritems():
+    parser.add_argument("--"+k, default=v)
+
+args = parser.parse_args(others)
+
+params = vars(args)
+
+itape = '''        nx     {0[nx]}
+       l_0     {0[l_0]}
+    int_pr     {0[int_pr]}
+    tar_lo     {0[tar_lo]}
+        dt     {0[dt]}
+        sc     {0[sc]}
+       vis     {0[vis]}
+       pre     {0[pre]}
+        ac     {0[ac]}
+        bc     {0[bc]}
+        u0     {0[u0]}
+        tf     {0[tf]}
+     dt_uf     {0[dt_uf]}
+      ctmp     {0[ctmp]}
+     c_cut     {0[c_cut]}
+     c_ref     {0[c_ref]}
+    min_wr     {0[min_wr]}
+   prof_wr     {0[prof_wr]}
+      rvis     {0[rvis]}
+     cstar     {0[cstar]}
+      sigw     {0[sigw]} '''
+
+print itape.format( params )
+
+# for i in range(15,60,5):
+    # with open("itape", 'w') as itape:
+        # itape.write(format.format(float(i)/10.0))
+    # os.system("../../ex &>"+str(i)+".txt")
--- a/code/test.f90
+++ b/code/test.f90
@ -6,5 +6,6 @@
    PROGRAM test
     USE ysolve
     IMPLICIT NONE
+     CALL parse
     CALL solve
    END PROGRAM test
--- a/code/test.sh
+++ b/code/test.sh
@ -1,4 +1,4 @@
 cd sample/4pi-IC1
 ls
 ../../ex
-diff sfield.dat sfield.diff
+diff -y sfield.dat sfield.diff
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@ -4,12 +4,15 @@
     PRIVATE
     REAL, PARAMETER :: pi=3.14159265358979323846
     REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
+     REAL :: lewis
     REAL :: c_cut,c_ref,refwr,minf,tar_lo,u0,ctmp,l_0,lo_flm=0.
     REAL :: er_lof=0., erdot=0.,min_wr,prof_wr
     REAL :: pflame,pflold,oldu
     INTEGER :: ncyc=0,int_pr,nx
-     REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,y1,y2,yf
+     REAL, DIMENSION(:), ALLOCATABLE :: u
+     REAL, DIMENSION(:,:), ALLOCATABLE :: y1,y2,yf
     REAL, DIMENSION(:), ALLOCATABLE :: uxt,duxt
+     REAL, DIMENSION(:), ALLOCATABLE :: dm

     INTEGER :: istage
     REAL, DIMENSION(5) :: a=(/ 970286171893./4311952581923.,    &
@ -23,8 +26,38 @@
                                2114624349019./3568978502595.,   &
                                5198255086312./14908931495163. /)

-     PUBLIC :: solve
+     LOGICAL :: read_itape, read_stdin
+     CHARACTER(100) :: itape_name
+     PUBLIC :: solve, parse
    CONTAINS
+!------------------------------------------------------------------------
+     SUBROUTINE parse
+        CHARACTER(100) :: num1char
+
+        !First, make sure the right number of inputs have been provided
+        IF(COMMAND_ARGUMENT_COUNT().EQ.0)THEN
+          ! read from itape
+          read_itape = .true.
+        ELSE IF(COMMAND_ARGUMENT_COUNT().EQ.1)THEN
+          CALL GET_COMMAND_ARGUMENT(1,num1char)   !first, read in the two values
+        ELSE
+          WRITE(*,*)'ERROR, TOO MANY COMMAND-LINE ARGUMENTS(MORE THAN ONE). STOPPING'
+          STOP
+        ENDIF
+
+        read_stdin = (num1char=="-")
+
+        IF(read_itape)THEN
+          itape_name = "itape"
+          WRITE(*,*) "Read from itape"
+        ELSE IF(read_stdin)THEN
+          WRITE(*,*) "Read from STDIN"
+        ELSE 
+          itape_name = num1char
+          WRITE(*,*) "Read from " // num1char
+        END IF
+
+     END SUBROUTINE parse
 !------------------------------------------------------------------------
     SUBROUTINE solve
      INTEGER :: i,j,k,savenum
@ -46,17 +79,17 @@
        pflold=pflame
        pflame=0.
        DO i=1,nx
-         pflame=pflame+y1(i,1,1)*hx
-         uxt(i)=y1(i,1,1)
+         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,1,1)
+        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,1,1),' // L_f:',delf
+         ' // 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(*,*)
@ -81,9 +114,9 @@
     END SUBROUTINE solve

     SUBROUTINE write_sd
-      REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
+      REAL :: c,yr,theta,wrate,dely,sdr,sdd,sd,uu,onelw,dd
      INTEGER :: i,j,k,nd
-      REAL, DIMENSION(1,nx) :: ux,dux,d2ux
+      REAL, DIMENSION(2,nx) :: ux,dux,d2ux
      REAL, DIMENSION(10,nx) :: sav
        
      sav=0.
@ -97,21 +130,22 @@
      WRITE(500,*) ' "Sdd","Sd","(1/C)/(dC/dx)","DIV(rho*Dmu*Gra(C))"'

      DO i=1,nx
-       ux(1,i)=y1(i,1,1)                               ! Yr
+       ux(1,i)=y1(i,1)                               ! Yr
+       ux(2,i)=y1(i,2)                               ! T
       IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
      ENDDO
          
-      nd=1
-      CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
-      nd=1 
-      CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) 
+      nd=2
+      CALL dfnonp(nx,hx,ux(:,:),dux(:,:),nd,1)
+      CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),nd,1)

      DO i=1,nx
       yr=ux(1,i)
+       theta=ux(2,i)
       c=1.-yr
       IF (c.lt.0.) c=0.

-       wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr)))   !wrate
+       wrate=pre*yr*exp(-ac/(1.+bc*(theta)))   !wrate

 !       IF (c.le.c_cut) THEN
 !        wrate=min_wr
@ -119,22 +153,22 @@
 !         ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
 !       ENDIF

-       IF (c.le.c_ref) THEN
+       IF (theta.le.c_ref) THEN
         wrate=min_wr
-        IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
+        IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-c_ref))+ &
        min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
       ENDIF

       dely=ABS(dux(1,i))
       sdr=wrate/dely
-       sdd=-diff*d2ux(1,i)/dely
+       sdd=-dm(i)*d2ux(1,i)/dely
       IF (dely.eq.0.) THEN
        sdr=0.; sdd=0.
       ENDIF
       sd=sdr+sdd
-       uu=u(1,1,1)
+       uu=u(1)
       onelw=(-dux(1,i))/c
-       dd=-diff*d2ux(1,i)
+       dd=-dm(i)*d2ux(1,i)
       if (c.eq.0.) onelw=0.

       sav(1,i)=sav(1,i)+yr
@ -159,38 +193,39 @@

      OPEN (305,FILE='sfield.bin',form='unformatted',status='unknown')
      DO i=1,nx
-       WRITE (305) y1(i,1,1)
+       WRITE (305) y1(i,1)
      ENDDO
      CLOSE (305)

      OPEN (305,FILE='sfield.dat')
      DO i=1,nx
-       WRITE (305,'(e30.20)') y1(i,1,1)
+       WRITE (305,'(2e30.20)') y1(i,1), y1(i,2)
      ENDDO
      CLOSE (305)

     END SUBROUTINE save_final_field

     SUBROUTINE write_pre
-      REAL :: yr,c,dy,maxdy=0.,del_f
+      REAL :: theta,yr,c,dy,maxdy=0.,del_f
      REAL :: S_L=0.,wrate
      INTEGER :: i
      REAL, DIMENSION(1,nx) :: ux, dux
      
      DO i=1,nx
-       yr=y1(i,1,1)
+       yr=y1(i,1)
+       theta=y1(i,2)
       c=1.-yr
       ux(1,i)=yr

-       wrate=pre*yr*exp(-ac/(1.+bc*(1.-yr)))   !wrate
+       wrate=pre*yr*exp(-ac/(1.+bc*theta))   !wrate
 !       IF (c.le.c_cut) THEN
 !        wrate=min_wr
 !        IF (c.gt.c_ref) wrate= &
 !         ((exp((c-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
 !       ENDIF
-       IF (c.le.c_ref) THEN
+       IF (theta.le.c_ref) THEN
         wrate=min_wr
-        IF (c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
+        IF (theta.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(theta-c_ref))+ &
        min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
       ENDIF

@ -207,7 +242,7 @@
      ENDDO
      del_f=1./maxdy
       
-      WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1,1,1)
+      WRITE(*,'(a13,e14.8,a25,e14.8)') ' Grid size : ',hx,' / Laminar flame speed : ',u(1)
      WRITE(*,'(a19,e14.8,a3,f9.5,a20)') &
       ' Flame thickness : ',del_f,' / ',del_f/hx,' grids in the flame.'
      WRITE(*,*)
@ -215,7 +250,7 @@

     SUBROUTINE SET_BC
      
-      y1(1,1,1)=1.-ctmp
+      y1(1,1)=1.-ctmp
     END SUBROUTINE SET_BC

     SUBROUTINE SET_IC
@ -236,12 +271,14 @@
       ENDIF

 !       y1(i,1,1)=(1.-xi)*1.
-       y1(i,1,1)=(1.-xi)   ! reactant mass fraction
+       y1(i,1)=(1.-xi)   ! reactant mass fraction
      ENDDO

+      y1(:,2) = 1.0 - y1(:,1)
+
      pflame=0.
      DO i=1,nx
-       pflame=pflame+y1(i,1,1)*hx
+       pflame=pflame+y1(i,1)*hx
      ENDDO

     END SUBROUTINE SET_IC
@ -249,9 +286,14 @@
     SUBROUTINE READ_INTRO
       CHARACTER(LEN=8) :: cdum
       INTEGER :: itape=300,otape=301,ierr
+       INTEGER :: i

+       IF (read_stdin) THEN
+         itape=5
+       ELSE
+         OPEN(itape,FILE=itape_name)
+       END IF

-       OPEN(itape,FILE='itape')
       OPEN(otape,FILE='otape')
        READ(itape,*) cdum,nx
        WRITE(otape,*) cdum,nx
@ -291,8 +333,12 @@
        WRITE(otape,*) cdum,min_wr
        READ(itape,*) cdum,prof_wr
        WRITE(otape,*) cdum,prof_wr
+        READ(itape,*) cdum,lewis
+        WRITE(otape,*) cdum,lewis

-       CLOSE(itape)
+       IF (.not.read_stdin) THEN
+         CLOSE(itape)
+       END IF

 !       hx=l_0*pi/REAL(nx)
       hx=l_0*pi/REAL(nx-1)
@ -312,12 +358,13 @@
       refwr=pre*1.*exp(-ac/(1.+bc*c_ref))          
       l_0=l_0*pi

-       ALLOCATE(u(nx,1,1),STAT=ierr)   ; u=0.
-       ALLOCATE(y1(nx,1,1),STAT=ierr)  ; y1=0.
-       ALLOCATE(y2(nx,1,1),STAT=ierr)  ; y2=0.
-       ALLOCATE(yf(nx,1,1),STAT=ierr)  ; yf=0.
+       ALLOCATE( u(nx),STAT=ierr)    ; u=0.
+       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(uxt(nx),STAT=ierr)     ; uxt=0.
       ALLOCATE(duxt(nx),STAT=ierr)    ; duxt=0.
+       ALLOCATE(dm(nx),STAT=ierr)      ; dm=diff

     END SUBROUTINE READ_INTRO

@ -330,8 +377,8 @@
     END SUBROUTINE RK4
 !------------------------------------------------------------------------
     SUBROUTINE substep(ri,r1,r2,f)
-      REAL, INTENT(INOUT),DIMENSION(:,:,:) :: ri,r1,r2
-      REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
+      REAL, INTENT(INOUT),DIMENSION(:,:) :: ri,r1,r2
+      REAL, INTENT(OUT),DIMENSION(:,:) :: f

      INTEGER :: i,j,k
      REAL :: at,bt
@ -340,118 +387,66 @@
       IF(istage<5) THEN
        at=a(istage)*dt
        bt=(b(istage)-a(istage))*dt
-        DO k=1,1  ! nz
-        DO j=1,1  ! ny
-        DO i=1,nx
-        r1(i,j,k)=r1(i,j,k)+at*f(i,j,k)
-        r2(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
-        ENDDO
-        ENDDO
-        ENDDO
+        r1=r1+at*f
+        r2=r1+bt*f
       ELSE
        bt=b(istage)*dt
-        DO k=1,1  ! nz
-        DO j=1,1  ! ny
-        DO i=1,nx
-        r1(i,j,k)=r1(i,j,k)+bt*f(i,j,k)
-        ENDDO
-        ENDDO
-        ENDDO
+        r1=r1+bt*f
       ENDIF
     END SUBROUTINE substep
 !------------------------------------------------------------------------
     SUBROUTINE fns(r1,f)
-      REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
-      REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
-      REAL, DIMENSION(3,nx) :: ux,dux,d2ux
+      REAL, INTENT(IN),DIMENSION(:,:) :: r1
+      REAL, INTENT(OUT),DIMENSION(:,:) :: f
+      REAL, DIMENSION(2,nx) :: ux, dux, d2ux
      INTEGER :: i,j,k
-      REAL :: wrate,Ly,Dy
+      REAL :: wrate,Ly,Dy,Lt,Dt

-      DO k=1,1   !nz
 ! x-direction
-       DO j=1,1   !ny
        DO i=1,nx
-         ux(1,i)=u(i,j,k)*r1(i,j,k)                     ! u*Y
-         ux(2,i)=u(i,j,k)                               ! u
-         ux(3,i)=r1(i,j,k)                              ! Y
+         ux(1,i)=r1(i,1)                              ! Y
+         ux(2,i)=r1(i,2)                              ! T
        ENDDO

-        CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
-        CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
+        CALL dfnonp(nx,hx,ux(:,:),dux(:,:),2,1)
+
+        CALL d2fnonp(nx,hx,ux(:,:),d2ux(:,:),2,1)

        DO i=1,nx
-         wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i))))   !wrate
-!         IF ((1.-ux(3,i)).le.c_cut) THEN
-!          wrate=min_wr
-!          IF ((1.-ux(3,i)).gt.c_ref) wrate= &
-!           ((exp(((1.-ux(3,i))-c_cut)*prof_wr)-minf)/(1.-minf)*(refwr-min_wr))+min_wr
-!         ENDIF
-         IF ((1.-ux(3,i)).le.c_ref) THEN
+         wrate=pre*ux(1,i)*exp(-ac/(1.+bc*(ux(2,i))))   !wrate
+
+         IF ((ux(2,i)).le.c_ref) THEN
           wrate=min_wr
-          IF ((1.-ux(3,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(3,i)-c_ref))+ &
+          IF ((ux(2,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(ux(2,i)-c_ref))+ &
          min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
         ENDIF

 ! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x
-         f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + &
-          ux(3,i)*dux(2,i) ) &
-                  + diff*d2ux(1,i) - wrate
+         f(i,1) = - ( u(i)*dux(1,i) ) + diff * d2ux(1,i) - wrate
+
+! -  u*dT/dx + D*d2T/d2x
+         f(i,2) = - ( u(i)*dux(2,i) ) + (lewis * diff) * d2ux(2,i) + wrate
+
         IF (i.eq.nx) THEN
-          Ly=ux(2,nx)*dux(3,nx)                           ! Ly = u*dYr/dx
-          f(nx,1,1)=-wrate
+          Ly=u(nx)*dux(1,nx)                           ! Ly = u*dYr/dx
+          f(nx,1)=-wrate
+
+          Lt=u(nx)*dux(2,nx)                           ! Ly = u*dYr/dx
+          f(nx,2)=wrate
         ENDIF
        ENDDO
-       ENDDO
-!! y-direction
-!       DO i=1,nx
-!        DO j=1,ny
-!         uy(1,j)=v(i,j,k)*r1(i,j,k)                     ! v*Y
-!         uy(2,j)=v(i,j,k)                               ! v
-!         uy(3,j)=r1(i,j,k)                              ! Y
-!        ENDDO
-!        CALL dfyz3(ny,hy,uy(1:3,:),duy)
-!        CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
-!        DO j=1,ny
-!        f(i,j,k)=f(i,j,k) + &
-!        -0.5*( duy(1,j) + uy(2,j)*duy(3,j) + uy(3,j)*duy(2,j) ) &
-!                   + diff*d2uy(j)
-!        ENDDO
-!       ENDDO
-      ENDDO
-!      f(1,1,1)=1.0
-!      f(nx,1,1)=-0.5*( dux(1,nx) + ux(2,nx)*dux(3,nx) + ux(3,nx)*dux(2,nx) )
-!      f(nx,1,1)=0.0
-!! z-direction
-!      DO j=1,ny
-!       DO i=1,nx
-!        DO k=1,nz
-!         uy(1,k)=w(i,j,k)*r1(i,j,k)                     ! w*Y
-!         uy(2,k)=w(i,j,k)                               ! w
-!         uy(3,k)=r1(i,j,k)                              ! Y
-!        ENDDO
-!        CALL dfyz3(ny,hy,uy(1:3,:),duy)
-!        CALL d2fyz1(ny,hy,uy(3:3,:),d2uy)
-!        DO k=1,nz
-!        f(i,j,k)=f(i,j,k) + &
-!        -0.5*( duy(1,k) + uy(2,k)*duy(3,k) + uy(3,k)*duy(2,k) ) &
-!                        + diff*d2uy(k)
-!        ENDDO
-!       ENDDO
-!      ENDDO
-
-!      DO j=1,ny
-!       DO k=1,nz
-!         f(1,j,k)=0.0
-!         f(nx,j,k)=-0.5*(dux(1,nx)+ux(2,nx)*dux(3,nx)+ux(3,nx)*&
-!           dux(2,nx))
-!       ENDDO
-!      ENDDO

 !    Boundary conditionS
-      f(1,1,1)=0.
+      f(1,1)=0.

      Dy=Ly
-      f(nx,1,1)=f(nx,1,1)-Dy
+      f(nx,1)=f(nx,1)-Dy
+
+
+      f(1,2)=0.
+
+      Dt=Lt
+      f(nx,2)=f(nx,2)-Dt

     END SUBROUTINE fns
 !------------------------------------------------------------------------
Author	SHA1	Message	Date
ignis	9f0647865c	solve T	2019-03-28 15:58:09 +09:00
ignis	358d1fea3c	added debug flag to makefile	2019-03-28 14:43:35 +09:00
ignis	c9c16c3f42	added solution variable T	2019-03-27 13:35:36 +09:00
ignis	a61e753b27	makefile fix	2019-03-27 13:34:32 +09:00
ignis	f16cf233c8	Revert "diffusion term contribution from variable diffusivity" This reverts commit `0d771077f6`. Conflicts: code/ysolve.f90	2019-03-27 10:39:43 +09:00
ignis	5b5f958c0f	T^0.76 diffusivity	2019-03-24 13:56:57 +09:00
ignis	805ca37b61	new input parameter sigw which controls sigmoid width	2019-03-24 12:06:25 +09:00
ignis	d74e63524f	input generator take all parameters	2019-03-22 07:04:57 +09:00
Yeongdo Park	c2200912cf	input generator	2019-03-21 04:54:09 -04:00
Yeongdo Park	41fd19c1d7	sigmoid diffusivity, command line options, removing comments	2019-03-21 04:40:27 -04:00
Yeongdo Park	b8ec45ad49	introduce diffusivity model functions	2019-03-19 20:21:08 -04:00
Yeongdo Park	0d771077f6	diffusion term contribution from variable diffusivity	2019-03-19 18:50:59 -04:00
Yeongdo Park	363942207e	diffusivity linear function of c	2019-03-19 03:18:20 -04:00
Yeongdo Park	df1b7fb30a	read input from stdin	2019-03-18 16:26:08 -04:00
Yeongdo Park	b03880aa48	update diffusivity array	2019-03-16 15:58:38 -04:00
Yeongdo Park	fa1d5bcec8	array for variable diffusivity	2019-03-16 15:07:59 -04:00
Yeongdo Park	d6261282a1	added input for variable diffusivity	2019-03-15 15:18:33 -04:00