Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
fluid_scheme_compressible.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 use field, only : field_t
37
40 use json_module, only : json_file
41 use num_types, only : rp
42 use mesh, only : mesh_t
44 use space, only : gll
45 use user_intf, only : user_t
47 use mpi_f08
48 use operators, only : cfl_compressible
54 use time_state, only : time_state_t
55 use logger, only : neko_log, log_size
56 use math, only : glsum
57 use num_types, only : i8
58 implicit none
59 private
60
62 type, public, abstract, extends(fluid_scheme_base_t) :: fluid_scheme_compressible_t
64 type(field_t), pointer :: m_x => null()
65 type(field_t), pointer :: m_y => null()
66 type(field_t), pointer :: m_z => null()
67 type(field_t), pointer :: e => null()
68 type(field_t), pointer :: max_wave_speed => null()
69 type(field_t), pointer :: s => null()
70
71 real(kind=rp) :: gamma
72
74 integer(kind=i8) :: glb_n_points
76 integer(kind=i8) :: glb_unique_points
77
78 type(scratch_registry_t) :: scratch
79
80 contains
82 procedure, pass(this) :: scheme_init => fluid_scheme_compressible_init
84 procedure, pass(this) :: scheme_free => fluid_scheme_compressible_free
85
87 procedure, pass(this) :: validate => fluid_scheme_compressible_validate
89 procedure, pass(this) :: compute_cfl &
92 procedure, pass(this) :: update_material_properties => &
95 procedure, pass(this) :: compute_entropy => &
98 procedure, pass(this) :: compute_max_wave_speed => &
101 procedure, pass(this) :: log_solver_info => &
103
105
106contains
114 subroutine fluid_scheme_compressible_init(this, msh, lx, params, scheme, user)
115 class(fluid_scheme_compressible_t), target, intent(inout) :: this
116 type(mesh_t), target, intent(inout) :: msh
117 integer, intent(in) :: lx
118 character(len=*), intent(in) :: scheme
119 type(json_file), target, intent(inout) :: params
120 type(user_t), target, intent(in) :: user
121
122 !
123 ! SEM simulation fundamentals
124 !
125
126 this%msh => msh
127
128 if (msh%gdim .eq. 2) then
129 call this%Xh%init(gll, lx, lx)
130 else
131 call this%Xh%init(gll, lx, lx, lx)
132 end if
133
134 call this%dm_Xh%init(msh, this%Xh)
135
136 call this%gs_Xh%init(this%dm_Xh)
137
138 call this%c_Xh%init(this%gs_Xh)
139
140 ! Local scratch registry
141 call this%scratch%init(this%dm_Xh, 10, 2)
142
143 ! Case parameters
144 this%params => params
145
146 ! Assign a name
147 call json_get_or_default(params, 'case.fluid.name', this%name, "fluid")
148
149 ! Fill mu and rho field with the physical value
150 call neko_field_registry%add_field(this%dm_Xh, this%name // "_mu")
151 call neko_field_registry%add_field(this%dm_Xh, this%name // "_rho")
152 this%mu => neko_field_registry%get_field(this%name // "_mu")
153 this%rho => neko_field_registry%get_field(this%name // "_rho")
154 call field_cfill(this%mu, 0.0_rp, this%mu%size())
155
156 ! Assign momentum fields
157 call neko_field_registry%add_field(this%dm_Xh, "m_x")
158 call neko_field_registry%add_field(this%dm_Xh, "m_y")
159 call neko_field_registry%add_field(this%dm_Xh, "m_z")
160 this%m_x => neko_field_registry%get_field("m_x")
161 this%m_y => neko_field_registry%get_field("m_y")
162 this%m_z => neko_field_registry%get_field("m_z")
163 call this%m_x%init(this%dm_Xh, "m_x")
164 call this%m_y%init(this%dm_Xh, "m_y")
165 call this%m_z%init(this%dm_Xh, "m_z")
166
167 ! Assign energy field
168 call neko_field_registry%add_field(this%dm_Xh, "E")
169 this%E => neko_field_registry%get_field("E")
170 call this%E%init(this%dm_Xh, "E")
171
172 ! Assign maximum wave speed field
173 call neko_field_registry%add_field(this%dm_Xh, "max_wave_speed")
174 this%max_wave_speed => neko_field_registry%get_field("max_wave_speed")
175 call this%max_wave_speed%init(this%dm_Xh, "max_wave_speed")
176
177 ! Assign entropy field
178 call neko_field_registry%add_field(this%dm_Xh, "S")
179 this%S => neko_field_registry%get_field("S")
180 call this%S%init(this%dm_Xh, "S")
181
182 ! ! Assign velocity fields
183 call neko_field_registry%add_field(this%dm_Xh, "u")
184 call neko_field_registry%add_field(this%dm_Xh, "v")
185 call neko_field_registry%add_field(this%dm_Xh, "w")
186 this%u => neko_field_registry%get_field("u")
187 this%v => neko_field_registry%get_field("v")
188 this%w => neko_field_registry%get_field("w")
189 call this%u%init(this%dm_Xh, "u")
190 call this%v%init(this%dm_Xh, "v")
191 call this%w%init(this%dm_Xh, "w")
192 call neko_field_registry%add_field(this%dm_Xh, 'p')
193 this%p => neko_field_registry%get_field('p')
194 call this%p%init(this%dm_Xh, "p")
195
196 !
197 ! Setup right-hand side fields.
198 !
199 allocate(this%f_x)
200 allocate(this%f_y)
201 allocate(this%f_z)
202 call this%f_x%init(this%dm_Xh, fld_name = "fluid_rhs_x")
203 call this%f_y%init(this%dm_Xh, fld_name = "fluid_rhs_y")
204 call this%f_z%init(this%dm_Xh, fld_name = "fluid_rhs_z")
205
206 ! Compressible parameters
207 call json_get_or_default(params, 'case.fluid.gamma', this%gamma, 1.4_rp)
208
209 ! Calculate global points for logging
210 this%glb_n_points = int(this%msh%glb_nelv, i8)*int(this%Xh%lxyz, i8)
211 this%glb_unique_points = int(glsum(this%c_Xh%mult, this%dm_Xh%size()), i8)
212
213 !
214 ! Log solver information
215 !
216 call this%log_solver_info(params, scheme, lx)
217 end subroutine fluid_scheme_compressible_init
218
222 class(fluid_scheme_compressible_t), intent(inout) :: this
223 call this%dm_Xh%free()
224 call this%gs_Xh%free()
225 call this%c_Xh%free()
226 call this%Xh%free()
227
228 if (associated(this%m_x)) then
229 call this%m_x%free()
230 end if
231
232 if (associated(this%m_y)) then
233 call this%m_y%free()
234 end if
235
236 if (associated(this%m_z)) then
237 call this%m_z%free()
238 end if
239
240 if (associated(this%E)) then
241 call this%E%free()
242 end if
243
244 if (associated(this%max_wave_speed)) then
245 call this%max_wave_speed%free()
246 end if
247
248 if (associated(this%S)) then
249 call this%S%free()
250 end if
251
252 nullify(this%m_x)
253 nullify(this%m_y)
254 nullify(this%m_z)
255 nullify(this%E)
256 nullify(this%max_wave_speed)
257 nullify(this%S)
258
259 nullify(this%u)
260 nullify(this%v)
261 nullify(this%w)
262 nullify(this%p)
263 nullify(this%rho)
264 nullify(this%mu)
265
266 end subroutine fluid_scheme_compressible_free
267
271 class(fluid_scheme_compressible_t), target, intent(inout) :: this
272 integer :: n
273 type(field_t), pointer :: temp
274 integer :: temp_indices(1)
275
276 n = this%dm_Xh%size()
277 call this%scratch%request_field(temp, temp_indices(1))
278
280 call field_col3(this%m_x, this%u, this%rho)
281 call field_col3(this%m_y, this%v, this%rho)
282 call field_col3(this%m_z, this%w, this%rho)
283
286 call field_cmult2(this%E, this%p, 1.0_rp/(this%gamma - 1.0_rp), n)
287 call field_col3(temp, this%u, this%u, n)
288 call field_addcol3(temp, this%v, this%v, n)
289 call field_addcol3(temp, this%w, this%w, n)
290 call field_col2(temp, this%rho, n)
291 call field_cmult(temp, 0.5_rp, n)
292 call field_add2(this%E, temp, n)
293
294 call this%scratch%relinquish_field(temp_indices)
295
297 call this%compute_max_wave_speed()
298
300
305 function fluid_scheme_compressible_compute_cfl(this, dt) result(c)
306 class(fluid_scheme_compressible_t), intent(in) :: this
307 real(kind=rp), intent(in) :: dt
308 real(kind=rp) :: c
309 integer :: n
310
311 associate(u => this%u, v => this%v, w => this%w, p => this%p, &
312 rho => this%rho, xh => this%Xh, c_xh => this%c_Xh, &
313 msh => this%msh, gamma => this%gamma, &
314 max_wave_speed => this%max_wave_speed)
315
316 n = xh%lx * xh%ly * xh%lz * msh%nelv
317
318 ! Use the compressible CFL function with precomputed maximum wave speed
319 c = cfl_compressible(dt, max_wave_speed%x, xh, c_xh, msh%nelv, msh%gdim)
320 end associate
321
323
327 class(fluid_scheme_compressible_t), intent(inout) :: this
328 type(time_state_t), intent(in) :: time
330
334 class(fluid_scheme_compressible_t), intent(inout) :: this
335 integer :: n
336
337 n = this%S%dof%size()
338
340 if (neko_bcknd_device .eq. 1) then
341 call compressible_ops_device_compute_entropy(this%S, this%p, this%rho, this%gamma, n)
342 else
343 call compressible_ops_cpu_compute_entropy(this%S%x, this%p%x, this%rho%x, this%gamma, n)
344 end if
345
347
351 class(fluid_scheme_compressible_t), intent(inout) :: this
352 integer :: n
353
354 n = this%u%dof%size()
355
357 if (neko_bcknd_device .eq. 1) then
358 call compressible_ops_device_compute_max_wave_speed(this%max_wave_speed, &
359 this%u, this%v, this%w, this%gamma, this%p, this%rho, n)
360 else
361 call compressible_ops_cpu_compute_max_wave_speed(this%max_wave_speed%x, &
362 this%u%x, this%v%x, this%w%x, this%gamma, this%p%x, this%rho%x, n)
363 end if
364
366
372 subroutine fluid_scheme_compressible_log_solver_info(this, params, scheme, lx)
373 class(fluid_scheme_compressible_t), intent(inout) :: this
374 type(json_file), intent(inout) :: params
375 character(len=*), intent(in) :: scheme
376 integer, intent(in) :: lx
377 character(len=LOG_SIZE) :: log_buf
378 logical :: logical_val
379 real(kind=rp) :: real_val
380 integer :: integer_val
381
382 call neko_log%section('Fluid')
383 write(log_buf, '(A, A)') 'Type : ', trim(scheme)
384 call neko_log%message(log_buf)
385 write(log_buf, '(A, A)') 'Name : ', trim(this%name)
386 call neko_log%message(log_buf)
387
388 ! Polynomial order
389 if (lx .lt. 10) then
390 write(log_buf, '(A, I1)') 'Poly order : ', lx-1
391 else if (lx .ge. 10) then
392 write(log_buf, '(A, I2)') 'Poly order : ', lx-1
393 else
394 write(log_buf, '(A, I3)') 'Poly order : ', lx-1
395 end if
396 call neko_log%message(log_buf)
397
398 ! Global points information
399 write(log_buf, '(A, I0)') 'GLL points : ', this%glb_n_points
400 call neko_log%message(log_buf)
401 write(log_buf, '(A, I0)') 'Unique pts.: ', this%glb_unique_points
402 call neko_log%message(log_buf)
403
404 ! Material properties
405 write(log_buf, '(A,ES13.6)') 'gamma :', this%gamma
406 call neko_log%message(log_buf)
407
408 ! Compressible-specific parameters
409 call json_get_or_default(params, 'case.numerics.c_avisc_low', real_val, 0.5_rp)
410 write(log_buf, '(A,ES13.6)') 'c_avisc_low:', real_val
411 call neko_log%message(log_buf)
412
413 call json_get_or_default(params, 'case.numerics.time_order', integer_val, 4)
414 write(log_buf, '(A, I0)') 'RK order : ', integer_val
415 call neko_log%message(log_buf)
416
418
Abstract interface to sets rho and mu.
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
CPU implementation of compressible flow operations.
subroutine, public compressible_ops_cpu_compute_entropy(s, p, rho, gamma, n)
Compute entropy field S = 1/(gamma-1) * rho * (log(p) - gamma * log(rho)) on CPU.
subroutine, public compressible_ops_cpu_compute_max_wave_speed(max_wave_speed, u, v, w, gamma, p, rho, n)
Compute maximum wave speed for compressible flows on CPU.
Device implementation of compressible flow operations.
subroutine, public compressible_ops_device_compute_entropy(s, p, rho, gamma, n)
Compute entropy field S = 1/(gamma-1) * rho * (log(p) - gamma * log(rho)) on device.
subroutine, public compressible_ops_device_compute_max_wave_speed(max_wave_speed, u, v, w, gamma, p, rho, n)
Compute maximum wave speed for compressible flows on device.
subroutine, public field_col2(a, b, n)
Vector multiplication .
subroutine, public field_cmult2(a, b, c, n)
Multiplication by constant c .
subroutine, public field_cfill(a, c, n)
Set all elements to a constant c .
subroutine, public field_addcol3(a, b, c, n)
Returns .
subroutine, public field_add2(a, b, n)
Vector addition .
subroutine, public field_col3(a, b, c, n)
Vector multiplication with 3 vectors .
subroutine, public field_cmult(a, c, n)
Multiplication by constant c .
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
subroutine fluid_scheme_compressible_log_solver_info(this, params, scheme, lx)
Log comprehensive solver information.
subroutine fluid_scheme_compressible_free(this)
Free allocated memory and cleanup resources.
subroutine fluid_scheme_compressible_update_material_properties(this, time)
Set rho and mu.
subroutine fluid_scheme_compressible_init(this, msh, lx, params, scheme, user)
Initialize common data for compressible fluid scheme.
real(kind=rp) function fluid_scheme_compressible_compute_cfl(this, dt)
Compute CFL number.
subroutine fluid_scheme_compressible_compute_max_wave_speed(this)
Compute maximum wave speed for compressible flows.
subroutine fluid_scheme_compressible_validate(this)
Validate field initialization and compute derived quantities.
subroutine fluid_scheme_compressible_compute_entropy(this)
Compute entropy field S = 1/(gamma-1) * rho * (log(p) - gamma * log(rho))
Utilities for retrieving parameters from the case files.
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
Definition math.f90:60
real(kind=rp) function, public glsum(a, n)
Sum a vector of length n.
Definition math.f90:499
Defines a mesh.
Definition mesh.f90:34
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public i8
Definition num_types.f90:7
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
real(kind=rp) function, public cfl_compressible(dt, max_wave_speed, xh, coef, nelv, gdim)
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
Defines a function space.
Definition space.f90:34
integer, parameter, public gll
Definition space.f90:48
Module with things related to the simulation time.
Interfaces for user interaction with NEKO.
Definition user_intf.f90:34
Base type of all fluid formulations.
A struct that contains all info about the time, expand as needed.
A type collecting all the overridable user routines and flag to suppress type injection from custom m...