Loading [MathJax]/extensions/tex2jax.js
Neko 0.9.99
A portable framework for high-order spectral element flow simulations
All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Pages
euler_res_cpu.f90
Go to the documentation of this file.
1! Copyright (c) 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!
34!? It handles the time advancement of primitive variables using Runge-Kutta methods
35!? and evaluates the right-hand side terms of the Euler equations including artificial viscosity.
37 use euler_residual, only : euler_rhs_t
38 use field, only : field_t
39 use ax_product, only : ax_t
40 use coefs, only : coef_t
41 use gather_scatter, only : gs_t
42 use num_types, only : rp
43 use operators, only: div
44 use math, only: subcol3, copy, sub2, add2, add3, &
46 use gs_ops, only : gs_op_add
49 use field_list, only : field_list_t
50 implicit none
51 private
52
53 type, public, extends(euler_rhs_t) :: euler_res_cpu_t
54 contains
55 procedure, nopass :: step => advance_primitive_variables_cpu
56 procedure, nopass :: evaluate_rhs_cpu
57 end type euler_res_cpu_t
58
59contains
74 subroutine advance_primitive_variables_cpu(rho_field, m_x, m_y, m_z, &
75 E, p, u, v, w, Ax, &
76 coef, gs, h, c_avisc_low, rk_scheme, dt)
77 type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E
78 type(field_t), intent(in) :: p, u, v, w, h
79 class(ax_t), intent(inout) :: Ax
80 type(coef_t), intent(inout) :: coef
81 type(gs_t), intent(inout) :: gs
82 real(kind=rp) :: c_avisc_low
83 class(runge_kutta_time_scheme_t), intent(in) :: rk_scheme
84 real(kind=rp), intent(in) :: dt
85 integer :: n, s, i, j, k
86 real(kind=rp) :: t, c
87 type(field_t), pointer :: k_rho_1, k_rho_2, k_rho_3, k_rho_4, &
88 k_m_x_1, k_m_x_2, k_m_x_3, k_m_x_4, &
89 k_m_y_1, k_m_y_2, k_m_y_3, k_m_y_4, &
90 k_m_z_1, k_m_z_2, k_m_z_3, k_m_z_4, &
91 k_e_1, k_e_2, k_e_3, k_e_4, &
92 temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_e
93 integer :: temp_indices(25)
94 type(field_list_t) :: k_rho, k_m_x, k_m_y, k_m_z, k_E
95
96 n = p%dof%size()
97 s = rk_scheme%order
98 call neko_scratch_registry%request_field(k_rho_1, temp_indices(1))
99 call neko_scratch_registry%request_field(k_rho_2, temp_indices(2))
100 call neko_scratch_registry%request_field(k_rho_3, temp_indices(3))
101 call neko_scratch_registry%request_field(k_rho_4, temp_indices(4))
102 call neko_scratch_registry%request_field(k_m_x_1, temp_indices(5))
103 call neko_scratch_registry%request_field(k_m_x_2, temp_indices(6))
104 call neko_scratch_registry%request_field(k_m_x_3, temp_indices(7))
105 call neko_scratch_registry%request_field(k_m_x_4, temp_indices(8))
106 call neko_scratch_registry%request_field(k_m_y_1, temp_indices(9))
107 call neko_scratch_registry%request_field(k_m_y_2, temp_indices(10))
108 call neko_scratch_registry%request_field(k_m_y_3, temp_indices(11))
109 call neko_scratch_registry%request_field(k_m_y_4, temp_indices(12))
110 call neko_scratch_registry%request_field(k_m_z_1, temp_indices(13))
111 call neko_scratch_registry%request_field(k_m_z_2, temp_indices(14))
112 call neko_scratch_registry%request_field(k_m_z_3, temp_indices(15))
113 call neko_scratch_registry%request_field(k_m_z_4, temp_indices(16))
114 call neko_scratch_registry%request_field(k_e_1, temp_indices(17))
115 call neko_scratch_registry%request_field(k_e_2, temp_indices(18))
116 call neko_scratch_registry%request_field(k_e_3, temp_indices(19))
117 call neko_scratch_registry%request_field(k_e_4, temp_indices(20))
118 call neko_scratch_registry%request_field(temp_rho, temp_indices(21))
119 call neko_scratch_registry%request_field(temp_m_x, temp_indices(22))
120 call neko_scratch_registry%request_field(temp_m_y, temp_indices(23))
121 call neko_scratch_registry%request_field(temp_m_z, temp_indices(24))
122 call neko_scratch_registry%request_field(temp_e, temp_indices(25))
123
124 ! Initialize Runge-Kutta stage variables for each conserved quantity
125 call k_rho%init(4)
126 call k_rho%assign(1, k_rho_1)
127 call k_rho%assign(2, k_rho_2)
128 call k_rho%assign(3, k_rho_3)
129 call k_rho%assign(4, k_rho_4)
130 call k_m_x%init(4)
131 call k_m_x%assign(1, k_m_x_1)
132 call k_m_x%assign(2, k_m_x_2)
133 call k_m_x%assign(3, k_m_x_3)
134 call k_m_x%assign(4, k_m_x_4)
135 call k_m_y%init(4)
136 call k_m_y%assign(1, k_m_y_1)
137 call k_m_y%assign(2, k_m_y_2)
138 call k_m_y%assign(3, k_m_y_3)
139 call k_m_y%assign(4, k_m_y_4)
140 call k_m_z%init(4)
141 call k_m_z%assign(1, k_m_z_1)
142 call k_m_z%assign(2, k_m_z_2)
143 call k_m_z%assign(3, k_m_z_3)
144 call k_m_z%assign(4, k_m_z_4)
145 call k_e%init(4)
146 call k_e%assign(1, k_e_1)
147 call k_e%assign(2, k_e_2)
148 call k_e%assign(3, k_e_3)
149 call k_e%assign(4, k_e_4)
150
151 ! Loop over Runge-Kutta stages
152 do i = 1, s
153 ! Copy current solution state to temporary arrays for this RK stage
154 call copy(temp_rho%x, rho_field%x, n)
155 call copy(temp_m_x%x, m_x%x, n)
156 call copy(temp_m_y%x, m_y%x, n)
157 call copy(temp_m_z%x, m_z%x, n)
158 call copy(temp_e%x, e%x, n)
159
160 ! Accumulate previous stage contributions using RK coefficients
161 do j = 1, i-1
162 do concurrent(k = 1:n)
163 temp_rho%x(k,1,1,1) = temp_rho%x(k,1,1,1) &
164 + dt * rk_scheme%coeffs_A(i, j) * k_rho%items(j)%ptr%x(k,1,1,1)
165 temp_m_x%x(k,1,1,1) = temp_m_x%x(k,1,1,1) &
166 + dt * rk_scheme%coeffs_A(i, j) * k_m_x%items(j)%ptr%x(k,1,1,1)
167 temp_m_y%x(k,1,1,1) = temp_m_y%x(k,1,1,1) &
168 + dt * rk_scheme%coeffs_A(i, j) * k_m_y%items(j)%ptr%x(k,1,1,1)
169 temp_m_z%x(k,1,1,1) = temp_m_z%x(k,1,1,1) &
170 + dt * rk_scheme%coeffs_A(i, j) * k_m_z%items(j)%ptr%x(k,1,1,1)
171 temp_e%x(k,1,1,1) = temp_e%x(k,1,1,1) &
172 + dt * rk_scheme%coeffs_A(i, j) * k_e%items(j)%ptr%x(k,1,1,1)
173 end do
174 end do
175
176 ! Evaluate RHS terms for current stage using intermediate solution values
177 call evaluate_rhs_cpu(k_rho%items(i)%ptr, k_m_x%items(i)%ptr, &
178 k_m_y%items(i)%ptr, k_m_z%items(i)%ptr, &
179 k_e%items(i)%ptr, &
180 temp_rho, temp_m_x, temp_m_y, temp_m_z, temp_e, &
181 p, u, v, w, ax, &
182 coef, gs, h, c_avisc_low)
183 end do
184
185 ! Update the solution
186 do i = 1, s
187 do concurrent(k = 1:n)
188 rho_field%x(k,1,1,1) = rho_field%x(k,1,1,1) &
189 + dt * rk_scheme%coeffs_b(i) * k_rho%items(i)%ptr%x(k,1,1,1)
190 m_x%x(k,1,1,1) = m_x%x(k,1,1,1) &
191 + dt * rk_scheme%coeffs_b(i) * k_m_x%items(i)%ptr%x(k,1,1,1)
192 m_y%x(k,1,1,1) = m_y%x(k,1,1,1) &
193 + dt * rk_scheme%coeffs_b(i) * k_m_y%items(i)%ptr%x(k,1,1,1)
194 m_z%x(k,1,1,1) = m_z%x(k,1,1,1) &
195 + dt * rk_scheme%coeffs_b(i) * k_m_z%items(i)%ptr%x(k,1,1,1)
196 e%x(k,1,1,1) = e%x(k,1,1,1) &
197 + dt * rk_scheme%coeffs_b(i) * k_e%items(i)%ptr%x(k,1,1,1)
198 end do
199 end do
200
201 call neko_scratch_registry%relinquish_field(temp_indices)
202
204
226 subroutine evaluate_rhs_cpu(rhs_rho_field, rhs_m_x, rhs_m_y, rhs_m_z, rhs_E, &
227 rho_field, m_x, m_y, m_z, E, p, u, v, w, Ax, &
228 coef, gs, h, c_avisc_low)
229 type(field_t), intent(inout) :: rhs_rho_field, &
230 rhs_m_x, rhs_m_y, rhs_m_z, rhs_E
231 type(field_t), intent(inout) :: rho_field, m_x, m_y, m_z, E
232 type(field_t), intent(in) :: p, u, v, w, h
233 class(ax_t), intent(inout) :: Ax
234 type(coef_t), intent(inout) :: coef
235 type(gs_t), intent(inout) :: gs
236 real(kind=rp) :: c_avisc_low
237 integer :: i, n
238 type(field_t), pointer :: temp, f_x, f_y, f_z, &
239 visc_rho, visc_m_x, visc_m_y, visc_m_z, visc_E
240 integer :: temp_indices(9)
241
242 n = coef%dof%size()
243 call neko_scratch_registry%request_field(temp, temp_indices(1))
244 call neko_scratch_registry%request_field(f_x, temp_indices(2))
245 call neko_scratch_registry%request_field(f_y, temp_indices(3))
246 call neko_scratch_registry%request_field(f_z, temp_indices(4))
247
249 ! Compute density flux divergence
250 call div(rhs_rho_field%x, m_x%x, m_y%x, m_z%x, coef)
251
253 ! Compute momentum flux divergences
254 ! m_x
255 do concurrent(i = 1:n)
256 f_x%x(i,1,1,1) = m_x%x(i,1,1,1) * m_x%x(i,1,1,1) / rho_field%x(i, 1, 1, 1) &
257 + p%x(i,1,1,1)
258 f_y%x(i,1,1,1) = m_x%x(i,1,1,1) * m_y%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
259 f_z%x(i,1,1,1) = m_x%x(i,1,1,1) * m_z%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
260 end do
261 call div(rhs_m_x%x, f_x%x, f_y%x, f_z%x, coef)
262 ! m_y
263 do concurrent(i = 1:n)
264 f_x%x(i,1,1,1) = m_y%x(i,1,1,1) * m_x%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
265 f_y%x(i,1,1,1) = m_y%x(i,1,1,1) * m_y%x(i,1,1,1) / rho_field%x(i, 1, 1, 1) &
266 + p%x(i,1,1,1)
267 f_z%x(i,1,1,1) = m_y%x(i,1,1,1) * m_z%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
268 end do
269 call div(rhs_m_y%x, f_x%x, f_y%x, f_z%x, coef)
270 ! m_z
271 do concurrent(i = 1:n)
272 f_x%x(i,1,1,1) = m_z%x(i,1,1,1) * m_x%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
273 f_y%x(i,1,1,1) = m_z%x(i,1,1,1) * m_y%x(i,1,1,1) / rho_field%x(i, 1, 1, 1)
274 f_z%x(i,1,1,1) = m_z%x(i,1,1,1) * m_z%x(i,1,1,1) / rho_field%x(i, 1, 1, 1) &
275 + p%x(i,1,1,1)
276 end do
277 call div(rhs_m_z%x, f_x%x, f_y%x, f_z%x, coef)
278
280 ! Compute energy flux divergence
281 do concurrent(i = 1:n)
282 f_x%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) * u%x(i,1,1,1)
283 f_y%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) * v%x(i,1,1,1)
284 f_z%x(i,1,1,1) = (e%x(i,1,1,1) + p%x(i,1,1,1)) * w%x(i,1,1,1)
285 end do
286 call div(rhs_e%x, f_x%x, f_y%x, f_z%x, coef)
287
288 ! gs
289 call gs%op(rhs_rho_field, gs_op_add)
290 call gs%op(rhs_m_x, gs_op_add)
291 call gs%op(rhs_m_y, gs_op_add)
292 call gs%op(rhs_m_z, gs_op_add)
293 call gs%op(rhs_e, gs_op_add)
294 do concurrent(i = 1:rhs_e%dof%size())
295 rhs_rho_field%x(i,1,1,1) = rhs_rho_field%x(i,1,1,1) * coef%mult(i,1,1,1)
296 rhs_m_x%x(i,1,1,1) = rhs_m_x%x(i,1,1,1) * coef%mult(i,1,1,1)
297 rhs_m_y%x(i,1,1,1) = rhs_m_y%x(i,1,1,1) * coef%mult(i,1,1,1)
298 rhs_m_z%x(i,1,1,1) = rhs_m_z%x(i,1,1,1) * coef%mult(i,1,1,1)
299 rhs_e%x(i,1,1,1) = rhs_e%x(i,1,1,1) * coef%mult(i,1,1,1)
300 end do
301
302 call neko_scratch_registry%request_field(visc_rho, temp_indices(5))
303 call neko_scratch_registry%request_field(visc_m_x, temp_indices(6))
304 call neko_scratch_registry%request_field(visc_m_y, temp_indices(7))
305 call neko_scratch_registry%request_field(visc_m_z, temp_indices(8))
306 call neko_scratch_registry%request_field(visc_e, temp_indices(9))
307
308 ! Calculate artificial diffusion
309 call ax%compute(visc_rho%x, rho_field%x, coef, p%msh, p%Xh)
310 call ax%compute(visc_m_x%x, m_x%x, coef, p%msh, p%Xh)
311 call ax%compute(visc_m_y%x, m_y%x, coef, p%msh, p%Xh)
312 call ax%compute(visc_m_z%x, m_z%x, coef, p%msh, p%Xh)
313 call ax%compute(visc_e%x, e%x, coef, p%msh, p%Xh)
314
315 ! gs
316 call gs%op(visc_rho, gs_op_add)
317 call gs%op(visc_m_x, gs_op_add)
318 call gs%op(visc_m_y, gs_op_add)
319 call gs%op(visc_m_z, gs_op_add)
320 call gs%op(visc_e, gs_op_add)
321
322 ! Move div to the rhs and apply artificial viscosity
323 ! i.e., calculate -div(grad(f)) + div(visc*grad(u))
324 do concurrent(i = 1:n)
325 rhs_rho_field%x(i,1,1,1) = -rhs_rho_field%x(i,1,1,1) &
326 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_rho%x(i,1,1,1)
327 rhs_m_x%x(i,1,1,1) = -rhs_m_x%x(i,1,1,1) &
328 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_x%x(i,1,1,1)
329 rhs_m_y%x(i,1,1,1) = -rhs_m_y%x(i,1,1,1) &
330 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_y%x(i,1,1,1)
331 rhs_m_z%x(i,1,1,1) = -rhs_m_z%x(i,1,1,1) &
332 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_m_z%x(i,1,1,1)
333 rhs_e%x(i,1,1,1) = -rhs_e%x(i,1,1,1) &
334 - c_avisc_low * h%x(i,1,1,1) * coef%Binv(i,1,1,1) * visc_e%x(i,1,1,1)
335 end do
336
337 call neko_scratch_registry%relinquish_field(temp_indices)
338 end subroutine evaluate_rhs_cpu
339
340end module euler_res_cpu
Defines a Matrix-vector product.
Definition ax.f90:34
Coefficients.
Definition coef.f90:34
This module implements CPU-based residual calculations for the Euler equations.
subroutine evaluate_rhs_cpu(rhs_rho_field, rhs_m_x, rhs_m_y, rhs_m_z, rhs_e, rho_field, m_x, m_y, m_z, e, p, u, v, w, ax, coef, gs, h, c_avisc_low)
Evaluates the right-hand side of the Euler equations Computes fluxes, divergence terms,...
subroutine advance_primitive_variables_cpu(rho_field, m_x, m_y, m_z, e, p, u, v, w, ax, coef, gs, h, c_avisc_low, rk_scheme, dt)
Advances the primitive variables (density, momentum, energy) in time using a Runge-Kutta scheme.
Defines a field.
Definition field.f90:34
Gather-scatter.
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 cmult(a, c, n)
Multiplication by constant c .
Definition math.f90:310
subroutine, public subcol3(a, b, c, n)
Returns .
Definition math.f90:755
subroutine, public add3(a, b, c, n)
Vector addition .
Definition math.f90:599
subroutine, public add2(a, b, n)
Vector addition .
Definition math.f90:586
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition math.f90:347
subroutine, public invcol3(a, b, c, n)
Invert a vector .
Definition math.f90:486
subroutine, public addcol3(a, b, c, n)
Returns .
Definition math.f90:800
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:728
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:238
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:741
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:628
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
subroutine, public div(res, ux, uy, uz, coef)
Compute the divergence of a vector field.
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.
Base type for a matrix-vector product providing .
Definition ax.f90:43
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Abstract type to compute rhs.
Definition euler_res.f90:47
field_list_t, To be able to group fields together