Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
wale_cpu.f90
Go to the documentation of this file.
1! Copyright (c) 2023-2025, 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!
35 use num_types, only : rp
36 use field_list, only : field_list_t
39 use field, only : field_t
40 use operators, only : dudxyz, strain_rate
41 use coefs, only : coef_t
42 use gs_ops, only : gs_op_add
43 use math, only : col2, neko_eps
44 implicit none
45 private
46
47 public :: wale_compute_cpu
48
49contains
50
59 subroutine wale_compute_cpu(if_ext, t, tstep, coef, nut, delta, c_w)
60 logical, intent(in) :: if_ext
61 real(kind=rp), intent(in) :: t
62 integer, intent(in) :: tstep
63 type(coef_t), intent(in) :: coef
64 type(field_t), intent(inout) :: nut
65 type(field_t), intent(in) :: delta
66 real(kind=rp), intent(in) :: c_w
67 type(field_t), pointer :: u, v, w
68 ! the velocity gradient tensor
69 type(field_t), pointer :: g11, g12, g13, g21, g22, g23, g31, g32, g33
70 ! strain rate tensor
71 type(field_t), pointer :: s11, s22, s33, s12, s13, s23
72
73 real(kind=rp) :: gsqr_11, gsqr_12, gsqr_13, gsqr_21, gsqr_22, gsqr_23, gsqr_31, gsqr_32, gsqr_33
74 real(kind=rp) :: sd11, sd22, sd33, sd12, sd13, sd23
75 real(kind=rp) :: sdij_sdij
76 real(kind=rp) :: sij_sij
77 real(kind=rp) :: op_wale
78 integer :: temp_indices(15)
79 integer :: e, i
80
81
82 if (if_ext .eqv. .true.) then
83 u => neko_field_registry%get_field_by_name("u_e")
84 v => neko_field_registry%get_field_by_name("v_e")
85 w => neko_field_registry%get_field_by_name("w_e")
86 else
87 u => neko_field_registry%get_field_by_name("u")
88 v => neko_field_registry%get_field_by_name("v")
89 w => neko_field_registry%get_field_by_name("w")
90 end if
91
92 call neko_scratch_registry%request_field(g11, temp_indices(1))
93 call neko_scratch_registry%request_field(g12, temp_indices(2))
94 call neko_scratch_registry%request_field(g13, temp_indices(3))
95 call neko_scratch_registry%request_field(g21, temp_indices(4))
96 call neko_scratch_registry%request_field(g22, temp_indices(5))
97 call neko_scratch_registry%request_field(g23, temp_indices(6))
98 call neko_scratch_registry%request_field(g31, temp_indices(7))
99 call neko_scratch_registry%request_field(g32, temp_indices(8))
100 call neko_scratch_registry%request_field(g33, temp_indices(9))
101
102 call neko_scratch_registry%request_field(s11, temp_indices(10))
103 call neko_scratch_registry%request_field(s22, temp_indices(11))
104 call neko_scratch_registry%request_field(s33, temp_indices(12))
105 call neko_scratch_registry%request_field(s12, temp_indices(13))
106 call neko_scratch_registry%request_field(s13, temp_indices(14))
107 call neko_scratch_registry%request_field(s23, temp_indices(15))
108
109
110 ! Compute the velocity gradient tensor g
111 call dudxyz(g11%x, u%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
112 call dudxyz(g12%x, u%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
113 call dudxyz(g13%x, u%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
114
115 call dudxyz(g21%x, v%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
116 call dudxyz(g22%x, v%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
117 call dudxyz(g23%x, v%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
118
119 call dudxyz(g31%x, w%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
120 call dudxyz(g32%x, w%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
121 call dudxyz(g33%x, w%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
122
123 call coef%gs_h%op(g11, gs_op_add)
124 call coef%gs_h%op(g12, gs_op_add)
125 call coef%gs_h%op(g13, gs_op_add)
126 call coef%gs_h%op(g21, gs_op_add)
127 call coef%gs_h%op(g22, gs_op_add)
128 call coef%gs_h%op(g23, gs_op_add)
129 call coef%gs_h%op(g31, gs_op_add)
130 call coef%gs_h%op(g32, gs_op_add)
131 call coef%gs_h%op(g33, gs_op_add)
132
133 ! Compute components of the strain rate tensor
134 call strain_rate(s11%x, s22%x, s33%x, s12%x, s13%x, s23%x, u, v, w, coef)
135
136 call coef%gs_h%op(s11, gs_op_add)
137 call coef%gs_h%op(s22, gs_op_add)
138 call coef%gs_h%op(s33, gs_op_add)
139 call coef%gs_h%op(s12, gs_op_add)
140 call coef%gs_h%op(s13, gs_op_add)
141 call coef%gs_h%op(s23, gs_op_add)
142
143
144 do concurrent(e = 1:coef%msh%nelv)
145 do concurrent(i = 1:coef%Xh%lxyz)
146 ! gij^2 = g_ik * g_kj
147 gsqr_11 = g11%x(i,1,1,e)*g11%x(i,1,1,e) + &
148 g12%x(i,1,1,e)*g21%x(i,1,1,e) + &
149 g13%x(i,1,1,e)*g31%x(i,1,1,e)
150 gsqr_12 = g11%x(i,1,1,e)*g12%x(i,1,1,e) + &
151 g12%x(i,1,1,e)*g22%x(i,1,1,e) + &
152 g13%x(i,1,1,e)*g32%x(i,1,1,e)
153 gsqr_13 = g11%x(i,1,1,e)*g13%x(i,1,1,e) + &
154 g12%x(i,1,1,e)*g23%x(i,1,1,e) + &
155 g13%x(i,1,1,e)*g33%x(i,1,1,e)
156 gsqr_21 = g21%x(i,1,1,e)*g11%x(i,1,1,e) + &
157 g22%x(i,1,1,e)*g21%x(i,1,1,e) + &
158 g23%x(i,1,1,e)*g31%x(i,1,1,e)
159 gsqr_22 = g21%x(i,1,1,e)*g12%x(i,1,1,e) + &
160 g22%x(i,1,1,e)*g22%x(i,1,1,e) + &
161 g23%x(i,1,1,e)*g32%x(i,1,1,e)
162 gsqr_23 = g21%x(i,1,1,e)*g13%x(i,1,1,e) + &
163 g22%x(i,1,1,e)*g23%x(i,1,1,e) + &
164 g23%x(i,1,1,e)*g33%x(i,1,1,e)
165 gsqr_31 = g31%x(i,1,1,e)*g11%x(i,1,1,e) + &
166 g32%x(i,1,1,e)*g21%x(i,1,1,e) + &
167 g33%x(i,1,1,e)*g31%x(i,1,1,e)
168 gsqr_32 = g31%x(i,1,1,e)*g12%x(i,1,1,e) + &
169 g32%x(i,1,1,e)*g22%x(i,1,1,e) + &
170 g33%x(i,1,1,e)*g32%x(i,1,1,e)
171 gsqr_33 = g31%x(i,1,1,e)*g13%x(i,1,1,e) + &
172 g32%x(i,1,1,e)*g23%x(i,1,1,e) + &
173 g33%x(i,1,1,e)*g33%x(i,1,1,e)
174
175 ! sdij components
176 sd11 = gsqr_11 - ( (gsqr_11 + gsqr_22 + gsqr_33) / 3.0_rp)
177 sd22 = gsqr_22 - ( (gsqr_11 + gsqr_22 + gsqr_33) / 3.0_rp)
178 sd33 = gsqr_33 - ( (gsqr_11 + gsqr_22 + gsqr_33) / 3.0_rp)
179 sd12 = 0.5_rp * (gsqr_12 + gsqr_21)
180 sd13 = 0.5_rp * (gsqr_13 + gsqr_31)
181 sd23 = 0.5_rp * (gsqr_23 + gsqr_32)
182
183 ! Sdij*Sdij
184 sdij_sdij = sd11*sd11 + sd22*sd22 + sd33*sd33 + &
185 2.0_rp * (sd12*sd12 + sd13*sd13 + sd23*sd23)
186 ! Sij*Sij
187 sij_sij = s11%x(i,1,1,e)*s11%x(i,1,1,e) + s22%x(i,1,1,e)*s22%x(i,1,1,e) + &
188 s33%x(i,1,1,e)*s33%x(i,1,1,e) + 2.0_rp * (s12%x(i,1,1,e)*s12%x(i,1,1,e) + &
189 s13%x(i,1,1,e)*s13%x(i,1,1,e) + s23%x(i,1,1,e)*s23%x(i,1,1,e))
190
191 ! Wale operator
192 op_wale = sdij_sdij**(3.0_rp / 2.0_rp) / &
193 max((sij_sij**(5.0_rp / 2.0_rp) + sdij_sdij**(5.0_rp / 4.0_rp)), neko_eps)
194
195 ! turbulent viscosity
196 nut%x(i,1,1,e) = c_w**2 * delta%x(i,1,1,e)**2 * op_wale * coef%mult(i,1,1,e)
197 end do
198 end do
199
200 call neko_scratch_registry%relinquish_field(temp_indices)
201 end subroutine wale_compute_cpu
202
203end module wale_cpu
204
Coefficients.
Definition coef.f90:34
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Defines a field.
Definition field.f90:34
Defines Gather-scatter operations.
Definition gs_ops.f90:34
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
Definition math.f90:60
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:860
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition math.f90:69
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
subroutine, public strain_rate(s11, s22, s33, s12, s13, s23, u, v, w, coef)
Compute the strain rate tensor, i.e 0.5 * du_i/dx_j + du_j/dx_i.
subroutine, public dudxyz(du, u, dr, ds, dt, coef)
Compute derivative of a scalar field along a single direction.
Definition operators.f90:80
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
type(scratch_registry_t), target, public neko_scratch_registry
Global scratch registry.
Implements the CPU kernel for the wale_t type.
Definition wale_cpu.f90:34
subroutine, public wale_compute_cpu(if_ext, t, tstep, coef, nut, delta, c_w)
Compute eddy viscosity on the CPU.
Definition wale_cpu.f90:60
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
field_list_t, To be able to group fields together
#define max(a, b)
Definition tensor.cu:40