Neko 1.99.9
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
ax_helm_svv_one_sided_full_cpu.f90
Go to the documentation of this file.
1! Copyright (c) 2025-2026, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
34module ax_helm_svv_one_sided_full_cpu
35 use ax_helm_svv_full, only : ax_helm_svv_full_t
36 use num_types, only : rp
37 use coefs, only : coef_t
38 use space, only : space_t
39 use mesh, only : mesh_t
40 use math, only : addcol4
41 use tensor, only : tnsr3d_el, tnsr3d
42 use spectral_vanishing_viscosity, only : svv_t
43 implicit none
44 private
45
47 type, public, extends(ax_helm_svv_full_t) :: ax_helm_svv_one_sided_full_cpu_t
48 contains
50 procedure, pass(this) :: compute_vector => &
51 ax_helm_svv_one_sided_full_compute_vector
52 end type ax_helm_svv_one_sided_full_cpu_t
53
54contains
55
67 subroutine ax_helm_svv_one_sided_full_compute_vector(this, au, av, aw, &
68 u, v, w, coef, msh, Xh)
69 class(ax_helm_svv_one_sided_full_cpu_t), intent(in) :: this
70 type(mesh_t), intent(in) :: msh
71 type(space_t), intent(in) :: Xh
72 type(coef_t), intent(in) :: coef
73 real(kind=rp), intent(in) :: u(xh%lx, xh%ly, xh%lz, msh%nelv)
74 real(kind=rp), intent(in) :: v(xh%lx, xh%ly, xh%lz, msh%nelv)
75 real(kind=rp), intent(in) :: w(xh%lx, xh%ly, xh%lz, msh%nelv)
76 real(kind=rp), intent(inout) :: au(xh%lx, xh%ly, xh%lz, msh%nelv)
77 real(kind=rp), intent(inout) :: av(xh%lx, xh%ly, xh%lz, msh%nelv)
78 real(kind=rp), intent(inout) :: aw(xh%lx, xh%ly, xh%lz, msh%nelv)
79
80 call ax_helm_svv_one_sided_full_lx(au, av, aw, u, v, w, &
81 xh%dx, xh%dy, xh%dz, xh%dxt, xh%dyt, xh%dzt, &
82 coef%h1, coef%drdx, coef%drdy, coef%drdz, coef%dsdx, coef%dsdy, &
83 coef%dsdz, coef%dtdx, coef%dtdy, coef%dtdz, &
84 coef%jacinv, xh%w3, this%svv%h1, this%svv%filter%fh, &
85 this%svv%filter%fht, this%svv%direction, this%svv%ident, &
86 msh%nelv, xh%lx)
87
88 if (coef%ifh2) then
89 call addcol4 (au, coef%h2, coef%B, u, coef%dof%size())
90 call addcol4 (av, coef%h2, coef%B, v, coef%dof%size())
91 call addcol4 (aw, coef%h2, coef%B, w, coef%dof%size())
92 end if
93
94
95 end subroutine ax_helm_svv_one_sided_full_compute_vector
96
129 subroutine ax_helm_svv_one_sided_full_lx(au, av, aw, u, v, w, &
130 Dx, Dy, Dz, Dxt, Dyt, Dzt, &
131 h1, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, &
132 jacinv, weights3, svv_h1, svv_Q, svv_Qt, svv_direction, ident, n, lx)
133 integer, intent(in) :: n, lx
134 real(kind=rp), intent(in) :: u(lx, lx, lx, n)
135 real(kind=rp), intent(in) :: v(lx, lx, lx, n)
136 real(kind=rp), intent(in) :: w(lx, lx, lx, n)
137 real(kind=rp), intent(inout) :: au(lx, lx, lx, n)
138 real(kind=rp), intent(inout) :: av(lx, lx, lx, n)
139 real(kind=rp), intent(inout) :: aw(lx, lx, lx, n)
140 real(kind=rp), intent(in) :: h1(lx, lx, lx, n)
141 real(kind=rp), intent(in) :: drdx(lx, lx, lx, n)
142 real(kind=rp), intent(in) :: drdy(lx, lx, lx, n)
143 real(kind=rp), intent(in) :: drdz(lx, lx, lx, n)
144 real(kind=rp), intent(in) :: dsdx(lx, lx, lx, n)
145 real(kind=rp), intent(in) :: dsdy(lx, lx, lx, n)
146 real(kind=rp), intent(in) :: dsdz(lx, lx, lx, n)
147 real(kind=rp), intent(in) :: dtdx(lx, lx, lx, n)
148 real(kind=rp), intent(in) :: dtdy(lx, lx, lx, n)
149 real(kind=rp), intent(in) :: dtdz(lx, lx, lx, n)
150 real(kind=rp), intent(in) :: jacinv(lx, lx, lx, n)
151 real(kind=rp), intent(in) :: weights3(lx, lx, lx)
152 real(kind=rp), intent(in) :: dx(lx,lx)
153 real(kind=rp), intent(in) :: dy(lx,lx)
154 real(kind=rp), intent(in) :: dz(lx,lx)
155 real(kind=rp), intent(in) :: dxt(lx,lx)
156 real(kind=rp), intent(in) :: dyt(lx,lx)
157 real(kind=rp), intent(in) :: dzt(lx,lx)
158 real(kind=rp), intent(in) :: svv_h1(lx, lx, lx, n)
159 real(kind=rp), intent(in) :: svv_q(lx, lx), svv_qt(lx, lx)
160 character(len=*), intent(in) :: svv_direction
161 real(kind=rp), intent(in) :: ident(lx, lx)
162
163 real(kind=rp) :: s11_h, s22_h, s33_h, s12_h, s13_h, s23_h
164
165 real(kind=rp) :: wur(lx, lx, lx)
166 real(kind=rp) :: wus(lx, lx, lx)
167 real(kind=rp) :: wut(lx, lx, lx)
168 real(kind=rp) :: wvr(lx, lx, lx)
169 real(kind=rp) :: wvs(lx, lx, lx)
170 real(kind=rp) :: wvt(lx, lx, lx)
171 real(kind=rp) :: wwr(lx, lx, lx)
172 real(kind=rp) :: wws(lx, lx, lx)
173 real(kind=rp) :: wwt(lx, lx, lx)
174
175 real(kind=rp) :: s11_svv(lx, lx, lx)
176 real(kind=rp) :: s22_svv(lx, lx, lx)
177 real(kind=rp) :: s33_svv(lx, lx, lx)
178 real(kind=rp) :: s12_svv(lx, lx, lx)
179 real(kind=rp) :: s13_svv(lx, lx, lx)
180 real(kind=rp) :: s23_svv(lx, lx, lx)
181 integer :: e, i, j, k, l
182
183 real(kind=rp) :: t1, t2, t3
184 real(kind=rp) :: s11(lx, lx, lx)
185 real(kind=rp) :: s22(lx, lx, lx)
186 real(kind=rp) :: s33(lx, lx, lx)
187 real(kind=rp) :: s12(lx, lx, lx)
188 real(kind=rp) :: s13(lx, lx, lx)
189 real(kind=rp) :: s23(lx, lx, lx)
190 real(kind=rp) :: u1, u2, u3, v1, v2, v3, w1, w2, w3
191 real(kind=rp) :: filter_r(lx, lx), filter_s(lx, lx), filter_t(lx, lx)
192
193 if (index(svv_direction, "r") > 0) then
194 filter_r = svv_q
195 else
196 filter_r = ident
197 end if
198 if (index(svv_direction, "s") > 0) then
199 filter_s = svv_qt
200 else
201 filter_s = ident
202 end if
203 if (index(svv_direction, "t") > 0) then
204 filter_t = svv_qt
205 else
206 filter_t = ident
207 end if
208
209 do e = 1, n
210 do j = 1, lx * lx
211 do i = 1, lx
212 t1 = 0.0_rp
213 t2 = 0.0_rp
214 t3 = 0.0_rp
215 do k = 1, lx
216 t1 = t1 + dx(i,k) * u(k,j,1,e)
217 t2 = t2 + dx(i,k) * v(k,j,1,e)
218 t3 = t3 + dx(i,k) * w(k,j,1,e)
219 end do
220 wur(i,j,1) = t1
221 wvr(i,j,1) = t2
222 wwr(i,j,1) = t3
223 end do
224 end do
225
226 do k = 1, lx
227 do j = 1, lx
228 do i = 1, lx
229 t1 = 0.0_rp
230 t2 = 0.0_rp
231 t3 = 0.0_rp
232 do l = 1, lx
233 t1 = t1 + dy(j,l) * u(i,l,k,e)
234 t2 = t2 + dy(j,l) * v(i,l,k,e)
235 t3 = t3 + dy(j,l) * w(i,l,k,e)
236 end do
237 wus(i,j,k) = t1
238 wvs(i,j,k) = t2
239 wws(i,j,k) = t3
240 end do
241 end do
242 end do
243
244 do k = 1, lx
245 do i = 1, lx*lx
246 t1 = 0.0_rp
247 t2 = 0.0_rp
248 t3 = 0.0_rp
249 do l = 1, lx
250 t1 = t1 + dz(k,l) * u(i,1,l,e)
251 t2 = t2 + dz(k,l) * v(i,1,l,e)
252 t3 = t3 + dz(k,l) * w(i,1,l,e)
253 end do
254 wut(i,1,k) = t1
255 wvt(i,1,k) = t2
256 wwt(i,1,k) = t3
257 end do
258 end do
259
260 do i = 1, lx*lx*lx
261 u1 = (drdx(i,1,1,e) * wur(i,1,1) &
262 + dsdx(i,1,1,e) * wus(i,1,1) &
263 + dtdx(i,1,1,e) * wut(i,1,1)) * jacinv(i,1,1,e)
264 u2 = (drdy(i,1,1,e) * wur(i,1,1) &
265 + dsdy(i,1,1,e) * wus(i,1,1) &
266 + dtdy(i,1,1,e) * wut(i,1,1)) * jacinv(i,1,1,e)
267 u3 = (drdz(i,1,1,e) * wur(i,1,1) &
268 + dsdz(i,1,1,e) * wus(i,1,1) &
269 + dtdz(i,1,1,e) * wut(i,1,1)) * jacinv(i,1,1,e)
270 v1 = (drdx(i,1,1,e) * wvr(i,1,1) &
271 + dsdx(i,1,1,e) * wvs(i,1,1) &
272 + dtdx(i,1,1,e) * wvt(i,1,1)) * jacinv(i,1,1,e)
273 v2 = (drdy(i,1,1,e) * wvr(i,1,1) &
274 + dsdy(i,1,1,e) * wvs(i,1,1) &
275 + dtdy(i,1,1,e) * wvt(i,1,1)) * jacinv(i,1,1,e)
276 v3 = (drdz(i,1,1,e) * wvr(i,1,1) &
277 + dsdz(i,1,1,e) * wvs(i,1,1) &
278 + dtdz(i,1,1,e) * wvt(i,1,1)) * jacinv(i,1,1,e)
279 w1 = (drdx(i,1,1,e) * wwr(i,1,1) &
280 + dsdx(i,1,1,e) * wws(i,1,1) &
281 + dtdx(i,1,1,e) * wwt(i,1,1)) * jacinv(i,1,1,e)
282 w2 = (drdy(i,1,1,e) * wwr(i,1,1) &
283 + dsdy(i,1,1,e) * wws(i,1,1) &
284 + dtdy(i,1,1,e) * wwt(i,1,1)) * jacinv(i,1,1,e)
285 w3 = (drdz(i,1,1,e) * wwr(i,1,1) &
286 + dsdz(i,1,1,e) * wws(i,1,1) &
287 + dtdz(i,1,1,e) * wwt(i,1,1)) * jacinv(i,1,1,e)
288 s11(i,1,1) = u1 + u1
289 s22(i,1,1) = v2 + v2
290 s33(i,1,1) = w3 + w3
291 s12(i,1,1) = u2 + v1
292 s13(i,1,1) = u3 + w1
293 s23(i,1,1) = v3 + w2
294 end do
295
296 ! spatial convolution for spectral vanishing (low pass filter (LPF))
297 call tnsr3d_el(s11_svv, lx, s11, lx, filter_r, filter_s, filter_t)
298 call tnsr3d_el(s22_svv, lx, s22, lx, filter_r, filter_s, filter_t)
299 call tnsr3d_el(s33_svv, lx, s33, lx, filter_r, filter_s, filter_t)
300 call tnsr3d_el(s12_svv, lx, s12, lx, filter_r, filter_s, filter_t)
301 call tnsr3d_el(s13_svv, lx, s13, lx, filter_r, filter_s, filter_t)
302 call tnsr3d_el(s23_svv, lx, s23, lx, filter_r, filter_s, filter_t)
303
304 do i = 1, lx*lx*lx
305 ! high pass filter from the LPF result
306 s11_svv(i,1,1) = s11(i,1,1) - s11_svv(i,1,1)
307 s22_svv(i,1,1) = s22(i,1,1) - s22_svv(i,1,1)
308 s33_svv(i,1,1) = s33(i,1,1) - s33_svv(i,1,1)
309 s12_svv(i,1,1) = s12(i,1,1) - s12_svv(i,1,1)
310 s13_svv(i,1,1) = s13(i,1,1) - s13_svv(i,1,1)
311 s23_svv(i,1,1) = s23(i,1,1) - s23_svv(i,1,1)
312
313 ! multiply the viscosity
314 s11_h = (svv_h1(i,1,1,e) * s11_svv(i,1,1) + &
315 h1(i,1,1,e) * s11(i,1,1)) * weights3(i,1,1)
316 s22_h = (svv_h1(i,1,1,e) * s22_svv(i,1,1) + &
317 h1(i,1,1,e) * s22(i,1,1)) * weights3(i,1,1)
318 s33_h = (svv_h1(i,1,1,e) * s33_svv(i,1,1) + &
319 h1(i,1,1,e) * s33(i,1,1)) * weights3(i,1,1)
320 s12_h = (svv_h1(i,1,1,e) * s12_svv(i,1,1) + &
321 h1(i,1,1,e) * s12(i,1,1)) * weights3(i,1,1)
322 s13_h = (svv_h1(i,1,1,e) * s13_svv(i,1,1) + &
323 h1(i,1,1,e) * s13(i,1,1)) * weights3(i,1,1)
324 s23_h = (svv_h1(i,1,1,e) * s23_svv(i,1,1) + &
325 h1(i,1,1,e) * s23(i,1,1)) * weights3(i,1,1)
326 ! utilize wur, wus, wut as work arrays again
327 wur(i,1,1) = drdx(i,1,1,e) * s11_h &
328 + drdy(i,1,1,e) * s12_h &
329 + drdz(i,1,1,e) * s13_h
330 wus(i,1,1) = dsdx(i,1,1,e) * s11_h &
331 + dsdy(i,1,1,e) * s12_h &
332 + dsdz(i,1,1,e) * s13_h
333 wut(i,1,1) = dtdx(i,1,1,e) * s11_h &
334 + dtdy(i,1,1,e) * s12_h &
335 + dtdz(i,1,1,e) * s13_h
336 wvr(i,1,1) = drdx(i,1,1,e) * s12_h &
337 + drdy(i,1,1,e) * s22_h &
338 + drdz(i,1,1,e) * s23_h
339 wvs(i,1,1) = dsdx(i,1,1,e) * s12_h &
340 + dsdy(i,1,1,e) * s22_h &
341 + dsdz(i,1,1,e) * s23_h
342 wvt(i,1,1) = dtdx(i,1,1,e) * s12_h &
343 + dtdy(i,1,1,e) * s22_h &
344 + dtdz(i,1,1,e) * s23_h
345 wwr(i,1,1) = drdx(i,1,1,e) * s13_h &
346 + drdy(i,1,1,e) * s23_h &
347 + drdz(i,1,1,e) * s33_h
348 wws(i,1,1) = dsdx(i,1,1,e) * s13_h &
349 + dsdy(i,1,1,e) * s23_h &
350 + dsdz(i,1,1,e) * s33_h
351 wwt(i,1,1) = dtdx(i,1,1,e) * s13_h &
352 + dtdy(i,1,1,e) * s23_h &
353 + dtdz(i,1,1,e) * s33_h
354 end do
355
356 do j = 1, lx*lx
357 do i = 1, lx
358 t1 = 0.0_rp
359 t2 = 0.0_rp
360 t3 = 0.0_rp
361 do k = 1, lx
362 t1 = t1 + dxt(i,k) * wur(k,j,1)
363 t2 = t2 + dxt(i,k) * wvr(k,j,1)
364 t3 = t3 + dxt(i,k) * wwr(k,j,1)
365 end do
366 au(i,j,1,e) = t1
367 av(i,j,1,e) = t2
368 aw(i,j,1,e) = t3
369 end do
370 end do
371
372 do k = 1, lx
373 do j = 1, lx
374 do i = 1, lx
375 t1 = 0.0_rp
376 t2 = 0.0_rp
377 t3 = 0.0_rp
378 do l = 1, lx
379 t1 = t1 + dyt(j,l) * wus(i,l,k)
380 t2 = t2 + dyt(j,l) * wvs(i,l,k)
381 t3 = t3 + dyt(j,l) * wws(i,l,k)
382 end do
383 au(i,j,k,e) = au(i,j,k,e) + t1
384 av(i,j,k,e) = av(i,j,k,e) + t2
385 aw(i,j,k,e) = aw(i,j,k,e) + t3
386 end do
387 end do
388 end do
389
390 do k = 1, lx
391 do i = 1, lx*lx
392 t1 = 0.0_rp
393 t2 = 0.0_rp
394 t3 = 0.0_rp
395 do l = 1, lx
396 t1 = t1 + dzt(k,l) * wut(i,1,l)
397 t2 = t2 + dzt(k,l) * wvt(i,1,l)
398 t3 = t3 + dzt(k,l) * wwt(i,1,l)
399 end do
400 au(i,1,k,e) = au(i,1,k,e) + t1
401 av(i,1,k,e) = av(i,1,k,e) + t2
402 aw(i,1,k,e) = aw(i,1,k,e) + t3
403 end do
404 end do
405
406 end do
407 end subroutine ax_helm_svv_one_sided_full_lx
408
409end module ax_helm_svv_one_sided_full_cpu
ax_helm_svv_full
Defines the full-stress SVV Helmholtz operator.
Definition ax_helm_svv_full.f90:34
ax_helm_svv_one_sided_full_cpu
CPU implementation of the full-stress one-sided SVV Helmholtz operator.
Definition ax_helm_svv_one_sided_full_cpu.f90:34
ax_helm_svv_one_sided_full_cpu::ax_helm_svv_one_sided_full_compute_vector
subroutine ax_helm_svv_one_sided_full_compute_vector(this, au, av, aw, u, v, w, coef, msh, xh)
Compute inside a Krylov method, taking 3 components of a vector field in a coupled manner.
Definition ax_helm_svv_one_sided_full_cpu.f90:69
ax_helm_svv_one_sided_full_cpu::ax_helm_svv_one_sided_full_lx
subroutine ax_helm_svv_one_sided_full_lx(au, av, aw, u, v, w, dx, dy, dz, dxt, dyt, dzt, h1, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, jacinv, weights3, svv_h1, svv_q, svv_qt, svv_direction, ident, n, lx)
Generic CPU kernel for the full-stress one-sided SVV product.
Definition ax_helm_svv_one_sided_full_cpu.f90:133
coefs
Coefficients.
Definition coef.f90:34
math
Definition math.f90:60
math::addcol4
subroutine, public addcol4(a, b, c, d, n)
Returns .
Definition math.f90:1183
mesh
Defines a mesh.
Definition mesh.f90:34
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:14
space
Defines a function space.
Definition space.f90:34
spectral_vanishing_viscosity
Data and filter construction for spectral vanishing viscosity.
Definition spectral_vanishing_viscosity.f90:33
tensor
Tensor operations.
Definition tensor.f90:61
tensor::tnsr3d_el
subroutine, public tnsr3d_el(v, nv, u, nu, a, bt, ct)
Tensor product performed on a single element.
Definition tensor.f90:174
tensor::tnsr3d
subroutine, public tnsr3d(v, nv, u, nu, a, bt, ct, nelv)
Tensor product performed on nelv elements.
Definition tensor.f90:234
ax_helm_svv_full::ax_helm_svv_full_t
Matrix-vector product for a Helmholtz problem.
Definition ax_helm_svv_full.f90:47
ax_helm_svv_one_sided_full_cpu::ax_helm_svv_one_sided_full_cpu_t
CPU matrix-vector product for a Helmholtz problem.
Definition ax_helm_svv_one_sided_full_cpu.f90:47
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:93
mesh::mesh_t
Definition mesh.f90:73
space::space_t
The function space for the SEM solution fields.
Definition space.f90:64
spectral_vanishing_viscosity::svv_t
Spectral vanishing viscosity configuration and coefficients.
Definition spectral_vanishing_viscosity.f90:61