A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
test_tl_ad.f90
Go to the documentation of this file.
1 
2 ! test_tl_ad.f90
3 !
4 !> Tests for the Tangent Linear (TL) and Adjoint (AD) model versions
5 !> of MAOOAM.
6 !
7 !> @copyright
8 !> 2016 Lesley De Cruz & Jonathan Demaeyer.
9 !> See LICENSE.txt for license information.
10 !
11 !---------------------------------------------------------------------------!
12 
13 
14 PROGRAM test_tl_ad
15  USE params, only:ndim,dt,t_trans
16  USE aotensor_def, only: init_aotensor
20  USE stoch_mod, only: gasdev
21  IMPLICIT NONE
22 
23  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: y0_IC,y0,y0prime,dy0,dy0_bis
24  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: y1,y1prime
25  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dy,dy_bis
26  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dy1,dy1_tl,dy1_bis_tl,dy1_ad,dy1_bis_ad
27  REAL(KIND=8) :: t=0.d0
28  REAL(KIND=8) :: norm1,norm2
29  INTEGER :: i,n
30 
31  ! Compute the tensors
32 
33  CALL init_aotensor
34  CALL init_tltensor
35  CALL init_adtensor
36 
37  CALL init_integrator ! Initialize the model integrator
38  CALL init_tl_ad_integrator ! Initialize the TL & AD integrator
39 
40  ALLOCATE(y0_ic(0:ndim),dy(0:ndim),y0(0:ndim),y0prime(0:ndim)&
41  &,y1(0:ndim))!, STAT=AllocStat)
42  ALLOCATE(y1prime(0:ndim),dy0(0:ndim),dy0_bis(0:ndim))!, STAT&
43  ! &=AllocStat)
44  ALLOCATE(dy1(0:ndim),dy1_tl(0:ndim),dy_bis(0:ndim)&
45  &,dy1_bis_tl(0:ndim),dy1_ad(0:ndim),dy1_bis_ad(0:ndim))!, STAT&
46  ! &=AllocStat)
47 
48 
49  ! Test Taylor property for the Tangent Linear.
50  ! lim(\lambda->0) M(x+\lambda dx) - M(x) / M'(\lambda dx) = 1
51 
52  y0_ic(0)=1.d0
53 
54  ! Set all values to random.
55 
56  DO i=1,ndim
57  y0_ic(i)=0.01*gasdev()
58  ENDDO
59 
60  ! Evolve during transient period
61 
62  ! PRINT*, 'Random values:',y0_IC(0:ndim)
63  DO WHILE (t<t_trans)
64  CALL step(y0_ic,t,dt,y0)
65  y0_ic=y0
66  END DO
67  print*, 'Initial values:',y0_ic(0:ndim)
68 
69  ! Test 1: Taylor test
70  ! Integrate the original model by one step, integrate with perturbed
71  ! IC, and test if the difference approximates the TL
72 
73  DO n=0,6
74  ! Small perturbation.
75  dy=2.d0**(-n)/sqrt(real(ndim))
76  dy(0)=0.d0
77  print*, "Perturbation size:",dot_product(dy,dy)
78 
79  y0 = y0_ic*1
80  y0prime = y0 + dy
81  CALL step(y0,t,dt,y1)
82  CALL step(y0prime,t,dt,y1prime)
83 
84  dy1 = y1prime - y1
85 
86  dy0 = dy*1
87  CALL tl_step(dy0,y0_ic,t,dt,dy1_tl)
88 
89  ! Don't forget to set 0'th component to 0...
90  dy1(0)=0.d0
91  dy1_tl(0)=0.d0
92 
93  print*, "Resulting difference in trajectory: (epsilon ~ 2^-",n,")"
94  print*, "diff: ",dot_product(dy1,dy1)
95  print*, "tl: ",dot_product(dy1_tl,dy1_tl)
96  print*, "ratio: ",dot_product(dy1,dy1)/dot_product(dy1_tl,dy1_tl)
97  END DO
98 
99  ! Test 2: Adjoint Identity: <M(TL).x,y> = <x,M(AD).y>
100 
101  DO i=1,100
102  ! Any perturbation.
103  DO n=1,ndim
104  dy(n)=gasdev()
105  END DO
106  dy(0)=0.d0
107 
108  DO n=1,ndim
109  dy_bis(n)=gasdev()
110  END DO
111  dy_bis(0)=0.d0
112 
113  ! Calculate M(TL).x in dy1_tl
114  dy0 = dy*1
115  CALL tl_step(dy0,y0_ic,t,dt,dy1_tl)
116 
117  ! Calculate M(AD).x in dy1_ad
118  CALL ad_step(dy0,y0_ic,t,dt,dy1_ad)
119 
120  ! Calculate M(TL).y in dy1_bis_tl
121  CALL tl_step(dy0_bis,y0_ic,t,dt,dy1_bis_tl)
122 
123  ! Calculate M(AD).y in dy1_bis_ad
124  CALL ad_step(dy0_bis,y0_ic,t,dt,dy1_bis_ad)
125 
126  ! Calculate norm <M(TL).x,y>
127  norm1 = dot_product(dy1_tl,dy0_bis)
128  ! Calculate norm <x,M(AD).y>
129  norm2 = dot_product(dy0,dy1_bis_ad)
130 
131  print*, "<M(TL).x,y> = ", norm1
132  print*, "<x,M(AD).y> = ", norm2
133  print*, "Ratio = ", norm1/norm2
134 
135  ! Calculate norm <M(TL).y,x>
136  norm1 = dot_product(dy1_bis_tl,dy0)
137  ! Calculate norm <y,M(AD).x>
138  norm2 = dot_product(dy0_bis,dy1_ad)
139 
140  print*, "<M(TL).y,x> = ", norm1
141  print*, "<y,M(AD).x> = ", norm2
142  print*, "Ratio = ", norm1/norm2
143  END DO
144 
145 END PROGRAM test_tl_ad
146 
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
subroutine, public step(y, t, dt, res)
Routine to perform an integration step (Heun algorithm). The incremented time is returned.
real(kind=8) function, public gasdev()
Definition: stoch_mod.f90:32
Tangent Linear (TL) and Adjoint (AD) model versions of MAOOAM. Integrators module.
subroutine, public ad_step(y, ystar, t, dt, res)
Routine to perform an integration step (Heun algorithm) of the adjoint model. The incremented time is...
subroutine, public init_tltensor
Routine to initialize the TL tensor.
program test_tl_ad
Tests for the Tangent Linear (TL) and Adjoint (AD) model versions of MAOOAM.
Definition: test_tl_ad.f90:14
The equation tensor for the coupled ocean-atmosphere model with temperature which allows for an exten...
subroutine, public init_adtensor
Routine to initialize the AD tensor.
Module with the integration routines.
The model parameters module.
Definition: params.f90:18
Utility module containing the stochastic related routines.
Definition: stoch_mod.f90:15
Tangent Linear (TL) and Adjoint (AD) model versions of MAOOAM. Tensors definition module...
subroutine, public init_tl_ad_integrator
Routine to initialise the integration buffers.
real(kind=8) dt
Integration time step.
Definition: params.f90:77
subroutine, public init_integrator
Routine to initialise the integration buffers.
real(kind=8) t_trans
Transient time period.
Definition: params.f90:75
subroutine, public init_aotensor
Subroutine to initialise the aotensor tensor.
subroutine, public tl_step(y, ystar, t, dt, res)
Routine to perform an integration step (Heun algorithm) of the tangent linear model. The incremented time is returned.