Merge branch 'master' of https://github.com/qiqi/fds

.travis.yml

@@ -5,10 +5,7 @@ dist: trusty
 
 python:
     - 2.7
-    - 3.3
-    - 3.4
     - 3.5
-    - nightly
 
 # command to install dependencies
 before_install:

tests/data/pisofoam_restart/0/U

@@ -731,7 +731,9 @@ boundaryField
     }
     outflow
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
     }
 }
 

tests/data/pisofoam_restart/0/U_0

@@ -731,7 +731,9 @@ boundaryField
     }
     outflow
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
     }
 }
 

tests/solvers/Lorenz96/Lorenz96.f90

@@ -0,0 +1,119 @@
+! Lorenz ' 96 system
+
+module Lorenz96
+	! system parameters:
+    real(kind=8), parameter :: M = 8.1d0
+contains
+
+subroutine Xnp1(X,D,Xnp1_res)
+
+
+    implicit none
+	integer, intent(in) :: D
+	real(kind=8), intent(in), dimension(D) :: X
+	real(kind=8), intent(out), dimension(D-3) :: Xnp1_res
+	integer :: i
+	real(kind=8) :: dt
+
+	dt = 0.01d0
+	do i=3,D-1
+		Xnp1_res(i-2) = (-X(i-2) + X(i+1))*X(i-1) - X(i) + M	
+	end do
+	!dXdt(1) = (-X(D-1) + X(2))*X(D) - X(1) + M
+  	!dXdt(2) = (-X(D) + X(3))*X(1) - X(2) + M
+	!dXdt(D) = (-X(D-2) + X(1))*X(D-1) - X(D) + M
+	Xnp1_res = X(3:D-1) + dt*Xnp1_res
+end subroutine Xnp1
+
+subroutine dfdX(X,D,dfdX_res)
+
+	implicit none
+	integer, intent(in) :: D
+	real(kind=8), intent(in), dimension(D) :: X
+	real(kind=8), intent(out), dimension(D-3,D-3) :: dfdX_res 	
+	integer :: i,j
+
+	dfdX_res = 0.d0
+	do j=3,D-4
+		i = j + 2	
+		dfdX_res(j,j) = -1.d0
+		dfdX_res(j,j-1) = -X(i-2) + X(i+1)
+		dfdX_res(j,j+1) = X(i-1)
+		dfdX_res(j,j-2) = -X(i-1)  					
+
+	enddo	
+
+	dfdX_res(1,1) = -1.d0
+	dfdX_res(1,2) = X(2)
+	dfdX_res(2,2) = -1.d0	 	
+	dfdX_res(2,1) = X(5) - X(2)
+	dfdX_res(2,3) = X(3)	
+
+	dfdX_res(D-3,D-3) = -1.d0
+	dfdX_res(D-3,D-4) = X(D) - X(D-3)
+	dfdX_res(D-3,D-5) = -X(D-2) 
+
+end subroutine dfdX  
+
+subroutine dvdt(X,D,v1,dvdt_res)
+
+		implicit none
+		integer, intent(in) :: D
+		real(kind=8), intent(in), dimension(D) :: X
+		real(kind=8), intent(out), dimension(D-3,1) :: dvdt_res
+		real(kind=8), dimension(D-3,D-3) :: dfdX_res
+		real(kind=8), intent(in), dimension(D,1) :: v1
+		real(kind=8), dimension(D-3,D) :: dfdX_ext
+		integer :: i		
+
+		call dfdX(X,D,dfdX_res)
+		dfdX_ext = 0.d0
+		dfdX_ext(1:D-3,3:D-1) = dfdX_res
+
+		dfdX_ext(1,2) = X(4) - X(1)
+		dfdX_ext(1,1) = -X(2)  				
+
+		dfdX_ext(2,1) = -X(3) 	
+		dfdX_ext(D-3,D) = X(D-2) 
+
+		dvdt_res = matmul(dfdX_ext,v1)	
+end subroutine dvdt
+subroutine rk45_full(X,D,v,vnp1)
+!Assumes full perturbation vector.
+	implicit none
+	integer, intent(in) :: D
+	real(kind=8) , intent(in), dimension(D) :: X,v
+	real(kind=8) , intent(out), dimension(D-3,1) :: vnp1
+	real(kind=8) , dimension(D,1) :: v1,k1,k2,k3,k4
+	real(kind=8) , dimension(D-3,1):: dvdt_res
+	real(kind=8) :: dt
+
+
+
+	dt = 0.001d0	
+	v1 = reshape(v,[D,1])
+	call dvdt(X,D,v1,dvdt_res)	
+	k1(3:D-1,1) = dt*dvdt_res(:,1)
+	k1(1,1) = k1(D-2,1)
+	k1(2,1) = k1(D-1,1)
+	k1(D,1) = k1(3,1)
+	call dvdt(X,D,v1 + 0.5d0*k1,dvdt_res)
+	k2(3:D-1,1) = dt*dvdt_res(:,1)
+	k2(1,1) = k2(D-2,1)
+	k2(2,1) = k2(D-1,1)
+	k2(D,1) = k2(3,1)
+	call dvdt(X,D,v1 + 0.5d0*k2,dvdt_res)
+	k3(3:D-1,1) = dt*dvdt_res(:,1)
+	k3(1,1) = k3(D-2,1)
+	k3(2,1) = k3(D-1,1)
+	k3(D,1) = k3(3,1)
+	call dvdt(X,D,v1 + k3,dvdt_res)
+	k4(3:D-1,1) = dt*dvdt_res(:,1)
+	vnp1(:,1) = v1(3:D-1,1) + 1.d0/6.d0*k1(3:D-1,1) + &
+		 1.d0/3.d0*k2(3:D-1,1) + 1.d0/3.d0*k3(3:D-1,1) + &
+		1.d0/6.d0*k4(3:D-1,1)
+
+
+end subroutine rk45_full 
+
+end module Lorenz96

tests/solvers/Lorenz96/Makefile

@@ -0,0 +1,32 @@
+
+OBJECTS1 = Lorenz96.o ensemble_tangent.o 
+OBJECTS2 = Lorenz96.o test1.o 
+MODULES = Lorenz96.mod
+LFLAGS = -lblas -llapack
+FFLAGS = -g
+NUM_PROCS ?= 4   
+
+FC = mpif90
+
+.PHONY: test1 enstan clean 
+enstan: enstan.exe
+	mpiexec -n $(NUM_PROCS) ./enstan.exe
+
+enstan.exe: $(MODULES) $(OBJECTS1)
+	$(FC) $(FFLAGS) $(LFLAGS) $(OBJECTS1) -o enstan.exe
+
+test1: test1.exe
+	mpiexec -n $(NUM_PROCS) ./test1.exe
+
+test1.exe: $(MODULES) $(OBJECTS2)
+	$(FC) $(FFLAGS) $(LFLAGS) $(OBJECTS2) -o test1.exe
+
+%.o : %.f90
+	$(FC) $(FFLAGS) $(LFLAGS) -c  $< 
+
+%.mod: %.f90
+	$(FC) $(FFLAGS) $(LFLAGS) -c $<
+
+clean:
+	rm -f *.o *.exe *.mod
+

tests/solvers/Lorenz96/ensemble_tangent.f90

@@ -0,0 +1,144 @@
+! Lorenz'96
+
+program ensemble_tangent
+
+    use Lorenz96
+    use mpi 
+
+    implicit none
+	real(kind=8), dimension(:), allocatable ::X,v,Xp,Xpnp1_res,Xnp1_res
+	real(kind=8), dimension(:), allocatable :: dXdt,g
+	real(kind=8), dimension(:,:), allocatable :: dfdX_res, vnp1_res
+	integer :: i, me, ierr, nprocs, Dproc, D, ns, ns_proc, j, Dext
+    integer :: istart, iend, lproc, rproc	
+	integer, allocatable :: seed(:)
+	integer :: rsize, req1, req2
+	integer, dimension(MPI_STATUS_SIZE) :: mpistatus
+ 	real(kind=8), pointer, dimension(:) :: p
+	real(kind=8) :: dt, dXavgds, dXavgds_avg_proc, L1, thetaEA_mean, thetaEA_var	
+	integer :: thefile, T, tau,k,k1,s
+	real(kind=8), dimension(:), allocatable :: dXavgds_all,thetaEA
+	integer, dimension(MPI_STATUS_SIZE) :: status
+    call mpi_init(ierr)
+    call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr)
+    call mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
+
+	tau = 1000
+	D = 40	
+	dt = 0.01d0
+	T = 100000
+	L1 = 0.d0
+	ns = 10000
+	ns_proc = ns/nprocs
+	Dext = D+3
+	s = T/tau
+	allocate(X(1:Dext),v(1:Dext),vnp1_res(1:D,1),Xnp1_res(1:D), &
+	Xpnp1_res(1:D), Xp(1:Dext), g(1:D), dXavgds_all(1:ns),thetaEA(1:ns/s))
+	if(me==0) then
+		k=0
+	endif
+
+	istart = 3
+	iend = Dext - 1
+	g = 0.d0
+	g(D) = 1.d0
+
+	dXavgds_avg_proc = 0.d0	
+	do j = 1, ns_proc
+
+
+			call RANDOM_SEED(SIZE=rsize)
+			allocate(seed(rsize))
+			call RANDOM_SEED(PUT=seed)	
+			call RANDOM_NUMBER(X)
+			v = 0.d0
+			Xp = X
+
+			v(istart) = 1.d0
+			Xp(istart) = Xp(istart) + 0.001d0
+
+
+		!	call MPI_FILE_OPEN(MPI_COMM_WORLD, 'test', &
+		!			MPI_MODE_WRONLY + MPI_MODE_CREATE, &
+		!			MPI_INFO_NULL, thefile, ierr)
+		!			open(unit=30+me, file='initialdata.dat')
+		!			write(30+me,*) X(istart:iend)	
+		!			close(30+me)	
+
+		    dXavgds = 0.d0
+
+			do i = 1, tau	
+
+
+
+				v(1) = v(Dext-2)
+                v(2) = v(Dext-1)
+                v(Dext) = v(istart)
+
+
+
+                X(1) = X(Dext-2)
+                X(2) = X(Dext-1)
+                X(Dext) = X(istart)
+
+                Xp(1) = Xp(Dext-2)
+                Xp(2) = Xp(Dext-1)
+                Xp(Dext) = Xp(istart)
+
+				call Xnp1(X,Dext,Xnp1_res)
+				call Xnp1(Xp,Dext,Xpnp1_res)
+				call rk45_full(X,Dext,v,vnp1_res)
+
+				if(me == 0) then
+				!Compute lift and drag.	
+				end if
+				X(istart:iend) = Xnp1_res
+				Xp(istart:iend) = Xpnp1_res
+				v(istart:iend) = vnp1_res(:,1)
+
+				if(me==0 .and. i==1000) then	
+					open(unit=20, file='X.dat')
+					write(20,*) X(istart:iend)	
+					open(unit=21, file='Xp.dat')
+					write(21,*) Xp(istart:iend)		
+				end if
+
+
+				dXavgds = dXavgds + DOT_PRODUCT(v(istart:iend),g)
+			enddo
+			close(20)
+			close(21)
+			dXavgds = dXavgds/tau
+
+			call MPI_SEND(dXavgds,1,MPI_DOUBLE_PRECISION, &
+							0, me, MPI_COMM_WORLD, ierr)
+
+			if(me==0) then
+				do k1 = 1,nprocs
+					k = k + 1
+					call MPI_RECV(dXavgds_all(k), &
+						1, MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE, &
+						MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
+				end do
+			end if      
+			deallocate(seed)
+	enddo
+
+
+
+	if(me==0) then
+		do k1 = 1,ns/s
+			thetaEA(k1) = sum(dXavgds_all((k1-1)*s+1:k1*s))/s
+		end do
+		thetaEA_mean = sum(thetaEA)/(ns/s)
+		thetaEA_var = 0.d0
+		do k1 = 1,ns/s
+			thetaEA_var = thetaEA_var + (thetaEA(k1)-thetaEA_mean)**2.0
+		end do
+		thetaEA_var = thetaEA_var/(ns/s)
+	end if
+
+
+	call mpi_finalize(ierr)	
+
+end program ensemble_tangent