-
Notifications
You must be signed in to change notification settings - Fork 6
/
mo_random_field.f90
289 lines (249 loc) · 13.5 KB
/
mo_random_field.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
!> \file mo_random_field.f90
!> \brief This module contains generating random fields with certain
!> statistical properties, e.g. correlation length.
!> \details This module contains generating random fields with certain
!> statistical properties, e.g. correlation length.\n
!> Right now there are the following implemented:
!> * Irrotational flow fields with Gauss-shape correlation functions
!> described by Attinger et al. (2008) 'Effective velocity for transport in
!> heterogeneous compressible flows with mean drift' (Eq. 31-32)
!> * some more...
!> \authors Juliane Mai
!> \date Jul 2014
module mo_random_field
! Written Juliane Mai, Jul 2014
! License
! -------
! This file is part of the JAMS Fortran package, distributed under the MIT License.
!
! Copyright (c) 2014 Juliane Mai
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
use mo_kind, only: i4, i8, dp
implicit none
public :: random_velocity_field_gauss ! Fields with gauss-shape correlation function
! ------------------------------------------------------------------
! NAME
! random_velocity_field_gauss
! PURPOSE
!> \brief Calculates random velocity fields by using a correlation function
!> which is a product of Gauss-shaped functions as described
!> by Attinger et al. (2008).
!
!> \details Calculates random velocity fields by using a correlation function
!> which is a product of Gauss-shaped functions as described
!> by Attinger et al. (2008).\n
!> A transport equation in a reference domain is considered where the
!> velocity field \f$ u \f$ splits into \f$ u = \bar{u} + \tilde{u}\f$.
!> \f$ \bar{u}\f$ is a constant mean and \f$ \tilde{u} \f$ is a random
!> compressible velocity field defined by a zero mean random Gaussian
!> correlation function.\n
!> The potential can be generated by a superposition of randomly
!> chosen cosine modes as follows:
!> \f[ u_i^N = \sigma \sqrt{\frac{2}{N}} \sum\limits_{j=1}^{N} k_i^j \cos(k^jx + \omega_j) \f]
!> where each component of \f$ k^j=(k_1^j,k_2^j)\f$ obeys a Gaussian
!> distribution with zero mean and variance given by \f$ 1/l_i^2 \f$.
!> \f$ l_i \f$ are the correlation lengths, \f$ N \f$ is the number of cosine modes used for approximation.
!
! INTENT(IN)
!> \param[in] "real(sp/dp) :: coord(n,m)" 2D-array of n coordinates in m dimensional space
!> where field has to be evaluated
!> \param[in] "real(sp/dp) :: corr_length(m)" correlation length for each of the m dimensions
!
! INTENT(INOUT)
! None
!
! INTENT(OUT)
! None
!
! INTENT(IN), OPTIONAL
!> \param[in] "real(sp/dp) :: sigma2" variance of correlation function\n
!> DEFAULT: 1.0_dp
!> \param[in] "integer(i4), optional :: ncosinemode" number of cosine modes summed up\n
!> DEFAULT: 100
!
! INTENT(INOUT), OPTIONAL
! None
!
! INTENT(OUT), OPTIONAL
!> \param[out] "real(sp/dp), optional :: cosine_modes(ncosinemode,m)"
!> returns <ncosinemode> random Gaussian modes
!> \f$ k^j=(k_1^j,...,k_m^j)\f$, i.e. all
!> information describing the field.
!> Therefore, additional positions
!> coord(:,:) can be evaluated afterwards.
!> \param[out] "real(sp/dp), optional :: potential(n)" returns the potential of the field at all
!> requested points
!
! RETURN
!> \return real(sp/dp) :: velocity(n,m) — velocity at given coordinates
!
! RESTRICTIONS
!> \note None.
!
! EXAMPLE
! -> see also example in test directory
! LITERATURE
! Attinger, S., & Abdulle, A. (2008).
! Effective velocity for transport in heterogeneous compressible flows with mean drift.
! Physics of Fluids
! HISTORY
!> \author Juliane Mai
!> \date Jul 2014
! Modified,
interface random_velocity_field_gauss
module procedure random_velocity_field_gauss_dp !, random_velocity_field_gauss_sp
end interface random_velocity_field_gauss
! ------------------------------------------------------------------
private
contains
function random_velocity_field_gauss_dp(coord, corr_length, ncosinemode, sigma2, seed, &
cosine_modes_in, cosine_modes_out, potential) result (velocity)
use mo_constants, only: pi_dp
use mo_xor4096, only: xor4096g, get_timeseed
use mo_xor4096_apps, only: xor4096g_array
implicit none
real(dp), dimension(:,:), intent(in) :: coord ! coordinates where field has to be evaluated
! ! dim1: number of requested points
! ! dim2: number of dimensions
real(dp), dimension(:), intent(in) :: corr_length ! correlation length for each dimension
! ! dim1: number of dimensions
integer(i4), intent(in), optional :: ncosinemode ! number of cosine modes summed up
! ! DEFAULT: 100
real(dp), intent(in), optional :: sigma2 ! variance of correlation function
! ! DEFAULT: 1.0_dp
integer(i8), intent(in), optional :: seed ! seed for random numbers
! ! DEFAULT: -9_i8 --> time_seed
real(dp), dimension(:,:), intent(in), optional :: cosine_modes_in ! cosine mode is already determined in prev.
! ! run and is now fixed there, i.e. field
! ! generation is skipped and only velocity
! ! and potential are determined
! ! dim1: number of cosine modes
! ! dim2: number of dimensions
real(dp), dimension(:,:), allocatable, intent(out), optional :: cosine_modes_out ! returns the random cosine modes which
! ! define the field
! ! dim1: number of cosine modes
! ! dim2: number of dimensions
real(dp), dimension(size(coord,1)), intent(out), optional :: potential ! potential of field at given point
! ! dim1: number of requested points
real(dp), dimension(size(coord,1),size(coord,2)) :: velocity ! velocity at given point
! ! dim1: number of requested points
! ! dim2: number of dimensions
! local variables
integer(i4) :: idim, ipoint, imode ! counter
real(dp) :: h, h_pot ! helper variables
real(dp) :: isigma2 ! variance of correlation function
integer(i4) :: incosinemode ! number of cosine modes summed up
integer(i4) :: npoints ! number of points given
integer(i4) :: dims ! number of dimensions
real(dp), dimension(size(coord,1)) :: ipotential ! potential of field at given point
real(dp), dimension(:,:), allocatable :: modes ! random cosine modes which define the field
integer(i8), dimension(:), allocatable :: iseed ! seed for random numbers
if (size(coord,2) .ne. size(corr_length,1)) then
stop 'random_velocity_field_gauss: correlation length for each dimension has to be given'
end if
if ( any(corr_length .lt. 0.0_dp) ) then
stop 'random_velocity_field_gauss: correlation length has to be positive'
end if
if (present(ncosinemode)) then
if ( ncosinemode .le. 0 ) then
stop 'random_velocity_field_gauss: ncosinemode is the number of cosine modes summed up and therefore has to be positive'
else
incosinemode = ncosinemode
end if
else
incosinemode = 100
end if
if (present(cosine_modes_in)) then
if ( size(cosine_modes_in,2) .ne. size(coord,2) ) then
stop 'random_velocity_field_gauss: shape mismatch in second dimension'
end if
incosinemode = size(cosine_modes_in,1)
! later: values handed over to "modes"
end if
if (present(sigma2)) then
if ( sigma2 .lt. 0.0_dp ) then
stop 'random_velocity_field_gauss: sigma2 is variance of correlation function and therefore has to be positive'
else
isigma2 = sigma2
end if
else
isigma2 = 1.0_dp
end if
npoints = size(coord,1)
dims = size(coord,2)
allocate(modes(incosinemode, dims))
modes = 0.0_dp
set_modes: if (present(cosine_modes_in)) then
modes = cosine_modes_in
else
allocate(iseed(dims))
if (present(seed)) then
if (seed .gt. 0_i8) then
iseed(1) = seed
do idim=2, dims
iseed(idim) = iseed(idim-1) + 1000_i8
end do
else
call get_timeseed(iseed)
end if
else
call get_timeseed(iseed)
end if
! initialize Random numbers
call xor4096g(iseed, modes(1,:))
! generate random modes
call xor4096g_array(modes)
! scale with correlation length
do idim=1, dims
modes(:,idim) = modes(:,idim) * 1.0_dp/corr_length(idim)
end do
deallocate(iseed)
end if set_modes
velocity(:,:) = 0.0_dp
ipotential(:) = 0.0_dp
do ipoint=1, npoints
do imode=1, incosinemode
h = 0.0_dp ! velocity
h_pot = 0.0_dp ! potential
do idim=1, dims
h = h + modes(imode,idim) * coord(ipoint,idim)
end do
h_pot = h
h = cos( h + real(imode,dp)*2.0_dp*pi_dp/real(incosinemode,dp) )
h_pot = sin( h_pot + real(imode,dp)*2.0_dp*pi_dp/real(incosinemode,dp) )
velocity(ipoint,:) = velocity(ipoint,:) + h * modes(imode,:)
ipotential(ipoint) = ipotential(ipoint) + h_pot
end do
end do
do idim=1,dims
velocity(:,idim) = velocity(:,idim) * (isigma2) * sqrt( 2.0_dp/real(incosinemode,dp) )
end do
ipotential(:) = ipotential(:) * (isigma2) * sqrt( 2.0_dp/real(incosinemode,dp) )
if (present(cosine_modes_out)) then
if (allocated(cosine_modes_out)) deallocate(cosine_modes_out)
allocate(cosine_modes_out(incosinemode, dims))
cosine_modes_out = modes
end if
if (present(potential)) then
potential = ipotential
end if
deallocate(modes)
end function random_velocity_field_gauss_dp
end module mo_random_field