A Modular Arbitrary-Order Ocean-Atmosphere Model -- Stochastic implementation
WL_tensor.f90
Go to the documentation of this file.
1 
2 ! WL_tensor.f90
3 !
4 !> The WL tensors used to integrate the model
5 !
6 !> @copyright
7 !> 2018 Jonathan Demaeyer.
8 !> See LICENSE.txt for license information.
9 !
10 !---------------------------------------------------------------------------!
11 !
12 !> @remark
13 !>
14 !>
15 !
16 !---------------------------------------------------------------------------!
17 
18 MODULE wl_tensor
19 
20  !-----------------------------------------------------!
21  ! !
22  ! Preamble and variables declaration !
23  ! !
24  !-----------------------------------------------------!
25 
26  USE tensor
27  USE dec_tensor
28  USE params, only:ndim
29  USE stoch_params, only:mems
30  USE util, only: mat_trace, mat_contract
31  USE corr_tensor
32  USE corrmod, only: corr_i_full,mean_full
33 
34 
35  IMPLICIT NONE
36 
37  PRIVATE
38 
39  PUBLIC :: init_wl_tensor
40 
41  ! M1 term vectors
42  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: m11 !< First component of the M1 term
43  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: m12 !< Second component of the M1 term
44  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: m13 !< Third component of the M1 term
45  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, PUBLIC :: m1tot !< Total \f$M_1\f$ vector
46 
47  ! M2 term tensors
48  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: m21 !< First tensor of the M2 term
49  TYPE(coolist), DIMENSION(:), ALLOCATABLE, PUBLIC :: m22 !< Second tensor of the M2 term
50 
51  ! M3 terms tensors
52  ! Tensors for the linear terms (no L3 for WL)
53  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: l1 !< First linear tensor
54  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: l2 !< Second linear tensor
55  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: l4 !< Fourth linear tensor
56  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: l5 !< Fifth linear tensor
57  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: ltot !< Total linear tensor
58 
59  ! Tensors for the quadratic terms
60  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: b1 !< First quadratic tensor
61  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: b2 !< Second quadratic tensor
62  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: b3 !< Third quadratic tensor
63  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: b4 !< Fourth quadratic tensor
64  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: b14 !< Joint 1st and 4th tensors
65  TYPE(coolist), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: b23 !< Joint 2nd and 3rd tensors
66 
67  TYPE(coolist4), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: mtot !< Tensor for the cubic terms
68 
69  REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dumb_vec !< Dummy vector
70  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat1 !< Dummy matrix
71  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat2 !< Dummy matrix
72  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat3 !< Dummy matrix
73  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dumb_mat4 !< Dummy matrix
74 
75  LOGICAL, PUBLIC :: m12def,m21def,m22def,ldef,b14def,b23def,mdef !< Boolean to (de)activate the computation of the terms
76 
77  !-----------------------------------------------------!
78  ! !
79  ! End of preamble !
80  ! !
81  !-----------------------------------------------------!
82 
83 CONTAINS
84 
85  !-----------------------------------------------------!
86  ! !
87  ! Initialisation routine !
88  ! !
89  !-----------------------------------------------------!
90 
91 
92  !> Subroutine to initialise the WL tensor
93  SUBROUTINE init_wl_tensor
94  INTEGER :: AllocStat,i,j,k,m
95 
96  print*, 'Initializing the decompostion tensors...'
97  CALL init_dec_tensor
98  print*, "Initializing the correlation matrices and tensors..."
99  CALL init_corr_tensor
100 
101  !M1 part
102  print*, "Computing the M1 terms..."
103 
104  ALLOCATE(m11(0:ndim), m12(ndim), m13(0:ndim), m1tot(0:ndim), stat=allocstat)
105  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
106 
107  ALLOCATE(dumb_mat1(ndim,ndim), dumb_mat2(ndim,ndim), stat=allocstat)
108  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
109 
110  ALLOCATE(dumb_mat3(ndim,ndim), dumb_mat4(ndim,ndim), stat=allocstat)
111  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
112 
113  ALLOCATE(dumb_vec(ndim), stat=allocstat)
114  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
115 
116  !M11
117  m11=0.d0
118  ! CALL coo_to_mat_ik(Lxy,dumb_mat1)
119  ! M11(1:ndim)=matmul(dumb_mat1,mean_full(1:ndim))
120 
121  !M12
122  ! dumb_mat2=0.D0
123  ! DO i=1,ndim
124  ! CALL coo_to_mat_i(i,Bxxy,dumb_mat1)
125  ! dumb_mat2(i,:)=matmul(dumb_mat1,mean_full(1:ndim))
126  ! ENDDO
127  ! CALL matc_to_coo(dumb_mat2,M12)
128 
129  m12def=.not.tensor_empty(m12)
130 
131  !M13
132  m13=0.d0
134 
135  !M1tot
136  m1tot=0.d0
137  m1tot=m11+m13
138 
139  print*, "Computing the M2 terms..."
140  ALLOCATE(m21(ndim), m22(ndim), stat=allocstat)
141  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
142 
143  !M21
144  CALL copy_coo(lxy,m21)
145  CALL add_to_tensor(bxxy,m21)
146 
147  m21def=.not.tensor_empty(m21)
148 
149  !M22
150  CALL copy_coo(bxyy,m22)
151 
152  m22def=.not.tensor_empty(m22)
153 
154  !M3 tensor
155  print*, "Computing the M3 terms..."
156  ! Linear terms
157  print*, "Computing the L subterms..."
158  ALLOCATE(l1(ndim,mems), l2(ndim,mems), l4(ndim,mems), l5(ndim,mems), stat=allocstat)
159  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
160 
161  !L1
164  DO m=1,mems
165  CALL coo_to_mat_ik(dy(:,m),dumb_mat3)
166  dumb_mat4=matmul(dumb_mat2,matmul(transpose(dumb_mat3),dumb_mat1))
167  CALL matc_to_coo(dumb_mat4,l1(:,m))
168  ENDDO
169 
170  !L2
171  DO m=1,mems
172  dumb_mat4=0.d0
173  DO i=1,ndim
174  CALL coo_to_mat_i(i,bxyy,dumb_mat1)
176  DO j=1,ndim
177  CALL coo_to_mat_j(j,byxy,dumb_mat1)
178  dumb_mat4(i,j)=mat_trace(matmul(dumb_mat1,dumb_mat2))
179  ENDDO
180  END DO
181  CALL matc_to_coo(dumb_mat4,l2(:,m))
182  ENDDO
183 
184  !L4
185  ! DO m=1,mems
186  ! dumb_mat4=0.D0
187  ! DO i=1,ndim
188  ! CALL coo_to_mat_i(i,Bxyy,dumb_mat1)
189  ! CALL sparse_mul3_with_mat(dYY(:,m),dumb_mat1,dumb_vec) ! Bxyy*dYY
190  ! CALL coo_to_mat_ik(Lyx,dumb_mat1)
191  ! dumb_mat4(i,:)=matmul(transpose(dumb_mat1),dumb_vec)
192  ! ENDDO
193  ! CALL matc_to_coo(dumb_mat4,L4(:,m))
194  ! ENDDO
195 
196  !L5
197 
198  ! CALL coo_to_mat_ik(Lxy,dumb_mat1)
199  ! DO m=1,mems
200  ! dumb_mat4=0.D0
201  ! DO i=1,ndim
202  ! CALL sparse_mul3_mat(YdY(:,m),dumb_mat1(i,:),dumb_mat2)
203  ! DO j=1,ndim
204  ! CALL coo_to_mat_j(j,Byxy,dumb_mat3)
205  ! dumb_mat4(i,j)=mat_trace(matmul(dumb_mat3,dumb_mat2))
206  ! ENDDO
207  ! END DO
208  ! CALL matc_to_coo(dumb_mat4,L5(:,m))
209  ! ENDDO
210 
211  !Ltot
212 
213  ALLOCATE(ltot(ndim,mems), stat=allocstat)
214  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
215 
216  DO m=1,mems
217  CALL add_to_tensor(l1(:,m),ltot(:,m))
218  CALL add_to_tensor(l2(:,m),ltot(:,m))
219  CALL add_to_tensor(l4(:,m),ltot(:,m))
220  CALL add_to_tensor(l5(:,m),ltot(:,m))
221  ENDDO
222 
223  ldef=.not.tensor_empty(ltot)
224 
225  print*, "Computing the B terms..."
226  ALLOCATE(b1(ndim,mems), b2(ndim,mems), b3(ndim,mems), b4(ndim,mems), stat=allocstat)
227  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
228 
229  ! B1
231  dumb_mat1=transpose(dumb_mat1)
232  DO m=1,mems
233  CALL coo_to_mat_ik(dy(:,m),dumb_mat2)
235  DO j=1,ndim
236  DO k=1,ndim
237  CALL coo_to_vec_jk(j,k,byxx,dumb_vec)
238  dumb_vec=matmul(dumb_vec,dumb_mat2)
239  CALL add_vec_jk_to_tensor(j,k,dumb_vec,b1(:,m))
240  ENDDO
241  ENDDO
242  ENDDO
243 
244  ! B2
246  dumb_mat3=transpose(dumb_mat3)
247  DO m=1,mems
248  DO i=1,ndim
249  CALL coo_to_mat_i(i,bxxy,dumb_mat1)
250  CALL coo_to_mat_ik(dy(:,m),dumb_mat2)
251  dumb_mat1=matmul(dumb_mat2,transpose(dumb_mat1))
253  CALL add_matc_to_tensor(i,dumb_mat1,b2(:,m))
254  ENDDO
255  ENDDO
256 
257  ! B3
258  ! DO m=1,mems
259  ! DO i=1,ndim
260  ! CALL coo_to_mat_i(i,Bxxy,dumb_mat1)
261  ! dumb_mat4=0.D0
262  ! DO j=1,ndim
263  ! CALL coo_to_mat_j(j,YdY(:,m),dumb_mat2)
264  ! CALL coo_to_mat_i(j,Byxy,dumb_mat3)
265  ! dumb_mat2=matmul(dumb_mat3,dumb_mat2)
266  ! dumb_mat4=dumb_mat4+dumb_mat2
267  ! ENDDO
268  ! dumb_mat4=matmul(dumb_mat4,transpose(dumb_mat1))
269  ! CALL add_matc_to_tensor(i,dumb_mat4,B3(:,m))
270  ! ENDDO
271  ! ENDDO
272 
273  ! B4
274  ! DO m=1,mems
275  ! DO i=1,ndim
276  ! CALL coo_to_mat_i(i,Bxyy,dumb_mat1)
277  ! CALL sparse_mul3_with_mat(dYY(:,m),dumb_mat1,dumb_vec) ! Bxyy*dYY
278  ! DO j=1,ndim
279  ! CALL coo_to_mat_j(j,Byxx,dumb_mat1)
280  ! dumb_mat4(j,:)=matmul(transpose(dumb_mat1),dumb_vec)
281  ! ENDDO
282  ! CALL add_matc_to_tensor(i,dumb_mat4,B4(:,m))
283  ! ENDDO
284  ! ENDDO
285 
286  ALLOCATE(b14(ndim,mems), b23(ndim,mems), stat=allocstat)
287  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
288 
289  DO m=1,mems
290  CALL add_to_tensor(b1(:,m),b14(:,m))
291  CALL add_to_tensor(b2(:,m),b23(:,m))
292  CALL add_to_tensor(b4(:,m),b14(:,m))
293  CALL add_to_tensor(b3(:,m),b23(:,m))
294  ENDDO
295 
296  b14def=.not.tensor_empty(b14)
297  b23def=.not.tensor_empty(b23)
298 
299  !M
300 
301  print*, "Computing the M term..."
302 
303  ALLOCATE(mtot(ndim,mems), stat=allocstat)
304  IF (allocstat /= 0) stop "*** Not enough memory ! ***"
305 
306  DO m=1,mems
307  DO i=1,ndim
308  CALL coo_to_mat_i(i,bxxy,dumb_mat1)
309  CALL coo_to_mat_ik(dy(:,m),dumb_mat2)
310  dumb_mat1=matmul(dumb_mat2,transpose(dumb_mat1))
311  DO j=1,ndim
312  DO k=1,ndim
313  CALL coo_to_vec_jk(j,k,byxx,dumb_vec)
314  dumb_vec=matmul(dumb_vec,dumb_mat1)
315  CALL add_vec_ijk_to_tensor4(i,j,k,dumb_vec,mtot(:,m))
316  ENDDO
317  END DO
318  END DO
319  END DO
320 
321  mdef=.not.tensor4_empty(mtot)
322 
323 
324  DEALLOCATE(dumb_mat1, dumb_mat2, stat=allocstat)
325  IF (allocstat /= 0) stop "*** Problem to deallocate ! ***"
326 
327  DEALLOCATE(dumb_mat3, dumb_mat4, stat=allocstat)
328  IF (allocstat /= 0) stop "*** Problem to deallocate ! ***"
329 
330  DEALLOCATE(dumb_vec, stat=allocstat)
331  IF (allocstat /= 0) stop "*** Problem to deallocate ! ***"
332 
333 
334  END SUBROUTINE init_wl_tensor
335 
336 
337 END MODULE wl_tensor
338 
339 
340 
type(coolist), dimension(:), allocatable, public lxy
Tensor holding the linear part of the resolved tendencies involving the unresolved variables...
Definition: dec_tensor.f90:38
The stochastic models parameters module.
integer ndim
Number of variables (dimension of the model)
Definition: params.f90:85
real(kind=8), dimension(:,:), allocatable, public corr_i_full
Covariance matrix of the unresolved variables (full version)
Definition: corrmod.f90:28
Utility module.
Definition: util.f90:12
subroutine, public add_vec_ijk_to_tensor4(i, j, k, src, dst)
Routine to add a vector to a rank-4 tensor.
Definition: tensor.f90:785
logical, public b14def
Definition: WL_tensor.f90:75
subroutine, public copy_coo(src, dst)
Routine to copy a coolist.
Definition: tensor.f90:72
logical, public b23def
Definition: WL_tensor.f90:75
Statistics accumulators.
Definition: stat.f90:14
logical, public ldef
Definition: WL_tensor.f90:75
subroutine, public add_vec_jk_to_tensor(j, k, src, dst)
Routine to add a vector to a rank-3 tensor.
Definition: tensor.f90:602
The resolved-unresolved components decomposition of the tensor.
Definition: dec_tensor.f90:16
type(coolist4), dimension(:,:), allocatable, public ydyy
Coolist holding the terms.
Definition: corr_tensor.f90:34
real(kind=8), dimension(:), allocatable, public mean_full
Vector holding the mean of the unresolved dynamics (full version)
Definition: corrmod.f90:27
type(coolist), dimension(:), allocatable, public m12
Second component of the M1 term.
Definition: WL_tensor.f90:43
type(coolist), dimension(:,:), allocatable, public l4
Fourth linear tensor.
Definition: WL_tensor.f90:55
logical, public m21def
Definition: WL_tensor.f90:75
type(coolist), dimension(:,:), allocatable, public l5
Fifth linear tensor.
Definition: WL_tensor.f90:56
subroutine, public coo_to_mat_i(i, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
Definition: tensor.f90:1112
type(coolist), dimension(:,:), allocatable, public b2
Second quadratic tensor.
Definition: WL_tensor.f90:61
subroutine, public coo_to_mat_j(j, src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix.
Definition: tensor.f90:1148
Tensor utility module.
Definition: tensor.f90:18
real(kind=8) function, public mat_trace(A)
Definition: util.f90:150
real(kind=8), dimension(:), allocatable, public m13
Third component of the M1 term.
Definition: WL_tensor.f90:44
type(coolist), dimension(:,:), allocatable, public b14
Joint 1st and 4th tensors.
Definition: WL_tensor.f90:64
real(kind=8) function, public mat_contract(A, B)
Definition: util.f90:162
logical, public m12def
Definition: WL_tensor.f90:75
subroutine, public add_matc_to_tensor(i, src, dst)
Routine to add a matrix to a rank-3 tensor.
Definition: tensor.f90:474
Coordinate list. Type used to represent the sparse tensor.
Definition: tensor.f90:38
type(coolist), dimension(:,:), allocatable, public dy
Coolist holding the terms.
Definition: corr_tensor.f90:31
logical function, public tensor_empty(t)
Test if a rank-3 tensor coolist is empty.
Definition: tensor.f90:1288
type(coolist), dimension(:), allocatable, public m22
Second tensor of the M2 term.
Definition: WL_tensor.f90:49
type(coolist), dimension(:), allocatable, public byxy
Tensor holding the quadratic part of the unresolved tendencies involving both resolved and unresolved...
Definition: dec_tensor.f90:47
subroutine, public sparse_mul4_with_mat_kl(coolist_ijkl, mat_kl, res)
Sparse multiplication of a rank-4 tensor coolist with a matrix : .
Definition: tensor.f90:1194
subroutine, public init_corr_tensor
Subroutine to initialise the correlations tensors.
Definition: corr_tensor.f90:45
type(coolist), dimension(:,:), allocatable, public b3
Third quadratic tensor.
Definition: WL_tensor.f90:62
The WL tensors used to integrate the model.
Definition: WL_tensor.f90:18
type(coolist4), dimension(:,:), allocatable, public mtot
Tensor for the cubic terms.
Definition: WL_tensor.f90:67
real(kind=8), dimension(:), allocatable dumb_vec
Dumb vector to be used in the calculation.
Definition: corr_tensor.f90:36
type(coolist), dimension(:,:), allocatable, public l1
First linear tensor.
Definition: WL_tensor.f90:53
subroutine, public coo_to_mat_ik(src, dst)
Routine to convert a rank-3 tensor coolist component into a matrix with i and k indices.
Definition: tensor.f90:1063
subroutine, public add_to_tensor(src, dst)
Routine to add a rank-3 tensor to another one.
Definition: tensor.f90:350
type(coolist), dimension(:,:), allocatable, public l2
Second linear tensor.
Definition: WL_tensor.f90:54
real(kind=8), dimension(:,:), allocatable dumb_mat4
Dummy matrix.
Definition: WL_tensor.f90:73
type(coolist), dimension(:), allocatable, public bxxy
Tensor holding the quadratic part of the resolved tendencies involving both resolved and unresolved v...
Definition: dec_tensor.f90:40
real(kind=8), dimension(:,:), allocatable dumb_mat3
Dummy matrix.
Definition: WL_tensor.f90:72
4d coordinate list. Type used to represent the rank-4 sparse tensor.
Definition: tensor.f90:44
subroutine, public matc_to_coo(src, dst)
Routine to convert a matrix to a rank-3 tensor.
Definition: tensor.f90:1244
subroutine, public sparse_mul3_with_mat(coolist_ijk, mat_jk, res)
Sparse multiplication of a rank-3 tensor coolist with a matrix: .
Definition: tensor.f90:1220
real(kind=8), dimension(:,:), allocatable dumb_mat2
Dumb matrix to be used in the calculation.
Definition: corr_tensor.f90:38
Module to initialize the correlation matrix of the unresolved variables.
Definition: corrmod.f90:16
real(kind=8), dimension(:,:), allocatable dumb_mat1
Dumb matrix to be used in the calculation.
Definition: corr_tensor.f90:37
The model parameters module.
Definition: params.f90:18
type(coolist), dimension(:), allocatable, public bxyy
Tensor holding the quadratic part of the resolved tendencies involving unresolved variables...
Definition: dec_tensor.f90:41
type(coolist), dimension(:,:), allocatable, public b23
Joint 2nd and 3rd tensors.
Definition: WL_tensor.f90:65
type(coolist), dimension(:,:), allocatable, public b4
Fourth quadratic tensor.
Definition: WL_tensor.f90:63
logical, public mdef
Boolean to (de)activate the computation of the terms.
Definition: WL_tensor.f90:75
logical, public m22def
Definition: WL_tensor.f90:75
logical function, public tensor4_empty(t)
Test if a rank-4 tensor coolist is empty.
Definition: tensor.f90:1304
type(coolist), dimension(:), allocatable, public lyx
Tensor holding the linear part of the unresolved tendencies involving the resolved variables...
Definition: dec_tensor.f90:44
type(coolist), dimension(:,:), allocatable, public ltot
Total linear tensor.
Definition: WL_tensor.f90:57
type(coolist), dimension(:), allocatable, public byxx
Tensor holding the quadratic part of the unresolved tendencies involving resolved variables...
Definition: dec_tensor.f90:46
Module to compute the correlations and derivatives used to compute the memory term of the WL paramete...
Definition: corr_tensor.f90:17
subroutine, public init_wl_tensor
Subroutine to initialise the WL tensor.
Definition: WL_tensor.f90:94
subroutine, public init_dec_tensor
Subroutine that initialize and compute the decomposed tensors.
Definition: dec_tensor.f90:195
real(kind=8), dimension(:), allocatable, public m11
First component of the M1 term.
Definition: WL_tensor.f90:42
type(coolist), dimension(:), allocatable, public m21
First tensor of the M2 term.
Definition: WL_tensor.f90:48
subroutine, public coo_to_vec_jk(j, k, src, dst)
Routine to convert a rank-3 tensor coolist component into a vector.
Definition: tensor.f90:1129
type(coolist), dimension(:,:), allocatable, public b1
First quadratic tensor.
Definition: WL_tensor.f90:60
real(kind=8), dimension(:), allocatable, public m1tot
Total vector.
Definition: WL_tensor.f90:45
integer mems
Number of steps in the memory kernel integral.