MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_simplex_noise.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/pre_process/m_simplex_noise.fpp"
2!>
3!! @file
4!! @brief Contains module m_simplex_noise
5
6!> @brief 2D and 3D simplex noise generation for procedural initial condition perturbations
7module m_simplex_noise
8
9 use m_constants
10
11 use m_precision_select
12
13 implicit none
14
15 private; public :: f_simplex3d, &
16 f_simplex2d
17
18 integer, parameter :: p_vec(0:511) = [ &
19 151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, &
20 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, &
21 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, &
22 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, &
23 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, &
24 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, &
25 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, &
26 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, &
27 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9, &
28 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, &
29 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, &
30 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, &
31 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, &
32 61, 156, 180, &
33 151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, &
34 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, &
35 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, &
36 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, &
37 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, &
38 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, &
39 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, &
40 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, &
41 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9, &
42 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, &
43 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, &
44 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, &
45 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, &
46 61, 156, 180]
47
48 real(wp), parameter :: grad3(12, 3) = reshape([ &
49 1._wp, 1._wp, 0._wp, &
50 -1._wp, 1._wp, 0._wp, &
51 1._wp, -1._wp, 0._wp, &
52 -1._wp, -1._wp, 0._wp, &
53 1._wp, 0._wp, 1._wp, &
54 -1._wp, 0._wp, 1._wp, &
55 1._wp, 0._wp, -1._wp, &
56 -1._wp, 0._wp, -1._wp, &
57 0._wp, 1._wp, 1._wp, &
58 0._wp, -1._wp, 1._wp, &
59 0._wp, 1._wp, -1._wp, &
60 0._wp, -1._wp, -1._wp], shape=[12, 3])
61
62 real(wp), parameter :: grad2(10, 2) = reshape([ &
63 1._wp, 1._wp, &
64 -1._wp, 1._wp, &
65 1._wp, -1._wp, &
66 -1._wp, -1._wp, &
67 1._wp, 0._wp, &
68 -1._wp, 0._wp, &
69 0._wp, 1._wp, &
70 0._wp, -1._wp, &
71 1._wp, 1._wp, &
72 -1._wp, 1._wp], shape=[10, 2])
73
74contains
75
76 !> @brief Evaluates 3D simplex noise at the given coordinates and returns a value in [-1, 1].
77 function f_simplex3d(xin, yin, zin) result(n)
78
79 real(wp), intent(in) :: xin, yin, zin
80 real(wp) :: n
81 real(wp) :: n0, n1, n2, n3
82 real(wp) :: f3, g3
83 real(wp) :: x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3
84 integer :: i, j, k, i1, j1, k1, i2, j2, k2
85 integer :: ii, jj, kk, gi0, gi1, gi2, gi3
86 real(wp) :: s, t, r, t0, t1, t2, t3
87 real(wp) :: g(3)
88 real(wp) :: x, y, z
89
90 f3 = 1._wp/3._wp
91 g3 = 1._wp/6._wp
92
93 s = (xin + yin + zin)*f3
94 i = floor(xin + s)
95 j = floor(yin + s)
96 k = floor(zin + s)
97
98 t = (i + j + k)*g3
99
100 x0 = xin - (i - t)
101 y0 = yin - (j - t)
102 z0 = zin - (k - t)
103
104 if (x0 >= y0) then
105 if (y0 >= z0) then
106 i1 = 1; j1 = 0; k1 = 0; i2 = 1; j2 = 1; k2 = 0
107 else if (x0 >= z0) then
108 i1 = 1; j1 = 0; k1 = 0; i2 = 1; j2 = 0; k2 = 1
109 else
110 i1 = 0; j1 = 0; k1 = 1; i2 = 1; j2 = 0; k2 = 1
111 end if
112 else
113 if (y0 < z0) then
114 i1 = 0; j1 = 0; k1 = 1; i2 = 0; j2 = 1; k2 = 1
115 else if (x0 < z0) then
116 i1 = 0; j1 = 1; k1 = 0; i2 = 0; j2 = 1; k2 = 1
117 else
118 i1 = 0; j1 = 1; k1 = 0; i2 = 1; j2 = 1; k2 = 0
119 end if
120 end if
121
122 x1 = x0 - i1 + g3
123 y1 = y0 - j1 + g3
124 z1 = z0 - k1 + g3
125 x2 = x0 - i2 + 2._wp*g3
126 y2 = y0 - j2 + 2._wp*g3
127 z2 = z0 - k2 + 2._wp*g3
128 x3 = x0 - 1._wp + 3._wp*g3
129 y3 = y0 - 1._wp + 3._wp*g3
130 z3 = z0 - 1._wp + 3._wp*g3
131
132 ii = iand(i, 255)
133 jj = iand(j, 255)
134 kk = iand(k, 255)
135
136 gi0 = mod(p_vec(ii + p_vec(jj + p_vec(kk) + 1) + 1), 12) + 1
137 gi1 = mod(p_vec(ii + i1 + p_vec(jj + j1 + p_vec(kk + k1) + 1) + 1), 12) + 1
138 gi2 = mod(p_vec(ii + i2 + p_vec(jj + j2 + p_vec(kk + k2) + 1) + 1), 12) + 1
139 gi3 = mod(p_vec(ii + 1 + p_vec(jj + 1 + p_vec(kk + 1) + 1) + 1), 12) + 1
140
141 t0 = 0.5_wp - x0*x0 - y0*y0 - z0*z0
142 if (t0 < 0._wp) then
143 n0 = 0._wp
144 else
145 t0 = t0*t0
146 n0 = t0*t0*dot_product(grad3(gi0, :), [x0, y0, z0])
147 end if
148
149 t1 = 0.5_wp - x1*x1 - y1*y1 - z1*z1
150 if (t1 < 0._wp) then
151 n1 = 0._wp
152 else
153 t1 = t1*t1
154 n1 = t1*t1*dot_product(grad3(gi1, :), [x1, y1, z1])
155 end if
156
157 t2 = 0.5_wp - x2*x2 - y2*y2 - z2*z2
158 if (t2 < 0._wp) then
159 n2 = 0._wp
160 else
161 t2 = t2*t2
162 n2 = t2*t2*dot_product(grad3(gi2, :), [x2, y2, z2])
163 end if
164
165 t3 = 0.5_wp - x3*x3 - y3*y3 - z3*z3
166 if (t3 < 0._wp) then
167 n3 = 0._wp
168 else
169 t3 = t3*t3
170 n3 = t3*t3*dot_product(grad3(gi3, :), [x3, y3, z3])
171 end if
172
173 n = 32._wp*(n0 + n1 + n2 + n3)
174
175 end function f_simplex3d
176
177 !> @brief Evaluates 2D simplex noise at the given coordinates and returns a value in [-1, 1].
178 function f_simplex2d(xin, yin) result(n)
179
180 real(wp), intent(in) :: xin, yin
181 real(wp) :: n
182 real(wp), parameter :: f2 = 0.5_wp*(sqrt(3._wp) - 1._wp)
183 real(wp), parameter :: g2 = (3._wp - sqrt(3._wp))/6._wp
184 integer :: i, j, ii, jj, gi0, gi1, gi2
185 real(wp) :: s, t, x0, y0, x1, y1, x2, y2
186 real(wp) :: t0, t1, t2, n0, n1, n2
187 integer :: i1, j1
188
189 s = (xin + yin)*f2
190 i = floor(xin + s)
191 j = floor(yin + s)
192
193 t = real(i + j, 8)*g2
194
195 x0 = xin - (i - t)
196 y0 = yin - (j - t)
197
198 if (x0 > y0) then
199 i1 = 1; j1 = 0
200 else
201 i1 = 0; j1 = 1
202 end if
203
204 x1 = x0 - i1 + g2
205 y1 = y0 - j1 + g2
206 x2 = x0 - 1._wp + 2._wp*g2
207 y2 = y0 - 1._wp + 2._wp*g2
208
209 ii = mod(i, 255)
210 jj = mod(j, 255)
211
212 gi0 = mod(p_vec(ii + p_vec(jj)), 10) + 1
213 gi1 = mod(p_vec(ii + i1 + p_vec(jj + j1)), 10) + 1
214 gi2 = mod(p_vec(ii + 1 + p_vec(jj + 1)), 10) + 1
215
216 t0 = 0.5_wp - x0*x0 - y0*y0
217 if (t0 < 0._wp) then
218 n0 = 0._wp
219 else
220 t0 = t0*t0
221 n0 = t0*t0*dot2(gi0, x0, y0)
222 end if
223
224 t1 = 0.5_wp - x1*x1 - y1*y1
225 if (t1 < 0._wp) then
226 n1 = 0._wp
227 else
228 t1 = t1*t1
229 n1 = t1*t1*dot2(gi1, x1, y1)
230 end if
231
232 t2 = 0.5_wp - x2*x2 - y2*y2
233 if (t2 < 0._wp) then
234 n2 = 0._wp
235 else
236 t2 = t2*t2
237 n2 = t2*t2*dot2(gi2, x2, y2)
238 end if
239
240 n = 70._wp*(n0 + n1 + n2)
241
242 end function f_simplex2d
243
244 !> @brief Computes the dot product of a 2D gradient vector with the given offset coordinates.
245 function dot2(g, x, y) result(dot)
246
247 integer, intent(in) :: g
248 real(wp), intent(in) :: x, y
249 real(wp) :: dot
250 dot = grad2(g + 1, 1)*x + grad2(g + 1, 2)*y
251
252 end function
253
254end module m_simplex_noise
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1121
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1121
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
m_precision_select
Working-precision kind selection (half/single/double) and corresponding MPI datatype parameters.
Definition m_precision_select.f90:5
m_precision_select::wp
integer, parameter wp
Definition m_precision_select.f90:26
m_simplex_noise
2D and 3D simplex noise generation for procedural initial condition perturbations
Definition m_simplex_noise.fpp.f90:7
m_simplex_noise::p_vec
integer, dimension(0:511), parameter p_vec
Definition m_simplex_noise.fpp.f90:18
m_simplex_noise::f_simplex2d
real(wp) function, public f_simplex2d(xin, yin)
Evaluates 2D simplex noise at the given coordinates and returns a value in [-1, 1].
Definition m_simplex_noise.fpp.f90:179
m_simplex_noise::grad3
real(wp), dimension(12, 3), parameter grad3
Definition m_simplex_noise.fpp.f90:48
m_simplex_noise::grad2
real(wp), dimension(10, 2), parameter grad2
Definition m_simplex_noise.fpp.f90:62
m_simplex_noise::dot2
real(wp) function dot2(g, x, y)
Computes the dot product of a 2D gradient vector with the given offset coordinates.
Definition m_simplex_noise.fpp.f90:246
m_simplex_noise::f_simplex3d
real(wp) function, public f_simplex3d(xin, yin, zin)
Evaluates 3D simplex noise at the given coordinates and returns a value in [-1, 1].
Definition m_simplex_noise.fpp.f90:78