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