.gitlab-ci.yml

@@ -16,6 +16,4 @@ test:
   script:
     - cd ${CI_PROJECT_DIR}
     - cd code
-    - make clean && make
+    - make test
-    - cd sample/4pi-IC1
-    - ../../ex

code/makefile

@@ -1,9 +1,5 @@
 flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8  -march=native 
 
-ifdef CHECK
-flags += -fcheck=all
-endif
-
 compiler = gfortran
 
 ex  : test.o ysolve.mod compact.mod

code/sample/4pi-IC1/itape

@@ -16,6 +16,4 @@
      c_ref     0.01  !0.003     ! 
     min_wr     0.        ! 5.0e-14
    prof_wr     1.0       
-      rvis     1.0
-     cstar     1.0
 

code/sample/itape.py

@@ -1,43 +0,0 @@
-#!/bin/python
-
-import os, sys
-
-import argparse
-
-parser = argparse.ArgumentParser()
-parser.add_argument("mu_ratio", help="mu1/mu0, ratio between reaction-diffusion layer mu1 and convection-diffusion layer mu0")
-parser.add_argument("c_star",  help="reference c value seperate reaction-diffusion and convection-diffusion layers")
-args = parser.parse_args()
-
-
-# with open("itape.format") as itapeform:
-    # format = itapeform.read()
-
-itape = '''        nx     512
-       l_0     4.   !3[pi]   
-    int_pr     400
-    tar_lo     0.60
-        dt     0.001 
-        sc     0.75
-       vis     0.020           !2pi(162)ref  !2pi(256)mid  !2pi(256)Max !2p(162)Low
-       pre     2.10E+4       !185.75          !185.75       !316.88       !440.5       !102.69
-        ac     26.7            !22.250          !22.25        !22.56        !22.3        !22.39
-        bc     3.              !7.000           !7.0          !7.0          !7.0         !7.0
-        u0     0.21
-        tf     100.
-     dt_uf     0.4
-      ctmp     0.        !1.0e-14
-     c_cut     0.001 !0.012     ! c < c_cut -> wrate = 0.
-     c_ref     0.01  !0.003     ! 
-    min_wr     0.        ! 5.0e-14
-   prof_wr     1.0       
-      rvis     {}
-     cstar     {}
-'''
-
-print itape.format( args.mu_ratio, args.c_star )
-
-# 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")

code/test.f90

@@ -6,6 +6,5 @@
     PROGRAM test
      USE ysolve
      IMPLICIT NONE
-     CALL parse
      CALL solve
     END PROGRAM test

code/test.sh

@@ -1,4 +1,4 @@
 cd sample/4pi-IC1
 ls
 ../../ex
-diff -y sfield.dat sfield.diff
+diff sfield.dat sfield.diff

code/ysolve.f90

@@ -4,14 +4,12 @@
      PRIVATE
      REAL, PARAMETER :: pi=3.14159265358979323846
      REAL :: hx,dt,vis,sc,diff,pre,ac,bc,tf,t_now,t_uf,dt_uf
-     REAL :: rvis,cstar
      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 :: uxt,duxt
-     REAL, DIMENSION(:), ALLOCATABLE :: dm
 
      INTEGER :: istage
      REAL, DIMENSION(5) :: a=(/ 970286171893./4311952581923.,    &
@@ -25,38 +23,8 @@
                                 2114624349019./3568978502595.,   &
                                 5198255086312./14908931495163. /)
 
-     LOGICAL :: read_itape, read_stdin
+     PUBLIC :: solve
-     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
@@ -116,7 +84,6 @@
       REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
       INTEGER :: i,j,k,nd
       REAL, DIMENSION(1,nx) :: ux,dux,d2ux
-      REAL, DIMENSION(1,nx) :: duy
       REAL, DIMENSION(10,nx) :: sav
         
       sav=0.
@@ -134,15 +101,11 @@
        IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
       ENDDO
           
-      dm = (/ (diff * diffusivity_sigmoid(1.-y1(i,1,1)), i=1,nx) /)
-
       nd=1
       CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
-      CALL dfnonp(nx,hx,dm,duy(1,:),nd,1)
       nd=1 
       CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) 
 
-
       DO i=1,nx
        yr=ux(1,i)
        c=1.-yr
@@ -164,14 +127,14 @@
 
        dely=ABS(dux(1,i))
        sdr=wrate/dely
-       sdd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i))/dely
+       sdd=-diff*d2ux(1,i)/dely
        IF (dely.eq.0.) THEN
         sdr=0.; sdd=0.
        ENDIF
        sd=sdr+sdd
        uu=u(1,1,1)
        onelw=(-dux(1,i))/c
-       dd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i))
+       dd=-diff*d2ux(1,i)
        if (c.eq.0.) onelw=0.
 
        sav(1,i)=sav(1,i)+yr
@@ -286,14 +249,9 @@
      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
@@ -333,18 +291,8 @@
         WRITE(otape,*) cdum,min_wr
         READ(itape,*) cdum,prof_wr
         WRITE(otape,*) cdum,prof_wr
-        READ(itape,*) cdum,rvis
-        WRITE(otape,*) cdum,rvis
-        READ(itape,*) cdum,cstar
-        WRITE(otape,*) cdum,cstar
 
-       DO i = 0,100
+       CLOSE(itape)
-        WRITE(*,*) diffusivity_sigmoid(i / 100.)
-       END DO
-
-       IF (.not.read_stdin) THEN
-         CLOSE(itape)
-       END IF
 
 !       hx=l_0*pi/REAL(nx)
        hx=l_0*pi/REAL(nx-1)
@@ -370,7 +318,6 @@
        ALLOCATE(yf(nx,1,1),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
 
@@ -417,7 +364,6 @@
       REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
       REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
       REAL, DIMENSION(3,nx) :: ux,dux,d2ux
-      REAL, DIMENSION(1,nx) :: duy
       INTEGER :: i,j,k
       REAL :: wrate,Ly,Dy
 
@@ -430,10 +376,7 @@
          ux(3,i)=r1(i,j,k)                              ! Y
         ENDDO
 
-      dm = (/ (diff * diffusivity_sigmoid(1.-r1(i,1,1)), i=1,nx) /)
-
         CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
-        CALL dfnonp(nx,hx,dm,duy(1,:),1,1)
         CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
 
         DO i=1,nx
@@ -452,14 +395,57 @@
 ! -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) ) &
-                  + dm(i)*d2ux(1,i) + duy(1,i)*dux(3,i) - wrate
+                  + diff*d2ux(1,i) - wrate
          IF (i.eq.nx) THEN
           Ly=ux(2,nx)*dux(3,nx)                           ! Ly = u*dYr/dx
           f(nx,1,1)=-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.
@@ -468,21 +454,5 @@
       f(nx,1,1)=f(nx,1,1)-Dy
 
      END SUBROUTINE fns
-!------------------------------------------------------------------------
-     REAL FUNCTION sigmoid(x)
-         REAL, INTENT(IN) :: x
-         sigmoid = 1. / (1 + exp(-x))
-     END FUNCTION sigmoid
-!------------------------------------------------------------------------
-     REAL FUNCTION diffusivity_sigmoid(c)
-         REAL, INTENT(IN) :: c
-         REAL :: a = 100.
-         diffusivity_sigmoid = 1. + (rvis - 1.) * sigmoid(a*(c-cstar))
-     END FUNCTION diffusivity_sigmoid
-!------------------------------------------------------------------------
-     REAL FUNCTION diffusivity_T076(c)
-         REAL, INTENT(IN) :: c
-         diffusivity_T076 = ((1.0 + bc * c) ** 0.76)
-     END FUNCTION diffusivity_T076
 !------------------------------------------------------------------------
     END MODULE ysolve