diff --git a/myf90/ESBGK/1D_Shock_wave.f90 b/myf90/ESBGK/ESBGK/1D_Shock_wave.f90
similarity index 100%
rename from myf90/ESBGK/1D_Shock_wave.f90
rename to myf90/ESBGK/ESBGK/1D_Shock_wave.f90
diff --git a/myf90/ESBGK/ESBGK MODI2.F90 b/myf90/ESBGK/ESBGK/ESBGK MODI2.F90
similarity index 100%
rename from myf90/ESBGK/ESBGK MODI2.F90
rename to myf90/ESBGK/ESBGK/ESBGK MODI2.F90
diff --git a/myf90/ESBGK/eig.f90 b/myf90/ESBGK/ESBGK/eig.f90
similarity index 100%
rename from myf90/ESBGK/eig.f90
rename to myf90/ESBGK/ESBGK/eig.f90
diff --git a/myf90/ESBGK/makefile b/myf90/ESBGK/ESBGK/makefile
similarity index 100%
rename from myf90/ESBGK/makefile
rename to myf90/ESBGK/ESBGK/makefile
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/IDmodule.f90 b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/IDmodule.f90
new file mode 100644
index 0000000..470f355
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/IDmodule.f90
@@ -0,0 +1,119 @@
+module IDmodule
+ !
+ ! Define Globals here
+ !
+contains
+ subroutine ID_name(name,theta,nx,P_deg,RK_stages,tau,IC_case,fmodel,f_case,method,IDn,IDf)
+ ! Define the inputs and outputs
+ character(len=*), intent(in) :: name
+ integer, intent(in) :: theta,nx,P_deg,RK_stages,IC_case,fmodel,f_case,method
+ real, intent(in) :: tau
+ character(len=*), intent(out) :: IDn,IDf
+ ! Define working variables
+ character(len=20) :: name1,name2,name3,statistic,feq,advec,p_degree,elements,RKs,IC,omega,f
+ ! Define format for ID files names
+ character(len=80) :: format_string,format_string2
+
+ ! Break Name into useful parts "SBBGK1d"
+ name1=name(1:2)
+ name2=name(3:5)
+ name3=name(6:7)
+
+ ! Define statistic
+ select case (theta)
+ case (-1)
+ statistic = "BE"
+ case (0)
+ statistic = "MB"
+ case (1)
+ statistic = "FD"
+ case default
+ print *, 'Not a valid statistic'
+ end select
+
+ ! Define equilibrium distirbution model useq
+ select case (fmodel)
+ case (1) ! UU
+ feq = "-UU"
+ case (2) ! ES
+ feq = "-ES"
+ case default
+ print *, 'Model not available'
+ end select
+
+ ! Define the Method to use
+ select case (method)
+ case (1)
+ advec = "Upwind"
+ P_degree = "1"
+ case (2)
+ advec = "*TVD**"
+ P_degree = "1"
+ case (3)
+ advec = "WENO_3"
+ P_degree = "3"
+ case (4)
+ advec = "WENO_5"
+ P_degree = "5"
+ case (5)
+ advec = "DGWENO"
+ P_degree = char(P_deg+48)
+ case default
+ print *, 'Advection method not available'
+ end select
+
+ ! Define the number of cells to be used,
+ write(elements,"(A1,I3)") 'X',nx
+
+ ! Define the number of RK stages
+ write(Rks,"(I1)") RK_stages
+
+ ! Define the ID of Initial Condition used
+ write(IC,"(I1)") IC_case
+
+ ! Define how to evolve f,
+ select case (f_case)
+ case (1) ! with Collision term-BGK approx
+ write(omega,"(A1,I5)") 'W',ceiling(1.0/tau)
+ case (2) ! no-collison-Euler limit
+ omega = "EL"
+ case default
+ print *, 'case not available'
+ end select
+
+ ! Define Format for file
+ f = ".plt"
+
+ format_string = "(A2,A3,A3,A1,A2,A6,A2,A2,I3,A1,I1,A2,I1,A1,A6,A2,A1,A4)"
+ write(IDf,format_string) trim(name1),trim(feq),trim(name2),'-', &
+ trim(statistic),trim(advec),trim(name3), &
+ '-X',nx,'P',P_deg,'RK',RK_stages,'-',omega,'IC',IC,trim(f)
+ format_string2 = "(A2,A3,A3,A1,A2,A6,A2,A2,I3,A1,I1,A2,I1,A1,A6,A2,A1)"
+ write(IDn,format_string2) trim(name1),trim(feq),trim(name2),' ', &
+ trim(statistic),trim(advec),trim(name3), &
+ ' X',nx,'P',P_deg,'RK',RK_stages,' ',omega,'IC',IC
+ end subroutine ID_name
+
+ subroutine output_names
+ ! I desided to preserve this idea for the future
+ implicit none
+
+ ! character names
+ character(len=12) :: filename
+ character(len=12) :: format_string
+ integer :: counter
+
+ do counter=1, 10
+ if (counter < 10) then
+ format_string = "(A5,I1)"
+ else
+ format_string = "(A5,I2)"
+ endif
+
+ write (filename,format_string) "hello", counter
+ print *, trim(filename)
+ end do
+
+ end subroutine output_names
+
+end module IDmodule
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.cbp b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.cbp
new file mode 100644
index 0000000..9875ee8
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.cbp
@@ -0,0 +1,42 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.depend b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.depend
new file mode 100644
index 0000000..d988863
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.depend
@@ -0,0 +1,5 @@
+# depslib dependency file v1.0
+1374927848 source:/home/manuel/git/aero-shock/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f95
+
+1374928209 source:/home/manuel/git/aero-shock/myf90/ESBGK/TestingDGWENO/TestingDGWENO/IDmodule.f90
+
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.layout b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.layout
new file mode 100644
index 0000000..2f7ca0a
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.layout
@@ -0,0 +1,14 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/com.f90 b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/com.f90
new file mode 100644
index 0000000..2a4ad1b
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/com.f90
@@ -0,0 +1,44 @@
+! the common block
+module data_module
+ implicit none
+ integer nx,md,nv
+
+ !real*8, allocatable :: ph_q_v(:), ph_q_w(:)
+ !real*8 pv_begin,pv_end
+
+ real*8 nh
+ parameter (nx=2560,md=5,nh=nx/2,nv=20)
+
+ real*8 u(0:md,-md:nx+md,1:nv),ueq(0:md,-md:nx+md,1:nv)
+ real*8 v(0:md,0:nx+1,1:nv),hg(0:md,0:nx+1,1:nv)
+ real*8 rki(0:md,-md:nx+md),uki(0:md,-md:nx+md),pki(0:md,-md:nx+md),tki(0:md,-md:nx+md)
+ real*8 eki(0:md,-md:nx+md)
+
+ real x(-md:nx+md),dx(-md:nx+md),bl(0:10),br(0:10),bi(0:md)
+ real*8 ai(0:10,0:10),am1(0:md),xx(200),uu(200),alim(3,3)
+ real*8 aintt(0:10,0:10,0:10)
+ real coef(10,10,10),rco(10,2,10),gauss(10,2),coef9(5,5,6)
+ real*8 rco9(5,5,6)
+ real sigma(2,6),gau(6,2),aii(0:10,0:10),sr
+ real*8 gh(1:nv),w(1:nv),cc(1:nv),vv(1:nv)
+
+ real*8 pi,cfld,dxmin,xm2
+ integer mp,kcmax,indexmax,indexmin
+ integer indexnum,indx1(0:nx), index
+ real*8 t,dt,mm
+
+ !! Input Param !!
+ integer mo !ooa in space
+ integer mt !ooa in time
+ integer kflux !common flux '1=Roe, 2=LF, 3=LLF'
+ real*8 cflc !CFL number, 0.3 for P1, 0.18 for P2 and 0.1 for P3
+ integer ierror !0=initial runs; 1=restart
+ real*8 tprint !input the terminal time
+ integer n !n cells
+ real*8 xmmm !TVB constant M
+ real*8 tau
+ integer phase_quadrature_rule
+ integer init_value !init values 0: 0.5+sin(pi*x) 1:sod problem
+
+ namelist/proj_list/mo,mt,kflux,cflc,tau,phase_quadrature_rule,init_value,ierror,tprint,n,xmmm
+end module data_module
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/configuration.in b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/configuration.in
new file mode 100644
index 0000000..a83bf86
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/configuration.in
@@ -0,0 +1,18 @@
+$Parameters_list
+ name="SBBGK1d",
+ theta=0,
+ nx=100,
+ P_deg=3,
+ RK_stages=3,
+ IC_case=1,
+ fmodel=1,
+ f_case=1,
+ method=1,
+ tau=0.0001,
+ kflux=1,
+ tEnd=0.1,
+ MM=1.0
+/
+!End of list
+!Coded by Manuel Diaz
+NTU, 2013.06.16
diff --git a/myf90/ESBGK/DGburweno.f90 b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/dg_weno_bgk_1d.f90
similarity index 54%
rename from myf90/ESBGK/DGburweno.f90
rename to myf90/ESBGK/TestingDGWENO/TestingDGWENO/dg_weno_bgk_1d.f90
index 17258eb..14a6ca8 100644
--- a/myf90/ESBGK/DGburweno.f90
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/dg_weno_bgk_1d.f90
@@ -1,15 +1,15 @@
-! use discontinuous Galerkin, polynomial truncation for the nonlinear
-! terms, to compute 1D conservation laws
+! Use discontinuous Galerkin, polynomial truncation for the nonlinear
+! terms, to compute 1D Boltzmann-BGK Equation (conservation laws) - July 18, 2013
+! The main modification is to add an outside loop with do is=1, nv (isigma from 1 to nv) enddo
program main
use data_module
integer kcount, lll, kk0, i
real burgex_val
character*32 yc,yc1
-
+
call setup
- open(109,file='wenopoint'//'_p'//char(48+mp)//'_flux'//char(kflux+48)//'_M'//char(mm+48)//'.plt')
-! open(109,file ='output.dat')
+ open(109,file='weno_pt'//'_p'//char(48+mp)//'_flux'//char(kflux+48)//'.plt')
indexmin=1000
indexmax=0
kcount=0
@@ -20,55 +20,91 @@ program main
! call smooth_coef
! begin time iteration
-100 continue
! determine dt
call setdt
-
- if(t.ge.tprint-1.e-10) goto 600
- if(kcount.ge.kcmax) goto 600
- if(kcount/1*1.eq.kcount) then
- ! write(6,*) kcount,' t= ',t
- ! pause
- end if
- ! Runge-Kutta
- call rk
- do lll=1,index
- write(109,*) x(indx1(lll)),t
- enddo
- t=t+dt
- kcount=kcount+1
-
- goto 100
-600 continue
+ do
+ if((t.ge.tprint-1.e-10).or.(kcount.ge.kcmax)) exit
-! print and compute errors
- call limit
- do lll=1,index
- write(109,*) x(indx1(lll)),t
- enddo
+ if(kcount/1*1.eq.kcount) then
+ ! write(6,*) kcount,' t= ',t
+ ! pause
+ end if
- yc='bur'//'_'//char(48+mo)//'_'//char(kflux+48)//'_'//char(mm+48)//'.plt'
- open(1,file=yc)
- do kk0=0,n
- burgex_val = burgex(tprint, x(kk0))
- write(1,123) x(kk0),u(0,kk0),burgex_val
+ do is=1,nv
+ ! Runge-Kutta for all discrete velocity v_{\sigma}
+ call rk(is)
+ do lll=1,index
+ write(109,*) is,x(indx1(lll)),t
+ enddo
+ enddo
+
+ call macrop
+ print *,kcount,t,dt
+
+ t=t+dt
+ kcount=kcount+1
+ call setdt
+ end do
+ call macrop
+
+ ! print and compute errors
+ do is=1,nv
+ call limit(is)
+ do lll=1,index
+ write(109,*) is,x(indx1(lll)),t
+ enddo
enddo
- close(1)
+
+ call tec_output
+
write(*,*) n, indexmin,indexmax
write(*,*) 'points:',index,(x(indx1(i)),i=1,index)
-
+
102 format(i6,1x,3('&',1x, es12.2e2,1x,'&' 1x,f8.2 ,1x)&
- ,'&',1x, i6,1x,'&' 1x,i6 ,1x&
- ,'\\',1x,'\hline')
+ ,'&',1x, i6,1x,'&' 1x,i6 ,1x&
+ ,'\\',1x,'\hline')
103 format(i6,1x,3('&',1x,es12.2E2,1x,'&',1x),'&',1x, i6,1x&
- ,'&' 1x,i6 ,1x,'\\',1x,'\hline')
-123 format(4(1x,f16.6))
+ ,'&' 1x,i6 ,1x,'\\',1x,'\hline')
124 format(i6,1x,3( es12.5e2,1x))
125 format(i6,1x,3( f12.5,1x))
stop
end program main
+
+!!-------------------------------------------
+!!
+!!-------------------------------------------
+subroutine tec_output
+ use data_module
+ real*8 rho,uuu,mx,me,tt,pp
+
+ !write sol
+ open(1,file='sol_cell_avg_'//'_p'//char(48+mo)//'_flux'//char(kflux+48)//'.dat')
+ WRITE(1,*) 'VARIABLES = "x","rho","p","t","u" '
+
+ do i=0,n
+ rho=0.
+ mx=0.
+ me=0.
+ do is=1,nv
+ rho = rho + cc(is) * u(0,i,is)
+ mx = mx + cc(is) * u(0,i,is) * vv(is)
+ me = me + cc(is) * u(0,i,is) * (0.5 * vv(is) * vv(is))
+ enddo
+
+ uuu = mx/rho
+ tt = 4.0d0*me/rho - 2.0d0*uuu**2
+ pp = 0.5d0*rho*tt
+
+ write(1,123) x(i),rho,pp,tt,uuu
+ enddo
+ close(1)
+
+123 format(4(1x,f16.6))
+end subroutine tec_output
+
+
!!-------------------------------------------
!!
!!-------------------------------------------
@@ -81,26 +117,6 @@ subroutine setup
open(unit= 24, file="proj.in", form="FORMATTED", action="READ")
read(24,nml=proj_list)
-
- ! write(6,*) 'input the space order of the scheme (2,3,4)'
- ! read(5,*) mo
- ! write(6,*) 'input the time order of the scheme (3,4)'
- ! read(5,*) mt
- ! write(6,*) '1=Roe, 2=LF, 3=LLF'
- ! read(5,*) kflux
- ! write(6,*) 'input CFL number, 0.3 for P1, 0.18 for P2 and 0.1 for P3'
- ! read(5,*) cflc
- ! ! cflc=0.3
- ! write(6,*) '0=initial runs; 1=restart'
- ! read(5,*) ierror
- ! ! ierror=0
- ! write(6,*) 'input the terminal time'
- ! read(5,*) tprint
- ! ! tprint=0.4
- ! write(6,*) 'input the number of cells in x-direction'
- ! read(5,*) n
- ! write(6,*) 'input TVB constant M'
- ! read(5,*) xmmm
mp=mo-1
if(xmmm.lt.0.10) then
@@ -118,7 +134,7 @@ subroutine setup
endif
kcmax=1000000
-
+
return
end subroutine setup
@@ -127,14 +143,45 @@ end subroutine setup
!!-------------------------------------------
subroutine init
use data_module
+
dimension a(0:10), temp(0:n)
! set up the initial condition
! u0(x)=.5+sin(pi*x)
-
+ gh(1:5) =(/-5.38748089001,-4.60368244955,-3.94476404012,-3.34785456738,-2.78880605843/)
+ gh(6:10) =(/-2.25497400209,-1.73853771212,-1.2340762154,-0.737473728545,- 0.245340708301/)
+ gh(11:15) =(/0.245340708301,0.737473728545,1.2340762154,1.73853771212,2.25497400209/)
+ gh(16:20) =(/2.78880605843,3.34785456738,3.94476404012,4.60368244955,5.38748089001/)
+ w(1:5) =(/0.898591961453,0.704332961176,0.62227869619,0.575262442852,0.544851742366/)
+ w(6:10) = (/0.524080350949,0.509679027117,0.499920871336,0.493843385272,0.490921500667/)
+ w(11:15) =(/0.490921500667,0.493843385272,0.499920871336,0.509679027117,0.524080350949/)
+ w(16:20) =(/0.544851742366,0.575262442852,0.62227869619,0.704332961176,0.898591961453/)
+
+ if (ighq .eq. 1) go to 10
+ do is=1,nv
+ cc(is) =w(is)
+ vv(is) = gh(is)
+ end do
+ go to 30
+10 continue
+ v1 = -10
+ v2 = 10
+ dv = (v2-v1)/(nv-1.)
+ do is =1,nv
+ vv(is) = v1 + dv * (is-1.)
+ end do
+ do is = 2, (nv-1)
+ cc(is) = 64./45. * dv
+ if(mod(is,4) .eq.1) cc(is) = 28./45. * dv
+ if (mod(is,4) .eq.3) cc(is) = 24./45. * dv
+ end do
+ cc(1) = 14./45. * dv
+ cc(nv) = cc(1)
+30 continue
+ pi=4.0*atan(1.0)
xleft= .0
- xright= 2.0
+ xright= 1.0
xlen=xright-xleft
-
+
dxuni=xlen/n
dxmin=1.e10
do i=0,n
@@ -152,26 +199,64 @@ subroutine init
x(n+i)=x(i)+xlen
dx(n+i)=dx(i)
enddo
-
- do i=0,n+1
- do kk=0,mp
- a(kk)=(u0(x(i)+gau(1,1)*dx(i))*fle(kk,gau(1,1))&
- +u0(x(i)+gau(2,1)*dx(i))*fle(kk,gau(2,1)))*gau(1,2)&
- +(u0(x(i)+gau(3,1)*dx(i))*fle(kk,gau(3,1))&
- +u0(x(i)+gau(4,1)*dx(i))*fle(kk,gau(4,1)))*gau(3,2)&
- +(u0(x(i)+gau(5,1)*dx(i))*fle(kk,gau(5,1))&
- +u0(x(i)+gau(6,1)*dx(i))*fle(kk,gau(6,1)))*gau(6,2)
- enddo
- ! take care of the mass matrix
- do kk=0,mp
- if(kk.eq.0) then
- u(kk,i)=ai(kk,kk)*a(kk)
- else
- u(kk,i)=ai(kk,kk)*a(kk)
- endif
+ !Initial condition for shock tube-Density: rho; Mean Velocity: uxm; Pressure: pre; Temperature:tem;
+ ! u(kk,i,is) will be the unknowns, the degree of freedom
+
+ rl=0.125
+ ul=0.
+ pl=0.1
+ rr=1.0
+ ur=0.
+ pr=1.
+
+ do i = 1, n
+ if (x(i) .le. 0.5)then
+ do kk=0,mp
+ rki(kk,i) = rl
+ uki(kk,i) = ul
+ pki(kk,i) = pl
+ tki(kk,i) = 2.0d0*pl/rl
+ end do
+ else
+ do kk=0,mp
+ rki(kk,i) = rr
+ uki(kk,i) = ur
+ pki(kk,i) = pr
+ tki(kk,i) = 2.0d0*pr/rr
+ end do
+ end if
+ do kk=0,mp
+ do is = 1, nv
+ pres = (vv(is)-uki(kk,i))**2/tki(kk,i)
+ u(kk,i,is) = rki(kk,i)*exp(-pres)/sqrt(pi*tki(kk,i))
+ end do
+ end do
+ end do
+
+ do is=1,nv
+ do i=0,n+1
+ do kk=0,mp
+ a(kk)=(u(kk,i,is)*fle(kk,gau(1,1)) + u(kk,i,is)*fle(kk,gau(2,1)))*gau(1,2)&
+ +(u(kk,i,is)*fle(kk,gau(3,1)) + u(kk,i,is)*fle(kk,gau(4,1)))*gau(3,2)&
+ +(u(kk,i,is)*fle(kk,gau(5,1)) + u(kk,i,is)*fle(kk,gau(6,1)))*gau(6,2)
+
+ ! a(kk)=(u0(x(i)+gau(1,1)*dx(i),is)*fle(kk,gau(1,1))&
+ ! +u0(x(i)+gau(2,1)*dx(i),is)*fle(kk,gau(2,1)))*gau(1,2)&
+ ! +(u0(x(i)+gau(3,1)*dx(i),is)*fle(kk,gau(3,1))&
+ ! +u0(x(i)+gau(4,1)*dx(i),is)*fle(kk,gau(4,1)))*gau(3,2)&
+ ! +(u0(x(i)+gau(5,1)*dx(i),is)*fle(kk,gau(5,1))&
+ ! +u0(x(i)+gau(6,1)*dx(i),is)*fle(kk,gau(6,1)))*gau(6,2)
+ enddo
+ ! take care of the mass matrix
+ do kk=0,mp
+ if(kk.eq.0) then
+ u(kk,i,is)=ai(kk,kk)*a(kk)
+ else
+ u(kk,i,is)=ai(kk,kk)*a(kk)
+ endif
+ enddo
enddo
enddo
-
kcount=0
t=0.
@@ -181,23 +266,23 @@ end subroutine init
!!-------------------------------------------
!!
!!-------------------------------------------
-function u0(x00)
- pi=4.0*atan(1.0)
- u0=sin(pi*x00)+0.5
- return
-end function u0
+!function u0(x00)
+! pi=4.0*atan(1.0)
+! u0=sin(pi*x00)+0.5
+! return
+!end function u0
!!-------------------------------------------
!!
!!-------------------------------------------
-subroutine res
+subroutine res(is)
use data_module
-
- dimension h(0:md,0:n+1),flx(0:n+1),a(0:md),b(0:md)&
+
+ dimension h(0:md,0:n+1,nv),flx(0:n+1),a(0:md),b(0:md)&
,un(-1:n+1),up(-1:n+1)
! compute the residue
- ! Burgers' equation
+ ! Boltzmann-BGK equation
! compute the maximum f'(u)
sr=sqrt(5.0)/10.0
@@ -205,69 +290,90 @@ subroutine res
if(kflux.eq.2) then
amax=0.0
do i=1,n
- amax=max(amax,abs(u(0,i)))
+ amax=max(amax,abs(vv(is)))
enddo
end if
-
+
! amax=1.0
! compute the contribution from the cell boundary
do i=0,n+1
- un(i)=eval(u(0,i),mp,0.5)
- up(i-1)=eval(u(0,i),mp,-0.5)
+ un(i)=eval(u(0,i,is),mp,0.5)
+ up(i-1)=eval(u(0,i,is),mp,-0.5)
enddo
do i=0,n
if(kflux.eq.1) then
! Roe flux
- if(un(i)+up(i).gt.0.) then
- flx(i)=flux(un(i))
- else
- flx(i)=flux(up(i))
- end if
+ ! if(un(i)+up(i).gt.0.) then
+ ! flx(i)=vv(is)*un(i)
+ ! flx(i)=flux(un(i))
+ flx(i)=0.5*(vv(is)*up(i) + vv(is)*un(i) - amax*(up(i)-un(i)))
+ !
+ ! flx(i)=vv(is)*up(i)
+ ! flx(i)=flux(up(i))
! flx(i)=un(i)
else if(kflux.eq.2) then
! global LF flux
- flx(i)=0.5*(flux(up(i))+flux(un(i))-amax*(up(i)-un(i)))
+ flx(i)=0.5*(vv(is)*up(i) + vv(is)*un(i) - amax*(up(i)-un(i)))
+ ! flx(i)=0.5*(flux(up(i))+flux(un(i))-amax*(up(i)-un(i)))
else
! local LF flux
amax=max(abs(un(i)),abs(up(i)))
- flx(i)=0.5*(flux(up(i))+flux(un(i))-amax*(up(i)-un(i)))
+ flx(i)=0.5*(vv(is)*up(i) + vv(is)*un(i) - amax*(up(i)-un(i)))
+ ! flx(i)=0.5*(flux(up(i))+flux(un(i))-amax*(up(i)-un(i)))
end if
enddo
do i=1,n
do k=0,mp
- h(k,i)=flx(i-1)*bl(k)-flx(i)*br(k)
+ h(k,i,is)=flx(i-1)*bl(k)-flx(i)*br(k)
enddo
enddo
! compute the contribution from the volume integral by Gauss-Lobatto quadrature
do i=1,n
if (mp.eq.1) then
- h(1,i)=h(1,i)+(flux(up(i-1))+flux(eval(u(0,i),mp,0.0))*4.0+flux(un(i)))/6.0
+ h(1,i,is)=h(1,i,is)+(vv(is)*up(i-1) + vv(is)*eval(u(0,i,is),mp,0.)*4.0 + vv(is)*un(i))/6.0
+ ! h(1,i,is)=h(1,i,is)+(flux(up(i-1))+flux(eval(u(0,i,is),mp,0))*4.0+flux(un(i)))/6.0
else if(mp.eq.2) then
do kk=1,mp
- h(kk,i)=h(kk,i)&
- +(flux(up(i-1))*fled(kk,-0.5)+flux(un(i))*fled(kk,0.5)&
- +(flux(eval(u(0,i),mp,-sr))*fled(kk,-sr)&
- +flux(eval(u(0,i),mp,sr))*fled(kk,sr))*5.0)/12.0
+ h(kk,i,is)=h(kk,i,is)&
+ +( vv(is)*up(i-1)*fled(kk,-0.5) + vv(is)*un(i)*fled(kk,0.5)&
+ +( vv(is)*eval(u(0,i,is),mp,-sr)*fled(kk,-sr)&
+ + vv(is)*eval(u(0,i,is),mp,sr)*fled(kk,sr) )*5.0 )/12.0
+ ! h(kk,i,is)=h(kk,i,is)&
+ ! +(flux(up(i-1))*fled(kk,-0.5)+flux(un(i))*fled(kk,0.5)&
+ ! +(flux(eval(u(0,i,is),mp,-sr))*fled(kk,-sr)&
+ ! +flux(eval(u(0,i,is),mp,sr))*fled(kk,sr))*5.0)/12.0
enddo
else if(mp.eq.3) then
do kk=1,mp
- h(kk,i)=h(kk,i)&
- +((flux(up(i-1))*fled(kk,-0.5)+flux(un(i))*fled(kk,0.5))*9.0&
- +(flux(eval(u(0,i),mp,-sr7))*fled(kk,-sr7)&
- +flux(eval(u(0,i),mp,sr7))*fled(kk,sr7))*49.0&
- +flux(eval(u(0,i),mp,0.0))*fled(kk,0.0)*64.0)/180.0
+ h(kk,i,is)=h(kk,i,is)&
+ +( (vv(is)*up(i-1)*fled(kk,-0.5) + vv(is)*un(i)*fled(kk,0.5) )*9.0&
+ +( vv(is)*eval(u(0,i,is),mp,-sr7)*fled(kk,-sr7)&
+ + vv(is)*eval(u(0,i,is),mp,sr7)*fled(kk,sr7) )*49.0&
+ + vv(is)*eval(u(0,i,is),mp,0.)*fled(kk,0.)*64.0)/180.0
+ ! h(kk,i,is)=h(kk,i,is)&
+ ! +((flux(up(i-1))*fled(kk,-0.5)+flux(un(i))*fled(kk,0.5))*9.0&
+ ! +(flux(eval(u(0,i,is),mp,-sr7))*fled(kk,-sr7)&
+ ! +flux(eval(u(0,i,is),mp,sr7))*fled(kk,sr7))*49.0&
+ ! +flux(eval(u(0,i,is),mp,0))*fled(kk,0.0)*64.0)/180.0
enddo
else
do kk=1,mp
- h(kk,i)=h(kk,i)&
- +(flux(up(i-1))*fled(kk,-.5)+flux(un(i))*fled(kk,.5))*gauss(1,2)&
- +(flux(eval(u(0,i),mp,gauss(3,1)))*fled(kk,gauss(3,1))&
- +flux(eval(u(0,i),mp,gauss(4,1)))*fled(kk,gauss(4,1)))*gauss(3,2)&
- +(flux(eval(u(0,i),mp,gauss(5,1)))*fled(kk,gauss(5,1))&
- +flux(eval(u(0,i),mp,gauss(6,1)))*fled(kk,gauss(6,1)))*gauss(5,2)
+ h(kk,i,is)=h(kk,i,is)&
+ +(vv(is)*up(i-1)*fled(kk,-.5) + vv(is)*un(i)*fled(kk,.5))*gauss(1,2)&
+ +( vv(is)*eval(u(0,i,is),mp,gauss(3,1))*fled(kk,gauss(3,1))&
+ + vv(is)*eval(u(0,i,is),mp,gauss(4,1))*fled(kk,gauss(4,1)) )*gauss(3,2)&
+ +( vv(is)*eval(u(0,i,is),mp,gauss(5,1))*fled(kk,gauss(5,1))&
+ + vv(is)*eval(u(0,i,is),mp,gauss(6,1))*fled(kk,gauss(6,1)) )*gauss(5,2)
+ ! h(kk,i,is)=h(kk,i,is)&
+ ! +(flux(up(i-1))*fled(kk,-.5)+flux(un(i))*fled(kk,.5))*gauss(1,2)&
+ ! +(flux(eval(u(0,i,is),mp,gauss(3,1)))*fled(kk,gauss(3,1))&
+ ! +flux(eval(u(0,i,is),mp,gauss(4,1)))*fled(kk,gauss(4,1)))*gauss(3,2)&
+ ! +(flux(eval(u(0,i,is),mp,gauss(5,1)))*fled(kk,gauss(5,1))&
+ ! +flux(eval(u(0,i,is),mp,gauss(6,1)))*fled(kk,gauss(6,1)))*gauss(5,2)
+
enddo
endif
enddo
@@ -275,13 +381,49 @@ subroutine res
! take care of the mass matrix
do i=1,n
do kk=0,mp
- hg(kk,i)=ai(kk,kk)*h(kk,i)/dx(i)
- ! if(kk.eq.4) hg(kk,i)=0.0
+ hg(kk,i,is)=ai(kk,kk)*h(kk,i,is)/dx(i)
+ ! if(kk.eq.4) hg(kk,i,is)=0.0
enddo
enddo
return
end subroutine res
+!!-------------------------------------------
+!! Find macroscopic properties and equilibrium distribution
+!!-------------------------------------------
+subroutine macrop
+ use data_module
+ do i = 1, n
+ do kk=0,mp
+ skr = 0.
+ sku = 0.
+ ske = 0.
+ do is = 1, nv
+ skr = skr + cc(is) * u(kk,i,is)
+ sku = sku + cc(is) * u(kk,i,is) * vv(is)
+ ske = ske + cc(is) * u(kk,i,is) * (0.5 * vv(is) * vv(is))
+ enddo
+ rki(kk,i) = skr
+ uki(kk,i) = sku/skr
+ eki(kk,i) = ske
+ tki(kk,i) = 4.0d0*eki(kk,i)/rki(kk,i) - 2.0d0*uki(kk,i)*uki(kk,i)
+ pki(kk,i) = 0.5d0*rki(kk,i)*tki(kk,i)
+
+ enddo
+ enddo
+ do is = 1, nv
+ do i = 1, n
+ do kk=0,mp
+ pres = (vv(is)-uki(kk,i))**2/tki(kk,i)
+ ueq(kk,i,is) = rki(kk,i)*exp(-pres)/sqrt(pi*tki(kk,i))
+ enddo
+ enddo
+ enddo
+
+end subroutine macrop
+!!-------------------------------------------
+!!Calculate macroscopic properties and equilibrium distribution
+!!-------------------------------------------
!!-------------------------------------------
!!
@@ -291,8 +433,8 @@ subroutine setdt
! set up dt
amax=0.
- do i=1,n
- amax=max(amax,abs(u(0,i)))
+ do iv=1,nv
+ amax=max(amax,abs(vv(iv)))
enddo
if(mp.eq.3) then
@@ -316,30 +458,29 @@ end subroutine setdt
!!-------------------------------------------
!!
!!-------------------------------------------
-subroutine bc
+subroutine bc(is)
use data_module
-
- ! set up the boundary condition
- ! periodic
+ ! set up the boundary condition keep constant state extrapolation
+ ! periodic
do k=0,mp
do j=0,md
- u(k,-j)=u(k,n-j)
- u(k,n+j+1)=u(k,j+1)
+ u(k,-j,is)=u(k,1,is)
+ u(k,n+j+1,is)=u(k,n,is)
enddo
enddo
-
+
return
end subroutine bc
!!-------------------------------------------
!!
!!-------------------------------------------
-function flux(x)
- ! the function of flux
- flux=x*x/2.0
- return
-end function flux
+!!function flux(x)
+! ! the function of flux
+! flux=x*vv(is)
+! return
+!end function flux
!!-------------------------------------------
!!
@@ -377,7 +518,7 @@ function burgex(te,xe)
if(abs(ay-1.0).lt.1.e-15) y=0.0
if(abs(te).lt.1.e-15) y=sin(pi*xe)
burgex=y*0.5+.250
-
+
return
end function burgex
@@ -447,7 +588,7 @@ subroutine initdata
use data_module
dimension a(0:md),b(0:md),c(0:md),aic(0:20,0:20,0:10),temp(0:n)
-
+
! set up the necessary data before setting the initial condition
do k=0,mp
bl(k)=fle(k,-0.5)
@@ -490,7 +631,7 @@ subroutine initdata
aic(3,3,4)=2800.0
aic(3,4,4)=0.0
aic(4,4,4)=44100.0
-
+
do i0=1,4
do i01=0,i0
do i02=i01+1,i0
@@ -563,7 +704,7 @@ subroutine initdata
gau(4,1)= sqrt(147.-42.*sqrt(7.))/42.0
gau(5,1)=-sqrt(147.+42.*sqrt(7.))/42.0
gau(6,1)= sqrt(147.+42.*sqrt(7.))/42.0
-
+
! coefficients of 10th order gau-Lobatto quadrature
gau(1,2)=1.0/30.0
gau(2,2)=gau(1,2)
@@ -580,7 +721,7 @@ subroutine initdata
gauss(1,2)=1.0/6.0
gauss(2,2)=gauss(1,2)
gauss(3,2)=2.0/3.0
-
+
endif
if(mp.eq.2) then
! the points of 6th order Gauss-Lobatto quadrature
@@ -617,7 +758,7 @@ subroutine initdata
gauss(6,1)= sqrt(147.-42.*sqrt(7.))/42.0
gauss(3,1)=-sqrt(147.+42.*sqrt(7.))/42.0
gauss(4,1)= sqrt(147.+42.*sqrt(7.))/42.0
-
+
! coefficients of 10th order Gauss-Lobatto quadrature
gauss(1,2)=1.0/30.0
gauss(2,2)=gauss(1,2)
@@ -636,10 +777,10 @@ end subroutine initdata
!!-------------------------------------------
!!
!!-------------------------------------------
-subroutine limit
+subroutine limit(is)
use data_module
- dimension b(0:md),du(-1:n+1),am(10),temp(10)
+ dimension b(0:md),du(-1:n+1,nv),am(10),temp(10)
! xmmm=-2.0/3.0
! limit the solution
index=0
@@ -648,21 +789,21 @@ subroutine limit
return
else
! compute the first order difference of the mean for limiting
+
do i=-1,n
- du(i)=u(0,i+1)-u(0,i)
+ du(i,is)=u(0,i+1,is)-u(0,i,is)
enddo
-
do i=0,n
xmm=xmmm*dx(i)**2
! compute the cell boundary values:
- uleft=eval(u(0,i),mp,0.5)
- uright=eval(u(0,i),mp,-0.5)
+ uleft=eval(u(0,i,is),mp,0.5)
+ uright=eval(u(0,i,is),mp,-0.5)
! limit the values uleft and uright
- am(1)=uleft-u(0,i)
- am(2)=du(i-1)
- am(3)=du(i)
+ am(1)=uleft-u(0,i,is)
+ am(2)=du(i-1,is)
+ am(3)=du(i,is)
sg=sign(1.,am(1))
amc=am(1)
if(abs(amc).le.xmm) goto 31
@@ -681,9 +822,9 @@ subroutine limit
31 continue
! if(abs(amc-am(1)).gt.1e-6) goto 100
! if(amc.ne.am(1)) goto 100
- am(1)=u(0,i)-uright
- am(2)=du(i-1)
- am(3)=du(i)
+ am(1)=u(0,i,is)-uright
+ am(2)=du(i-1,is)
+ am(3)=du(i,is)
amc=am(1)
if(abs(amc).le.xmm) cycle
do j=2,3
@@ -700,13 +841,13 @@ subroutine limit
cycle
! 100 if(abs(amc-am(1)).gt.1e-6) then
330 index=index+1
- call wenorecon(i)
+ call wenorecon(i,is)
indexmin=min(indexmin,index)
indexmax=max(indexmax,index)
indx1(index)=i
end do
end if
-
+
if(index.gt.n/20) then
! write(*,*) index
endif
@@ -716,124 +857,122 @@ end subroutine limit
!!-------------------------------------------
!!
!!-------------------------------------------
-subroutine rk
+subroutine rk(is)
use data_module
dimension temp1(0:mp,n,3)
! third order Runge-Kutta
-
! set up boundary condition for u
if(mt.eq.3) then
- call bc
+ call bc(is)
! limit the solution
- call limit
- call bc
+ call limit(is)
+ call bc(is)
do i=1,n
- ! uold(i)=u(0,i)
+ ! uold(i)=u(0,i,is)
do k=0,mp
- v(k,i)=u(k,i)
+ v(k,i,is)=u(k,i,is)
enddo
enddo
! compute the residue
- call res
+ call res(is)
! call wenorecon
do i=1,n
do k=0,mp
- u(k,i)=v(k,i)+dt*hg(k,i)
- temp1(k,i,1)=u(k,i)
+ u(k,i,is)=v(k,i,is)+dt*hg(k,i,is)
+ temp1(k,i,1)=u(k,i,is)
enddo
enddo
! call bc
! call limit
- call bc
- call res
+ call bc(is)
+ call res(is)
! call wenorecon
do i=1,n
do k=0,mp
- u(k,i)=0.75*v(k,i)+0.25*(temp1(k,i,1)+dt*hg(k,i))
- temp1(k,i,2)=u(k,i)
+ u(k,i,is)=0.75*v(k,i,is)+0.25*(temp1(k,i,1)+dt*hg(k,i,is))
+ temp1(k,i,2)=u(k,i,is)
enddo
enddo
! call bc
! call limit
- call bc
- call res
+ call bc(is)
+ call res(is)
! call wenorecon
do i=1,n
do k=0,mp
- u(k,i)=(v(k,i)+2.0*(temp1(k,i,2)+dt*hg(k,i)))/3.0
+ u(k,i,is)=(v(k,i,is)+2.0*(temp1(k,i,2)+dt*hg(k,i,is)))/3.0
enddo
enddo
else
! set up boundary condition for u
- call bc
+ call bc(is)
! limit the solution
- call limit
- call bc
+ call limit(is)
+ call bc(is)
do i=1,n
! uold(i)=u(0,i)
do k=0,mp
- v(k,i)=u(k,i)
+ v(k,i,is)=u(k,i,is)
enddo
enddo
! compute the residue
- call res
+ call res(is)
! call wenorecon
do i=1,n
do k=0,mp
- u(k,i)=v(k,i)+0.5*dt*hg(k,i)
- temp1(k,i,1)=u(k,i)
+ u(k,i,is)=v(k,i,is)+0.5*dt*hg(k,i,is)
+ temp1(k,i,1)=u(k,i,is)
enddo
enddo
- call bc
- call limit
- call bc
- call res
+ call bc(is)
+ call limit(is)
+ call bc(is)
+ call res(is)
! call wenorecon
do i=1,n
do k=0,mp
- u(k,i)=v(k,i)+0.5*dt*hg(k,i)
- temp1(k,i,2)=u(k,i)
+ u(k,i,is)=v(k,i,is)+0.5*dt*hg(k,i,is)
+ temp1(k,i,2)=u(k,i,is)
enddo
enddo
- call bc
- call limit
- call bc
- call res
+ call bc(is)
+ call limit(is)
+ call bc(is)
+ call res(is)
! call wenorecon
do i=1,n
do k=0,mp
- u(k,i)=v(k,i)+dt*hg(k,i)
- temp1(k,i,3)=u(k,i)
+ u(k,i,is)=v(k,i,is)+dt*hg(k,i,is)
+ temp1(k,i,3)=u(k,i,is)
enddo
enddo
- call bc
- call limit
- call bc
- call res
+ call bc(is)
+ call limit(is)
+ call bc(is)
+ call res(is)
! call wenorecon
do i=1,n
do k=0,mp
- u(k,i)=(-v(k,i)+temp1(k,i,1)+2*temp1(k,i,2)&
- +temp1(k,i,3)+0.5*dt*hg(k,i))/3.
+ u(k,i,is)=(-v(k,i,is)+temp1(k,i,1)+2*temp1(k,i,2)&
+ +temp1(k,i,3)+0.5*dt*hg(k,i,is))/3.
enddo
enddo
endif
-
return
end subroutine rk
@@ -869,149 +1008,149 @@ end function qll
!!-------------------------------------------
!! polynomial function of reconstruction with cells: kk+k1---kk+k2
!!-------------------------------------------
-function pll(x0,k1,k2,i)
+function pll(x0,k1,k2,i,is)
use data_module
! dimension a(-10:10)
pll=0.0
do l=k1,k2
- pll=pll+qll(x0,k1,k2,l,i)*u(0,l+i)
+ pll=pll+qll(x0,k1,k2,l,i)*u(0,l+i,is)
enddo
return
end function pll
!!-------------------------------------------
-!! set the coefficent of function pll, and compute the smooth indicators
+!! set the coefficent of function pll, and compute the smooth indicators (is indicates isigma)
!!-------------------------------------------
-subroutine smooth(i,s)
+subroutine smooth(i,is,s)
use data_module
- dimension s(0:10),aa(0:5,2)
+ dimension s(0:10,nv),aa(0:5,2)
if(mp.eq.1) then
- aa(1,1)=(pll(gau(1,1),-1,0,i)*gau(1,1)&
- +pll(gau(2,1),-1,0,i)*gau(2,1))*gau(1,2)&
- +(pll(gau(3,1),-1,0,i)*gau(3,1)&
- +pll(gau(4,1),-1,0,i)*gau(4,1))*gau(3,2)&
- +(pll(gau(5,1),-1,0,i)*gau(5,1)&
- +pll(gau(6,1),-1,0,i)*gau(6,1))*gau(5,2)
- s(1)=aa(1,1)**2
- aa(1,1)=(pll(gau(1,1),0,1,i)*gau(1,1)&
- +pll(gau(2,1),0,1,i)*gau(2,1))*gau(1,2)&
- +(pll(gau(3,1),0,1,i)*gau(3,1)&
- +pll(gau(4,1),0,1,i)*gau(4,1))*gau(3,2)&
- +(pll(gau(5,1),0,1,i)*gau(5,1)&
- +pll(gau(6,1),0,1,i)*gau(6,1))*gau(5,2)
- s(2)=aa(1,1)**2
+ aa(1,1)=(pll(gau(1,1),-1,0,i,is)*gau(1,1)&
+ +pll(gau(2,1),-1,0,i,is)*gau(2,1))*gau(1,2)&
+ +(pll(gau(3,1),-1,0,i,is)*gau(3,1)&
+ +pll(gau(4,1),-1,0,i,is)*gau(4,1))*gau(3,2)&
+ +(pll(gau(5,1),-1,0,i,is)*gau(5,1)&
+ +pll(gau(6,1),-1,0,i,is)*gau(6,1))*gau(5,2)
+ s(1,is)=aa(1,1)**2
+ aa(1,1)=(pll(gau(1,1),0,1,i,is)*gau(1,1)&
+ +pll(gau(2,1),0,1,i,is)*gau(2,1))*gau(1,2)&
+ +(pll(gau(3,1),0,1,i,is)*gau(3,1)&
+ +pll(gau(4,1),0,1,i,is)*gau(4,1))*gau(3,2)&
+ +(pll(gau(5,1),0,1,i,is)*gau(5,1)&
+ +pll(gau(6,1),0,1,i,is)*gau(6,1))*gau(5,2)
+ s(2,is)=aa(1,1)**2
elseif(mp.eq.2) then
do l=0,2
- aa(0,1)=(pll(gau(1,1),l-2,l,i)&
- +pll(gau(2,1),l-2,l,i))*gau(1,2)&
- +(pll(gau(3,1),l-2,l,i)&
- +pll(gau(4,1),l-2,l,i))*gau(3,2)&
- +(pll(gau(5,1),l-2,l,i)&
- +pll(gau(6,1),l-2,l,i))*gau(5,2)
- aa(1,1)=(pll(gau(1,1),l-2,l,i)*gau(1,1)&
- +pll(gau(2,1),l-2,l,i)*gau(2,1))*gau(1,2)&
- +(pll(gau(3,1),l-2,l,i)*gau(3,1)&
- +pll(gau(4,1),l-2,l,i)*gau(4,1))*gau(3,2)&
- +(pll(gau(5,1),l-2,l,i)*gau(5,1)&
- +pll(gau(6,1),l-2,l,i)*gau(6,1))*gau(5,2)
- aa(2,1)=(pll(gau(1,1),l-2,l,i)*gau(1,1)**2&
- +pll(gau(2,1),l-2,l,i)*gau(2,1)**2)*gau(1,2)&
- +(pll(gau(3,1),l-2,l,i)*gau(3,1)**2&
- +pll(gau(4,1),l-2,l,i)*gau(4,1)**2)*gau(3,2)&
- +(pll(gau(5,1),l-2,l,i)*gau(5,1)**2&
- +pll(gau(6,1),l-2,l,i)*gau(6,1)**2)*gau(5,2)
+ aa(0,1)=(pll(gau(1,1),l-2,l,i,is)&
+ +pll(gau(2,1),l-2,l,i,is))*gau(1,2)&
+ +(pll(gau(3,1),l-2,l,i,is)&
+ +pll(gau(4,1),l-2,l,i,is))*gau(3,2)&
+ +(pll(gau(5,1),l-2,l,i,is)&
+ +pll(gau(6,1),l-2,l,i,is))*gau(5,2)
+ aa(1,1)=(pll(gau(1,1),l-2,l,i,is)*gau(1,1)&
+ +pll(gau(2,1),l-2,l,i,is)*gau(2,1))*gau(1,2)&
+ +(pll(gau(3,1),l-2,l,i,is)*gau(3,1)&
+ +pll(gau(4,1),l-2,l,i,is)*gau(4,1))*gau(3,2)&
+ +(pll(gau(5,1),l-2,l,i,is)*gau(5,1)&
+ +pll(gau(6,1),l-2,l,i,is)*gau(6,1))*gau(5,2)
+ aa(2,1)=(pll(gau(1,1),l-2,l,i,is)*gau(1,1)**2&
+ +pll(gau(2,1),l-2,l,i,is)*gau(2,1)**2)*gau(1,2)&
+ +(pll(gau(3,1),l-2,l,i,is)*gau(3,1)**2&
+ +pll(gau(4,1),l-2,l,i,is)*gau(4,1)**2)*gau(3,2)&
+ +(pll(gau(5,1),l-2,l,i,is)*gau(5,1)**2&
+ +pll(gau(6,1),l-2,l,i,is)*gau(6,1)**2)*gau(5,2)
do ll=1,2
aa(ll,2)=aa(0,1)*aii(ll,0)+aa(1,1)*aii(ll,1)+aa(2,1)*aii(ll,2)
enddo
- s(l+1)=aa(1,2)**2+13.0/3.0*aa(2,2)**2
+ s(l+1,is)=aa(1,2)**2+13.0/3.0*aa(2,2)**2
enddo
elseif(mp.eq.3) then
do l=0,3
- aa(0,1)=(pll(gau(1,1),l-3,l,i)&
- +pll(gau(2,1),l-3,l,i))*gau(1,2)&
- +(pll(gau(3,1),l-3,l,i)&
- +pll(gau(4,1),l-3,l,i))*gau(3,2)&
- +(pll(gau(5,1),l-3,l,i)&
- +pll(gau(6,1),l-3,l,i))*gau(5,2)
- aa(1,1)=(pll(gau(1,1),l-3,l,i)*gau(1,1)&
- +pll(gau(2,1),l-3,l,i)*gau(2,1))*gau(1,2)&
- +(pll(gau(3,1),l-3,l,i)*gau(3,1)&
- +pll(gau(4,1),l-3,l,i)*gau(4,1))*gau(3,2)&
- +(pll(gau(5,1),l-3,l,i)*gau(5,1)&
- +pll(gau(6,1),l-3,l,i)*gau(6,1))*gau(5,2)
- aa(2,1)=(pll(gau(1,1),l-3,l,i)*gau(1,1)**2&
- +pll(gau(2,1),l-3,l,i)*gau(2,1)**2)*gau(1,2)&
- +(pll(gau(3,1),l-3,l,i)*gau(3,1)**2&
- +pll(gau(4,1),l-3,l,i)*gau(4,1)**2)*gau(3,2)&
- +(pll(gau(5,1),l-3,l,i)*gau(5,1)**2&
- +pll(gau(6,1),l-3,l,i)*gau(6,1)**2)*gau(5,2)
- aa(3,1)=(pll(gau(1,1),l-3,l,i)*gau(1,1)**3&
- +pll(gau(2,1),l-3,l,i)*gau(2,1)**3)*gau(1,2)&
- +(pll(gau(3,1),l-3,l,i)*gau(3,1)**3&
- +pll(gau(4,1),l-3,l,i)*gau(4,1)**3)*gau(3,2)&
- +(pll(gau(5,1),l-3,l,i)*gau(5,1)**3&
- +pll(gau(6,1),l-3,l,i)*gau(6,1)**3)*gau(5,2)
+ aa(0,1)=(pll(gau(1,1),l-3,l,i,is)&
+ +pll(gau(2,1),l-3,l,i,is))*gau(1,2)&
+ +(pll(gau(3,1),l-3,l,i,is)&
+ +pll(gau(4,1),l-3,l,i,is))*gau(3,2)&
+ +(pll(gau(5,1),l-3,l,i,is)&
+ +pll(gau(6,1),l-3,l,i,is))*gau(5,2)
+ aa(1,1)=(pll(gau(1,1),l-3,l,i,is)*gau(1,1)&
+ +pll(gau(2,1),l-3,l,i,is)*gau(2,1))*gau(1,2)&
+ +(pll(gau(3,1),l-3,l,i,is)*gau(3,1)&
+ +pll(gau(4,1),l-3,l,i,is)*gau(4,1))*gau(3,2)&
+ +(pll(gau(5,1),l-3,l,i,is)*gau(5,1)&
+ +pll(gau(6,1),l-3,l,i,is)*gau(6,1))*gau(5,2)
+ aa(2,1)=(pll(gau(1,1),l-3,l,i,is)*gau(1,1)**2&
+ +pll(gau(2,1),l-3,l,i,is)*gau(2,1)**2)*gau(1,2)&
+ +(pll(gau(3,1),l-3,l,i,is)*gau(3,1)**2&
+ +pll(gau(4,1),l-3,l,i,is)*gau(4,1)**2)*gau(3,2)&
+ +(pll(gau(5,1),l-3,l,i,is)*gau(5,1)**2&
+ +pll(gau(6,1),l-3,l,i,is)*gau(6,1)**2)*gau(5,2)
+ aa(3,1)=(pll(gau(1,1),l-3,l,i,is)*gau(1,1)**3&
+ +pll(gau(2,1),l-3,l,i,is)*gau(2,1)**3)*gau(1,2)&
+ +(pll(gau(3,1),l-3,l,i,is)*gau(3,1)**3&
+ +pll(gau(4,1),l-3,l,i,is)*gau(4,1)**3)*gau(3,2)&
+ +(pll(gau(5,1),l-3,l,i,is)*gau(5,1)**3&
+ +pll(gau(6,1),l-3,l,i,is)*gau(6,1)**3)*gau(5,2)
do ll=1,3
aa(ll,2)=aa(0,1)*aii(ll,0)+aa(1,1)*aii(ll,1)&
+aa(2,1)*aii(ll,2)+aa(3,1)*aii(ll,3)
enddo
- s(l+1)=3129./80.*aa(3,2)**2+0.5*aa(1,2)*aa(3,2)&
+ s(l+1,is)=3129./80.*aa(3,2)**2+0.5*aa(1,2)*aa(3,2)&
+13./3.*aa(2,2)**2+aa(1,2)**2
enddo
else
do l=0,4
- aa(0,1)=(pll(gau(1,1),l-4,l,i)&
- +pll(gau(2,1),l-4,l,i))*gau(1,2)&
- +(pll(gau(3,1),l-4,l,i)&
- +pll(gau(4,1),l-4,l,i))*gau(3,2)&
- +(pll(gau(5,1),l-4,l,i)&
- +pll(gau(6,1),l-4,l,i))*gau(5,2)
- aa(1,1)=(pll(gau(1,1),l-4,l,i)*gau(1,1)&
- +pll(gau(2,1),l-4,l,i)*gau(2,1))*gau(1,2)&
- +(pll(gau(3,1),l-4,l,i)*gau(3,1)&
- +pll(gau(4,1),l-4,l,i)*gau(4,1))*gau(3,2)&
- +(pll(gau(5,1),l-4,l,i)*gau(5,1)&
- +pll(gau(6,1),l-4,l,i)*gau(6,1))*gau(5,2)
- aa(2,1)=(pll(gau(1,1),l-4,l,i)*gau(1,1)**2&
- +pll(gau(2,1),l-4,l,i)*gau(2,1)**2)*gau(1,2)&
- +(pll(gau(3,1),l-4,l,i)*gau(3,1)**2&
- +pll(gau(4,1),l-4,l,i)*gau(4,1)**2)*gau(3,2)&
- +(pll(gau(5,1),l-4,l,i)*gau(5,1)**2&
- +pll(gau(6,1),l-4,l,i)*gau(6,1)**2)*gau(5,2)
- aa(3,1)=(pll(gau(1,1),l-4,l,i)*gau(1,1)**3&
- +pll(gau(2,1),l-4,l,i)*gau(2,1)**3)*gau(1,2)&
- +(pll(gau(3,1),l-4,l,i)*gau(3,1)**3&
- +pll(gau(4,1),l-4,l,i)*gau(4,1)**3)*gau(3,2)&
- +(pll(gau(5,1),l-4,l,i)*gau(5,1)**3&
- +pll(gau(6,1),l-4,l,i)*gau(6,1)**3)*gau(5,2)
- aa(4,1)=(pll(gau(1,1),l-4,l,i)*gau(1,1)**4&
- +pll(gau(2,1),l-4,l,i)*gau(2,1)**4)*gau(1,2)&
- +(pll(gau(3,1),l-4,l,i)*gau(3,1)**4&
- +pll(gau(4,1),l-4,l,i)*gau(4,1)**4)*gau(3,2)&
- +(pll(gau(5,1),l-4,l,i)*gau(5,1)**4&
- +pll(gau(6,1),l-4,l,i)*gau(6,1)**4)*gau(5,2)
+ aa(0,1)=(pll(gau(1,1),l-4,l,i,is)&
+ +pll(gau(2,1),l-4,l,i,is))*gau(1,2)&
+ +(pll(gau(3,1),l-4,l,i,is)&
+ +pll(gau(4,1),l-4,l,i,is))*gau(3,2)&
+ +(pll(gau(5,1),l-4,l,i,is)&
+ +pll(gau(6,1),l-4,l,i,is))*gau(5,2)
+ aa(1,1)=(pll(gau(1,1),l-4,l,i,is)*gau(1,1)&
+ +pll(gau(2,1),l-4,l,i,is)*gau(2,1))*gau(1,2)&
+ +(pll(gau(3,1),l-4,l,i,is)*gau(3,1)&
+ +pll(gau(4,1),l-4,l,i,is)*gau(4,1))*gau(3,2)&
+ +(pll(gau(5,1),l-4,l,i,is)*gau(5,1)&
+ +pll(gau(6,1),l-4,l,i,is)*gau(6,1))*gau(5,2)
+ aa(2,1)=(pll(gau(1,1),l-4,l,i,is)*gau(1,1)**2&
+ +pll(gau(2,1),l-4,l,i,is)*gau(2,1)**2)*gau(1,2)&
+ +(pll(gau(3,1),l-4,l,i,is)*gau(3,1)**2&
+ +pll(gau(4,1),l-4,l,i,is)*gau(4,1)**2)*gau(3,2)&
+ +(pll(gau(5,1),l-4,l,i,is)*gau(5,1)**2&
+ +pll(gau(6,1),l-4,l,i,is)*gau(6,1)**2)*gau(5,2)
+ aa(3,1)=(pll(gau(1,1),l-4,l,i,is)*gau(1,1)**3&
+ +pll(gau(2,1),l-4,l,i,is)*gau(2,1)**3)*gau(1,2)&
+ +(pll(gau(3,1),l-4,l,i,is)*gau(3,1)**3&
+ +pll(gau(4,1),l-4,l,i,is)*gau(4,1)**3)*gau(3,2)&
+ +(pll(gau(5,1),l-4,l,i,is)*gau(5,1)**3&
+ +pll(gau(6,1),l-4,l,i,is)*gau(6,1)**3)*gau(5,2)
+ aa(4,1)=(pll(gau(1,1),l-4,l,i,is)*gau(1,1)**4&
+ +pll(gau(2,1),l-4,l,i,is)*gau(2,1)**4)*gau(1,2)&
+ +(pll(gau(3,1),l-4,l,i,is)*gau(3,1)**4&
+ +pll(gau(4,1),l-4,l,i,is)*gau(4,1)**4)*gau(3,2)&
+ +(pll(gau(5,1),l-4,l,i,is)*gau(5,1)**4&
+ +pll(gau(6,1),l-4,l,i,is)*gau(6,1)**4)*gau(5,2)
do ll=1,4
aa(ll,2)=aa(0,1)*aii(ll,0)+aa(1,1)*aii(ll,1)&
+aa(2,1)*aii(ll,2)+aa(3,1)*aii(ll,3)+aa(4,1)*aii(ll,4)
enddo
- s(l+1)=87617./140.*aa(4,2)**2+4.25*aa(2,2)*aa(4,2)&
+ s(l+1,is)=87617./140.*aa(4,2)**2+4.25*aa(2,2)*aa(4,2)&
+3129./80.*aa(3,2)**2+0.5*aa(1,2)*aa(3,2)&
+13./3.*aa(2,2)**2+aa(1,2)**2
enddo
endif
return
end subroutine smooth
-
+
!!-------------------------------------------
!! reconstruction polynomial on trouble cell by WENO
!!-------------------------------------------
-subroutine wenorecon(i)
+subroutine wenorecon(i,is)
use data_module
- dimension s(0:10),ww(6),pp1(6),a(-10:10),temp(6)
+ dimension s(0:10,nv),ww(6),pp1(6),a(-10:10),temp(6)
weps=1.0e-6
sr=sqrt(5.0)
- call smooth(i,s)
+ call smooth(i,is,s)
if(mp.eq.1) then
do l=1,2
ko=mp
@@ -1032,30 +1171,30 @@ subroutine wenorecon(i)
rco9(kk+1,l,1)=(qll(x1,-ko,ko,kk,i)-ttc)/coef9(ko+1,kk+1,l)
enddo
enddo
-
+
do l=1,2
temp(l)=0.0
sum=0.0
do j=1,ko+1
- sum=rco9(j,l,1)/(weps+s(j))**2+sum
+ sum=rco9(j,l,1)/(weps+s(j,is))**2+sum
enddo
do j=1,ko+1
- ww(j)=rco9(j,l,1)/(weps+s(j))**2/sum
+ ww(j)=rco9(j,l,1)/(weps+s(j,is))**2/sum
enddo
-
+
do ll=1,ko+1
pp1(ll)=0.0
do j=1,ko+1
- pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll)
+ pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll,is)
enddo
enddo
do ll=1,ko+1
temp(l)=temp(l)+ww(ll)*pp1(ll)
enddo
-
+
enddo
- u(1,i)=(temp(2)-temp(1))
-
+ u(1,i,is)=(temp(2)-temp(1))
+
elseif(mp.eq.2) then
do l=1,4
ko=mp
@@ -1075,36 +1214,36 @@ subroutine wenorecon(i)
enddo
rco9(kk+1,l,1)=(qll(x1,-ko,ko,kk,i)-ttc)/coef9(ko+1,kk+1,l)
enddo
-! write(*,*) x1
-! write(*,*) coef9,rco9
-! pause
+ ! write(*,*) x1
+ ! write(*,*) coef9,rco9
+ ! pause
enddo
do l=1,4
temp(l)=0.0
sum=0.0
do j=1,ko+1
- sum=rco9(j,l,1)/(weps+s(j))**2+sum
+ sum=rco9(j,l,1)/(weps+s(j,is))**2+sum
enddo
do j=1,ko+1
- ww(j)=rco9(j,l,1)/(weps+s(j))**2/sum
+ ww(j)=rco9(j,l,1)/(weps+s(j,is))**2/sum
enddo
-
+
do ll=1,ko+1
pp1(ll)=0.0
do j=1,ko+1
- pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll)
+ pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll,is)
enddo
enddo
do ll=1,ko+1
temp(l)=temp(l)+ww(ll)*pp1(ll)
enddo
-
+
enddo
-
- u(1,i)=(temp(2)-temp(1)+sr*(temp(4)-temp(3)))/2.0
- u(2,i)=(temp(2)+temp(1)-(temp(4)+temp(3)))*2.5
-
+
+ u(1,i,is)=(temp(2)-temp(1)+sr*(temp(4)-temp(3)))/2.0
+ u(2,i,is)=(temp(2)+temp(1)-(temp(4)+temp(3)))*2.5
+
elseif(mp.eq.3) then
do l=1,4
ko=mp
@@ -1124,25 +1263,25 @@ subroutine wenorecon(i)
enddo
rco9(kk+1,l,1)=(qll(x1,-ko,ko,kk,i)-ttc)/coef9(ko+1,kk+1,l)
enddo
-! write(*,*) x1
-! write(*,*) coef9,rco9
-! pause
+ ! write(*,*) x1
+ ! write(*,*) coef9,rco9
+ ! pause
enddo
-
+
do l=1,4
temp(l)=0.0
sum=0.0
do j=1,ko+1
- sum=rco9(j,l,1)/(weps+s(j))**2+sum
+ sum=rco9(j,l,1)/(weps+s(j,is))**2+sum
enddo
do j=1,ko+1
- ww(j)=rco9(j,l,1)/(weps+s(j))**2/sum
+ ww(j)=rco9(j,l,1)/(weps+s(j,is))**2/sum
enddo
do ll=1,ko+1
pp1(ll)=0.0
do j=1,ko+1
- pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll)
+ pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll,is)
enddo
enddo
do ll=1,ko+1
@@ -1151,20 +1290,20 @@ subroutine wenorecon(i)
enddo
- temp(5)=(180.0*u(0,i)-9.0*(temp(1)+temp(2))&
+ temp(5)=(180.0*u(0,i,is)-9.0*(temp(1)+temp(2))&
-(temp(3)+temp(4))*49.0)/64.0
- ! reconstruction of u(1,i), u(2,i),u(3,i)
- u(1,i)=((fle(1,gauss(1,1))*temp(1)&
+ ! reconstruction of u(1,i,is), u(2,i,is),u(3,i,is)
+ u(1,i,is)=((fle(1,gauss(1,1))*temp(1)&
+fle(1,gauss(2,1))*temp(2))*gauss(1,2)&
+(fle(1,gauss(3,1))*temp(3)&
+fle(1,gauss(4,1))*temp(4))*gauss(3,2)&
+(fle(1,gauss(5,1))*temp(5))*gauss(5,2))*12.0
- u(2,i)=((fle(2,gauss(1,1))*temp(1)&
+ u(2,i,is)=((fle(2,gauss(1,1))*temp(1)&
+fle(2,gauss(2,1))*temp(2))*gauss(1,2)&
+(fle(2,gauss(3,1))*temp(3)&
+fle(2,gauss(4,1))*temp(4))*gauss(3,2)&
+(fle(2,gauss(5,1))*temp(5))*gauss(5,2))*180.0
- u(3,i)=((fle(3,gauss(1,1))*temp(1)&
+ u(3,i,is)=((fle(3,gauss(1,1))*temp(1)&
+fle(3,gauss(2,1))*temp(2))*gauss(1,2)&
+(fle(3,gauss(3,1))*temp(3)&
+fle(3,gauss(4,1))*temp(4))*gauss(3,2)&
@@ -1199,16 +1338,16 @@ subroutine wenorecon(i)
temp(l)=0.0
sum=0.0
do j=1,ko+1
- sum=rco9(j,l,1)/(weps+s(j))**2+sum
+ sum=rco9(j,l,1)/(weps+s(j,is))**2+sum
enddo
do j=1,ko+1
- ww(j)=rco9(j,l,1)/(weps+s(j))**2/sum
+ ww(j)=rco9(j,l,1)/(weps+s(j,is))**2/sum
enddo
do ll=1,ko+1
pp1(ll)=0.0
do j=1,ko+1
- pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll)
+ pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll,is)
enddo
enddo
do ll=1,ko+1
@@ -1235,16 +1374,16 @@ subroutine wenorecon(i)
temp(l)=0.0
sum=0.0
do j=1,ko+1
- sum=rco9(j,l,1)/(weps+s(j))**2+sum
+ sum=rco9(j,l,1)/(weps+s(j,is))**2+sum
enddo
do j=1,ko+1
- ww(j)=rco9(j,l,1)/(weps+s(j))**2/sum
+ ww(j)=rco9(j,l,1)/(weps+s(j,is))**2/sum
enddo
do ll=1,ko+1
pp1(ll)=0.0
do j=1,ko+1
- pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll)
+ pp1(ll)=pp1(ll)+coef9(j,ll,l)*u(0,i-ko-2+j+ll,is)
enddo
enddo
do ll=1,ko+1
@@ -1252,36 +1391,36 @@ subroutine wenorecon(i)
enddo
sum=0.0
do j=1,ko+1
- sum=rco9(j,l,2)/(weps+s(j))**2+sum
+ sum=rco9(j,l,2)/(weps+s(j,is))**2+sum
enddo
do j=1,ko+1
- ww(j)=rco9(j,l,2)/(weps+s(j))**2/sum
+ ww(j)=rco9(j,l,2)/(weps+s(j,is))**2/sum
enddo
temp(l)=-sigma(1,l)*temp(l)
do ll=1,ko+1
temp(l)=temp(l)+ww(ll)*pp1(ll)*sigma(2,l)
enddo
enddo
- ! reconstruction of u(1,i), u(2,i),u(3,i),u(4,i)
- u(1,i)=((fle(1,gauss(1,1))*temp(1)&
+ ! reconstruction of u(1,i,is), u(2,i,is),u(3,i,is),u(4,i,is)
+ u(1,i,is)=((fle(1,gauss(1,1))*temp(1)&
+fle(1,gauss(2,1))*temp(2))*gauss(1,2)&
+(fle(1,gauss(3,1))*temp(3)&
+fle(1,gauss(4,1))*temp(4))*gauss(3,2)&
+(fle(1,gauss(5,1))*temp(5)&
+fle(1,gauss(6,1))*temp(6))*gauss(5,2))*12.0
- u(2,i)=((fle(2,gauss(1,1))*temp(1)&
+ u(2,i,is)=((fle(2,gauss(1,1))*temp(1)&
+fle(2,gauss(2,1))*temp(2))*gauss(1,2)&
+(fle(2,gauss(3,1))*temp(3)&
+fle(2,gauss(4,1))*temp(4))*gauss(3,2)&
+(fle(2,gauss(5,1))*temp(5)&
+fle(2,gauss(6,1))*temp(6))*gauss(5,2))*180.0
- u(3,i)=((fle(3,gauss(1,1))*temp(1)&
+ u(3,i,is)=((fle(3,gauss(1,1))*temp(1)&
+fle(3,gauss(2,1))*temp(2))*gauss(1,2)&
+(fle(3,gauss(3,1))*temp(3)&
+fle(3,gauss(4,1))*temp(4))*gauss(3,2)&
+(fle(3,gauss(5,1))*temp(5)&
+fle(3,gauss(6,1))*temp(6))*gauss(5,2))*2800.0
- u(4,i)=((fle(4,gauss(1,1))*temp(1)&
+ u(4,i,is)=((fle(4,gauss(1,1))*temp(1)&
+fle(4,gauss(2,1))*temp(2))*gauss(1,2)&
+(fle(4,gauss(3,1))*temp(3)&
+fle(4,gauss(4,1))*temp(4))*gauss(3,2)&
@@ -1291,4 +1430,3 @@ subroutine wenorecon(i)
endif
return
end subroutine wenorecon
-
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/filename.f95 b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/filename.f95
new file mode 100644
index 0000000..41d8209
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/filename.f95
@@ -0,0 +1,41 @@
+! For Testing the DG&WENO subroutines
+! By Manuel Diaz
+
+program main
+ implicit none
+ ! Define
+ integer :: mo,mt,kflux,phase_quadrature_rule,init_value,ierror,n
+ integer :: k = 50
+ real :: cflc,tau,tprint,xmmm
+ ! Namelist
+ namelist /proj_list/ mt,kflux,cflc,tau,phase_quadrature_rule,init_value,ierror,tprint,n,xmmm,mo
+ ! character names
+ character(len=12) :: filename
+ character(len=12) :: format_string
+ character(len=12) :: name
+ integer :: counter
+
+ do counter=1, 10
+ if (counter < 10) then
+ format_string = "(A5,I1)"
+ else
+ format_string = "(A5,I2)"
+ endif
+
+ write (filename,format_string) "hello", counter
+ print *, trim(filename)
+ end do
+
+ write(name,"(A5,I2,A5)") "hello",10,"folks" ! <-- this is the key!!!
+ print *,'this is a filename: ', name
+
+ ! Prepare to read file
+ open(10, file="proj.in", form="FORMATTED", action="READ")
+ read(10, nml=proj_list)
+
+ ! Write namelist output
+ write(*, nml=proj_list)
+
+ print *, 'this is mt: ',mt
+
+end
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f95 b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f95
new file mode 100644
index 0000000..ec2c578
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f95
@@ -0,0 +1,50 @@
+! For Testing the DG&WENO Subroutines
+! By Manuel Diaz, NTU, 2013.06.16
+
+program main
+ ! Load modules
+ use IDmodule
+
+ ! No implicit definitions allowed
+ implicit none
+
+ ! Name of simulation
+ character(len=7) :: name = "SBBGK1d"
+
+ ! Define manually inputs
+ !integer :: theta = 0
+ !integer :: nx = 100
+ !integer :: P_deg = 3
+ !integer :: RK_stages = 3
+ !integer :: IC_case = 1
+ !integer :: fmodel = 1
+ !integer :: f_case = 1
+ !integer :: method = 1
+ !real :: tau = 1.0/10000
+
+ ! Expected Inputs
+ integer :: theta,nx,P_deg,RK_stages,IC_case,fmodel,f_case,method,kflux
+ real :: tau,tEnd,MM
+
+ ! IDs for files
+ character(len=100) :: IDn, IDf
+
+ ! Name list (comment if parameters where to be setup manually)
+ namelist /parameters_list/ name,theta,nx,P_deg,RK_stages,IC_case,fmodel,f_case,method,kflux,tau,tEnd,MM
+ ! this list need not to be in order
+
+ ! Read file with parameters
+ open(10, file="configuration.in", form="FORMATTED", action="READ")
+ read(10, nml=parameters_list)
+
+ ! Checking parametes read
+ write(*, nml=parameters_list)
+
+ ! Create IDs for simulation and result file
+ call ID_name(name,theta,nx,P_deg,RK_stages,tau,IC_case,fmodel,f_case,method,IDn,IDf)
+
+ ! Checking IDs
+ print *, 'IDn: ',IDn
+ print *, 'IDf: ',IDf
+
+end program
diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/proj.in b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/proj.in
new file mode 100644
index 0000000..aa57d47
--- /dev/null
+++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/proj.in
@@ -0,0 +1,13 @@
+$proj_list
+ mo=4,
+ mt=3,
+ kflux=1,
+ cflc=0.1,
+ tau=1e-4,
+ phase_quadrature_rule=0,
+ init_value=1,
+ ierror=0,
+ tprint=0.1,
+ n=100,
+ xmmm=1.0
+/ !End of list
diff --git a/myf90/ESBGK/com.f90 b/myf90/ESBGK/com.f90
deleted file mode 100644
index 5afdefe..0000000
--- a/myf90/ESBGK/com.f90
+++ /dev/null
@@ -1,33 +0,0 @@
-! the common block
-module data_module
- implicit none
- integer nx,md
- real nh
- parameter (nx=2560,md=5,nh=nx/2)
-
- real u(0:md,-md:nx+md),v(0:md,0:nx+1),hg(0:md,0:nx+1)
- real x(-md:nx+md),dx(-md:nx+md),bl(0:10),br(0:10),bi(0:md)
- real ai(0:10,0:10),am1(0:md),xx(200),uu(200),alim(3,3)
- real aintt(0:10,0:10,0:10)
- real coef(10,10,10),rco(10,2,10),gauss(10,2),coef9(5,5,6)
- real rco9(5,5,6)
- real sigma(2,6),gau(6,2),aii(0:10,0:10),sr
-
- real pi,cfld,dxmin,xm2
- integer mp,kcmax,indexmax,indexmin
- integer indexnum,indx1(0:nx), index
- integer mm
- real t,dt
-
- !! Input Param !!
- integer mo !ooa in space
- integer mt !ooa in time
- integer kflux !common flux '1=Roe, 2=LF, 3=LLF'
- real cflc !CFL number, 0.3 for P1, 0.18 for P2 and 0.1 for P3
- integer ierror !0=initial runs; 1=restart
- real tprint !input the terminal time
- integer n !n cells
- real xmmm !TVB constant M
-
- namelist/proj_list/mo,mt,kflux,cflc,ierror,tprint,n,xmmm
-end module data_module
diff --git a/myf90/ESBGK/inverse.f90 b/myf90/ESBGK/inverse.f90
deleted file mode 100644
index b7c88ff..0000000
--- a/myf90/ESBGK/inverse.f90
+++ /dev/null
@@ -1,112 +0,0 @@
-program main
-!====================================================================
-! Computing Inverse matrix
-! Method: Based on the Doolittle LU method
-!====================================================================
-implicit none
-integer, parameter :: n=3
-double precision a(n,n), c(n,n)
-integer i,j
-! matrix A
- data (a(1,i), i=1,3) / 3.0, 2.0, 4.0 /
- data (a(2,i), i=1,3) / 2.0, -3.0, 1.0 /
- data (a(3,i), i=1,3) / 1.0, 1.0, 2.0 /
-
-! print a header and the original matrix
- write (*,200)
- do i=1,n
- write (*,201) (a(i,j),j=1,n)
- end do
-
- call inverse(a,c,n)
-
-! print the inverse matrix C = A^{-1}
- write (*,202)
- do i = 1,n
- write (*,201) (c(i,j),j=1,n)
- end do
-200 format (' Computing Inverse matrix ',/,/, &
- ' Matrix A')
-201 format (6f12.6)
-202 format (/,' Inverse matrix A^{-1}')
-end
-
- subroutine inverse(a,c,n)
-!============================================================
-! Inverse matrix
-! Method: Based on Doolittle LU factorization for Ax=b
-! Alex G. December 2009
-!-----------------------------------------------------------
-! input ...
-! a(n,n) - array of coefficients for matrix A
-! n - dimension
-! output ...
-! c(n,n) - inverse matrix of A
-! comments ...
-! the original matrix a(n,n) will be destroyed
-! during the calculation
-!===========================================================
-implicit none
-integer n
-double precision a(n,n), c(n,n)
-double precision L(n,n), U(n,n), b(n), d(n), x(n)
-double precision coeff
-integer i, j, k
-
-! step 0: initialization for matrices L and U and b
-! Fortran 90/95 aloows such operations on matrices
-L=0.0
-U=0.0
-b=0.0
-
-! step 1: forward elimination
-do k=1, n-1
- do i=k+1,n
- coeff=a(i,k)/a(k,k)
- L(i,k) = coeff
- do j=k+1,n
- a(i,j) = a(i,j)-coeff*a(k,j)
- end do
- end do
-end do
-
-! Step 2: prepare L and U matrices
-! L matrix is a matrix of the elimination coefficient
-! + the diagonal elements are 1.0
-do i=1,n
- L(i,i) = 1.0
-end do
-! U matrix is the upper triangular part of A
-do j=1,n
- do i=1,j
- U(i,j) = a(i,j)
- end do
-end do
-
-! Step 3: compute columns of the inverse matrix C
-do k=1,n
- b(k)=1.0
- d(1) = b(1)
-! Step 3a: Solve Ld=b using the forward substitution
- do i=2,n
- d(i)=b(i)
- do j=1,i-1
- d(i) = d(i) - L(i,j)*d(j)
- end do
- end do
-! Step 3b: Solve Ux=d using the back substitution
- x(n)=d(n)/U(n,n)
- do i = n-1,1,-1
- x(i) = d(i)
- do j=n,i+1,-1
- x(i)=x(i)-U(i,j)*x(j)
- end do
- x(i) = x(i)/u(i,i)
- end do
-! Step 3c: fill the solutions x(n) into column k of C
- do i=1,n
- c(i,k) = x(i)
- end do
- b(k)=0.0
-end do
-end subroutine inverse
diff --git a/myf90/ESBGK/main.f90 b/myf90/ESBGK/main.f90
deleted file mode 100644
index 51d33fb..0000000
--- a/myf90/ESBGK/main.f90
+++ /dev/null
@@ -1,85 +0,0 @@
-program matrix_operations
-! Example code on how to deal with array operations in Fortran 95. I'm also creating some functions so that I can imitate Matlab stile of programming. ;D
-! coded by Manuel Diaz, NTU, 2013.07.11
-! manuel.ade'at'gmail.com
-
-!PACKAGES
-use dispmodule
-use mathmodule!, only: ones, zeros, inverse, transposev
-
-!Preamble
-implicit none
-
-integer, parameter :: n=3, m=3
-real, dimension(n,m) :: A,B,C !Matrix Arrays
-real, dimension(1,m) :: x,z,ty !Row vectors
-real, dimension(n,1) :: y,w,tz !Column vectors
-real, dimension(:,:), allocatable :: xx,zz !unknown dimension
-
-! allocate variables in memory
-allocate (zz(n,1))
-allocate (xx(n,1))
-
-! Build arrays
-A = transpose( reshape( (/1,2,3,4,5,6,7,8,9/),shape(A)))
-B = transpose( reshape( (/2,0,0,0,2,0,0,0,2/),shape(B)))
-call disp('A = ',A)
-call disp('B = ',B)
-
-! Build vector arrays
-x = reshape( (/1,2,3/),shape(x))
-y = reshape( (/2,2,2/),shape(y))
-w = reshape(x,(/size(x,2),size(x,1)/) )
-zz= reshape(x,(/size(x,1),size(x,2)/) )
-xx= transpose(x)
-call disp('x = ',x)
-call disp('y = ',y)
-call disp('w = ',w)
-call disp('zz= ',zz)
-call disp('xx= ',xx)
-
-! Vector Addition
-call disp(' A + B = ',A+B)
-call disp('x+T(y) = ',x+transpose(y))
-
-! Dot product
-call disp(' A * B = ',A*B)
-call disp('A*B*T(A)=',A*B*transpose(A))
-call disp('x*T(y) = ',x*transpose(y))
-call disp(' A*A*B = ',A*A*B)
-C = sum(B*A*transpose(B))
-call disp(' C = ',C)
-
-! matrix multiplication and inner product
-call disp('matmul(tA,B)= ',matmul(transpose(A),B))
-call disp('matmul(A,B) = ',matmul(A,B))
-call disp('matmul(x,y) = ',matmul(x,y))
-call disp('matmul(A,x) = ',matmul(A,y))
-call disp('matmul(ty,B)= ',matmul(transpose(y),B))
-
-! Sum of elements in arrays
-call disp('sum(A,1) = ',sum(A,dim=1))
-z = reshape(sum(A,dim=1),shape(z))
-call disp('sum(A,1) = ',z)
-!call disp('sum(A,1) = ',reshape(sum(A,dim=1),shape(z)) <-can't
-y = reshape(sum(A,dim=1),shape(y))
-call transposev(y,ty)
-call disp('tsum(A,1)= ',ty)
-call disp('sum(A,2) = ',sum(A,dim=2))
-
-! Create zeros array
-call disp('ones(3,3) = ',ones(3,3))
-call disp('ones(1,3) = ',ones(1,3))
-
-! Create ones array
-call disp('zeros(4,4)= ',zeros(4,4))
-call disp('zeros(4,1)= ',zeros(4,1))
-
-! Invert array
-call disp('inv(B^2) = ',inverse(matmul(B,B)))
-
-! Print array using dispmodule.f90
-!call disp('C = ',C)
-!call disp('z = ',z)
-
-end program matrix_operations
diff --git a/myf90/ESBGK/proj.in b/myf90/ESBGK/proj.in
deleted file mode 100644
index 0485469..0000000
--- a/myf90/ESBGK/proj.in
+++ /dev/null
@@ -1,10 +0,0 @@
-$proj_list
-mo=4,
-mt=3,
-kflux=1,
-cflc=0.1,
-ierror=0,
-tprint=1.,
-n=100,
-xmmm=1.
-/