Skip to content

Commit

Permalink
Merge 43a9dc3 into 67fec74
Browse files Browse the repository at this point in the history
  • Loading branch information
Nisha Chandramoorthy committed Oct 18, 2016
2 parents 67fec74 + 43a9dc3 commit 5c3e786
Show file tree
Hide file tree
Showing 5 changed files with 393 additions and 1 deletion.
111 changes: 111 additions & 0 deletions tests/solvers/Lorenz96/Lorenz96.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
! 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(X,D,v,vnp1)

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-3,1) :: v1,dvdt_res,k1,k2,k3,k4
real(kind=8) :: dt



dt = 0.001d0
v1 = reshape(v,[D-3,1])
call dvdt(X,D,v1,dvdt_res)
k1 = dt*dvdt_res
call dvdt(X,D,v1 + 0.5d0*k1,dvdt_res)
k2 = dt*dvdt_res
call dvdt(X,D,v1 + 0.5d0*k2,dvdt_res)
k3 = dt*dvdt_res
call dvdt(X,D,v1 + k3,dvdt_res)
k4 = dt*dvdt_res

vnp1 = v1 + 1.d0/6.d0*k1 + &
1.d0/3.d0*k2 + 1.d0/3.d0*k3 + &
1.d0/6.d0*k4


end subroutine rk45

end module Lorenz96
32 changes: 32 additions & 0 deletions tests/solvers/Lorenz96/Makefile
View file
Original file line number Diff line number Diff line change
@@ -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

107 changes: 107 additions & 0 deletions tests/solvers/Lorenz96/ensemble_tangent.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
! 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, ncyc, 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, Xavg
integer :: thefile, ns, T
call mpi_init(ierr)
call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr)
call mpi_comm_rank(MPI_COMM_WORLD, me, ierr)

ncyc = 2000
D = 40
dt = 0.01d0
T = 2000

ns = T/ncyc
ns_proc = ns/nprocs
Dext = D+3
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))

istart = 3
iend = D
g = 0.d0
g(D) = 1.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

if(me == 0) then
v(istart) = 1.d0
Xp(istart) = Xp(istart) + 0.001d0
end if

! 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)

Xavg = 0.d0
do i = 1, ncyc



X(1) = X(Dext-1)
X(2) = X(Dext)
X(Dext) = X(istart)

Xp(1) = Xp(Dext-1)
Xp(2) = Xp(Dext)
Xp(Dext) = Xp(istart)

call Xnp1(X,Dext,Xnp1_res)
call Xnp1(Xp,Dext,Xpnp1_res)
call rk45(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


Xavg = Xavg + DOT_PRODUCT(v(istart:iend),g)
enddo
close(20)
close(21)
Xavg = Xavg/ncyc
if(me==0) then

print *, 1.d0/dt/ncyc*log(abs(Xp(istart)-X(istart))/0.001d0)
endif
enddo

call mpi_finalize(ierr)

end program ensemble_tangent
139 changes: 139 additions & 0 deletions tests/solvers/Lorenz96/test1.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
! Lorenz'96

program test1

use Lorenz96
use mpi

implicit none
real(kind=8), dimension(:), allocatable ::X,v,Xp,Xpnp1_res,Xnp1_res
real(kind=8), dimension(:), allocatable :: dXdt
real(kind=8), dimension(:,:), allocatable :: dfdX_res, vnp1_res
integer :: i, me, ierr, nprocs, Dproc, D
integer :: istart, iend, ncyc, 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
integer :: thefile
call mpi_init(ierr)
call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr)
call mpi_comm_rank(MPI_COMM_WORLD, me, ierr)

ncyc = 2000
D = 40
Dproc = D/nprocs
dt = 0.01d0
istart = me*Dproc + 1
iend = min((me + 1)*Dproc, D)
lproc = MOD(me + nprocs - 1, nprocs)
rproc = MOD(me + 1, nprocs)

allocate(X(istart-2:iend+1),v(istart-2:iend+1),vnp1_res(istart:iend,1),Xnp1_res(istart:iend), &
Xpnp1_res(istart:iend), Xp(istart-2:iend+1))

call RANDOM_SEED(SIZE=rsize)
allocate(seed(rsize))
call RANDOM_SEED(PUT=seed)
call RANDOM_NUMBER(X)
v = 0.d0
if(me == 0) then
v(istart) = 1.d0
Xp(istart) = Xp(istart) + 0.001d0
end if

! 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)


do i = 1, ncyc

call mpi_isend(X(istart), 1, &
MPI_DOUBLE_PRECISION, &
lproc, &
1, MPI_COMM_WORLD, req1, ierr)

call mpi_isend(X(iend-1:iend), &
2, MPI_DOUBLE_PRECISION, &
rproc, &
2, MPI_COMM_WORLD, req2, ierr)


call mpi_recv(X(istart-2:istart-1), &
2, MPI_DOUBLE_PRECISION, lproc, &
2, MPI_COMM_WORLD, mpistatus, ierr)

call mpi_recv(X(iend+1), &
1, MPI_DOUBLE_PRECISION, rproc, &
1, MPI_COMM_WORLD, mpistatus, ierr)

call mpi_isend(v(istart), 1, &
MPI_DOUBLE_PRECISION, &
lproc, &
3, MPI_COMM_WORLD, req1, ierr)

call mpi_isend(v(iend-1:iend), &
2, MPI_DOUBLE_PRECISION, &
rproc, &
4, MPI_COMM_WORLD, req2, ierr)


call mpi_recv(v(istart-2:istart-1), &
2, MPI_DOUBLE_PRECISION, lproc, &
4, MPI_COMM_WORLD, mpistatus, ierr)

call mpi_recv(v(iend+1), &
1, MPI_DOUBLE_PRECISION, rproc, &
3, MPI_COMM_WORLD, mpistatus, ierr)

call mpi_isend(Xp(istart), 1, &
MPI_DOUBLE_PRECISION, &
lproc, &
5, MPI_COMM_WORLD, req1, ierr)

call mpi_isend(Xp(iend-1:iend), &
2, MPI_DOUBLE_PRECISION, &
rproc, &
6, MPI_COMM_WORLD, req2, ierr)


call mpi_recv(Xp(istart-2:istart-1), &
2, MPI_DOUBLE_PRECISION, lproc, &
6, MPI_COMM_WORLD, mpistatus, ierr)

call mpi_recv(Xp(iend+1), &
1, MPI_DOUBLE_PRECISION, rproc, &
5, MPI_COMM_WORLD, mpistatus, ierr)

call Xnp1(X,Dproc+3,Xnp1_res)
call Xnp1(Xp,Dproc+3,Xpnp1_res)
call rk45(X,Dproc+3,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
enddo
close(20)
close(21)
if(me==0) then
print *,"direct sensitivity: ", 1.d0/dt/ncyc*log(abs(Xp(istart)-X(istart))/0.001d0)

endif
call mpi_finalize(ierr)

end program test1
Loading

0 comments on commit 5c3e786

Please sign in to comment.