update diffusivity array

.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

code/makefile

@@ -1,5 +1,9 @@
 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,4 +16,6 @@
      c_ref     0.01  !0.003     ! 
     min_wr     0.        ! 5.0e-14
    prof_wr     1.0       
+      rvis     1.0
+     cstar     1.0
 

code/test.sh

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

code/ysolve.f90

@@ -4,12 +4,14 @@
      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.,    &
@@ -106,6 +108,9 @@
       nd=1 
       CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) 
 
+      dm = diff
+      where (y1(:,1,1)<(1-cstar)) dm = diff * rvis
+
       DO i=1,nx
        yr=ux(1,i)
        c=1.-yr
@@ -127,14 +132,14 @@
 
        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)
        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
@@ -291,6 +296,10 @@
         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
 
        CLOSE(itape)
 
@@ -318,6 +327,7 @@
        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
 
@@ -379,6 +389,9 @@
         CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
         CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
 
+      dm = diff
+      where (r1(:,1,1)<(1.0-cstar)) dm = diff * rvis
+
         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
@@ -395,7 +408,7 @@
 ! -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
+                  + dm(i)*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