Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
les_model.f90
Go to the documentation of this file.
1! Copyright (c) 2023-2024, 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!
33!
36 use num_types, only : rp
38 use fluid_pnpn, only : fluid_pnpn_t
40 use rhs_maker, only : rhs_maker_sumab_t, rhs_maker_sumab_fctry
41 use field, only : field_t
43 use json_module, only : json_file
45 use dofmap, only : dofmap_t
46 use coefs, only : coef_t
47 use gs_ops, only : gs_op_add
50 use math, only : col2
51 use device_math, only : device_col2
53 use comm, only : pe_rank
54 implicit none
55 private
56
57 ! List of all possible types created by the factory routine
58 character(len=20) :: delta_known_types(3) = [character(len=20) :: &
59 "pointwise", &
60 "elementwise_average", &
61 "elementwise_max"]
62
64 type, abstract, public :: les_model_t
66 type(time_scheme_controller_t), pointer :: ext_bdf => null()
68 type(field_series_t), pointer :: ulag => null()
69 type(field_series_t), pointer :: vlag => null()
70 type(field_series_t), pointer :: wlag => null()
72 class(rhs_maker_sumab_t), allocatable :: sumab
74 logical :: if_ext = .false.
76 type(field_t), pointer :: nut => null()
78 character(len=:), allocatable :: delta_type
80 type(field_t), pointer :: delta => null()
82 type(coef_t), pointer :: coef => null()
83 contains
85 procedure, pass(this) :: init_base => les_model_init_base
87 procedure, pass(this) :: free_base => les_model_free_base
89 procedure, pass(this) :: compute_delta => les_model_compute_delta
91 procedure(les_model_init), pass(this), deferred :: init
93 procedure(les_model_free), pass(this), deferred :: free
95 procedure(les_model_compute), pass(this), deferred :: compute
96 end type les_model_t
97
98 abstract interface
99
102 subroutine les_model_compute(this, t, tstep)
103 import les_model_t, rp
104 class(les_model_t), intent(inout) :: this
105 real(kind=rp), intent(in) :: t
106 integer, intent(in) :: tstep
107 end subroutine les_model_compute
108 end interface
109
110 abstract interface
111
114 subroutine les_model_init(this, fluid, json)
115 import les_model_t, json_file, fluid_scheme_base_t
116 class(les_model_t), intent(inout) :: this
117 class(fluid_scheme_base_t), intent(inout), target :: fluid
118 type(json_file), intent(inout) :: json
119 end subroutine les_model_init
120 end interface
121
122 abstract interface
123
124 subroutine les_model_free(this)
125 import les_model_t
126 class(les_model_t), intent(inout) :: this
127 end subroutine les_model_free
128 end interface
129
130 interface
131
134 module subroutine les_model_allocator(object, type_name)
135 class(les_model_t), allocatable, intent(inout) :: object
136 character(len=*), intent(in) :: type_name
137 end subroutine les_model_allocator
138 end interface
139
140 interface
141
148 module subroutine les_model_factory(object, type_name, fluid, json)
149 class(les_model_t), allocatable, intent(inout) :: object
150 character(len=*), intent(in) :: type_name
151 class(fluid_scheme_base_t), intent(inout) :: fluid
152 type(json_file), intent(inout) :: json
153 end subroutine les_model_factory
154 end interface
155
156 !
157 ! Machinery for injecting user-defined types
158 !
159
163 abstract interface
164 subroutine les_model_allocate(obj)
165 import les_model_t
166 class(les_model_t), allocatable, intent(inout) :: obj
167 end subroutine les_model_allocate
168 end interface
169
170 interface
171
172 module subroutine register_les_model(type_name, allocator)
173 character(len=*), intent(in) :: type_name
174 procedure(les_model_allocate), pointer, intent(in) :: allocator
175 end subroutine register_les_model
176 end interface
177
178 ! A name-allocator pair for user-defined types. A helper type to define a
179 ! registry of custom allocators.
180 type allocator_entry
181 character(len=20) :: type_name
182 procedure(les_model_allocate), pointer, nopass :: allocator
183 end type allocator_entry
184
186 type(allocator_entry), allocatable :: les_model_registry(:)
187
189 integer :: les_model_registry_size = 0
190
191 public :: les_model_factory, les_model_allocator, register_les_model, &
192 les_model_allocate
193
194
195contains
201 subroutine les_model_init_base(this, fluid, nut_name, delta_type, if_ext)
202 class(les_model_t), intent(inout) :: this
203 class(fluid_scheme_base_t), intent(inout), target :: fluid
204 character(len=*), intent(in) :: nut_name
205 character(len=*), intent(in) :: delta_type
206 logical, intent(in) :: if_ext
207
208 associate(dofmap => fluid%dm_Xh, &
209 coef => fluid%c_Xh)
210
211 call neko_field_registry%add_field(dofmap, trim(nut_name), .true.)
212 call neko_field_registry%add_field(dofmap, "les_delta", .true.)
213 this%nut => neko_field_registry%get_field(trim(nut_name))
214 this%delta => neko_field_registry%get_field("les_delta")
215 this%coef => coef
216 this%delta_type = delta_type
217 this%if_ext = if_ext
218
219 if (pe_rank .eq. 0) then
220 if (if_ext .eqv. .true.) then
221 call neko_warning("Extrapolation of the velocity in eddy &
222 &viscosity estimation might be unstable.")
223 else
224 call neko_warning("The time integration for eddy viscosity &
225 &estimation is only first-order accurate")
226 end if
227 end if
228
229 call this%compute_delta()
230
231 select type (fluid)
232 type is (fluid_pnpn_t)
233 this%ulag => fluid%ulag
234 this%vlag => fluid%vlag
235 this%wlag => fluid%wlag
236 ! Setup backend dependent summation of AB/BDF
237 this%ext_bdf => fluid%ext_bdf
238 call rhs_maker_sumab_fctry(this%sumab)
239 class default
240 if (this%if_ext .eqv. .true.) then
241 call neko_error("Fluid scheme does not support &
242 &velocity extrapolation")
243 end if
244 end select
245
246 end associate
247 end subroutine les_model_init_base
248
250 subroutine les_model_free_base(this)
251 class(les_model_t), intent(inout) :: this
252
253 nullify(this%nut)
254 nullify(this%delta)
255 nullify(this%coef)
256 if (allocated(this%sumab)) then
257 deallocate(this%sumab)
258 end if
259 end subroutine les_model_free_base
260
267 subroutine les_model_compute_delta(this)
268 class(les_model_t), intent(inout) :: this
269 integer :: e, i, j, k
270 integer :: im, ip, jm, jp, km, kp
271 real(kind=rp) :: di, dj, dk, ndim_inv, volume_element
272 integer :: lx_half, ly_half, lz_half
273 character(len=:), allocatable :: type_string
274
275 lx_half = this%coef%Xh%lx / 2
276 ly_half = this%coef%Xh%ly / 2
277 lz_half = this%coef%Xh%lz / 2
278
279 if (this%delta_type .eq. "elementwise_max") then
280 ! use a same length scale throughout an entire element
281 ! the length scale is based on maximum GLL spacing
282 do e = 1, this%coef%msh%nelv
283 di = (this%coef%dof%x(lx_half, 1, 1, e) &
284 - this%coef%dof%x(lx_half + 1, 1, 1, e))**2 &
285 + (this%coef%dof%y(lx_half, 1, 1, e) &
286 - this%coef%dof%y(lx_half + 1, 1, 1, e))**2 &
287 + (this%coef%dof%z(lx_half, 1, 1, e) &
288 - this%coef%dof%z(lx_half + 1, 1, 1, e))**2
289
290 dj = (this%coef%dof%x(1, ly_half, 1, e) &
291 - this%coef%dof%x(1, ly_half + 1, 1, e))**2 &
292 + (this%coef%dof%y(1, ly_half, 1, e) &
293 - this%coef%dof%y(1, ly_half + 1, 1, e))**2 &
294 + (this%coef%dof%z(1, ly_half, 1, e) &
295 - this%coef%dof%z(1, ly_half + 1, 1, e))**2
296
297 dk = (this%coef%dof%x(1, 1, lz_half, e) &
298 - this%coef%dof%x(1, 1, lz_half + 1, e))**2 &
299 + (this%coef%dof%y(1, 1, lz_half, e) &
300 - this%coef%dof%y(1, 1, lz_half + 1, e))**2 &
301 + (this%coef%dof%z(1, 1, lz_half, e) &
302 - this%coef%dof%z(1, 1, lz_half + 1, e))**2
303 di = sqrt(di)
304 dj = sqrt(dj)
305 dk = sqrt(dk)
306 this%delta%x(:,:,:,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
307 end do
308 else if (this%delta_type .eq. "elementwise_average") then
309 ! use a same length scale throughout an entire element
310 ! the length scale is based on (volume)^(1/3)/(N+1)
311 do e = 1, this%coef%msh%nelv
312 volume_element = 0.0_rp
313 do k = 1, this%coef%Xh%lx * this%coef%Xh%ly * this%coef%Xh%lz
314 volume_element = volume_element + this%coef%B(k, 1, 1, e)
315 end do
316 this%delta%x(:,:,:,e) = (volume_element / &
317 (this%coef%Xh%lx - 1.0_rp) / &
318 (this%coef%Xh%ly - 1.0_rp) / &
319 (this%coef%Xh%lz - 1.0_rp) ) ** (1.0_rp / 3.0_rp)
320 end do
321 else if (this%delta_type .eq. "pointwise") then
322 do e = 1, this%coef%msh%nelv
323 do k = 1, this%coef%Xh%lz
324 km = max(1, k-1)
325 kp = min(this%coef%Xh%lz, k+1)
326
327 do j = 1, this%coef%Xh%ly
328 jm = max(1, j-1)
329 jp = min(this%coef%Xh%ly, j+1)
330
331 do i = 1, this%coef%Xh%lx
332 im = max(1, i-1)
333 ip = min(this%coef%Xh%lx, i+1)
334
335 di = (this%coef%dof%x(ip, j, k, e) - &
336 this%coef%dof%x(im, j, k, e))**2 &
337 + (this%coef%dof%y(ip, j, k, e) - &
338 this%coef%dof%y(im, j, k, e))**2 &
339 + (this%coef%dof%z(ip, j, k, e) - &
340 this%coef%dof%z(im, j, k, e))**2
341
342 dj = (this%coef%dof%x(i, jp, k, e) - &
343 this%coef%dof%x(i, jm, k, e))**2 &
344 + (this%coef%dof%y(i, jp, k, e) - &
345 this%coef%dof%y(i, jm, k, e))**2 &
346 + (this%coef%dof%z(i, jp, k, e) - &
347 this%coef%dof%z(i, jm, k, e))**2
348
349 dk = (this%coef%dof%x(i, j, kp, e) - &
350 this%coef%dof%x(i, j, km, e))**2 &
351 + (this%coef%dof%y(i, j, kp, e) - &
352 this%coef%dof%y(i, j, km, e))**2 &
353 + (this%coef%dof%z(i, j, kp, e) - &
354 this%coef%dof%z(i, j, km, e))**2
355
356 di = sqrt(di) / (ip - im)
357 dj = sqrt(dj) / (jp - jm)
358 dk = sqrt(dk) / (kp - km)
359 this%delta%x(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
360
361 end do
362 end do
363 end do
364 end do
365 else
366 call neko_type_error("delta_type for LES model", &
367 this%delta_type, delta_known_types)
368 stop
369 end if
370
371 if (neko_bcknd_device .eq. 1) then
372 call device_memcpy(this%delta%x, this%delta%x_d, this%delta%dof%size(),&
373 host_to_device, sync = .false.)
374 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
375 call device_col2(this%delta%x_d, this%coef%mult_d, this%delta%dof%size())
376 else
377 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
378 call col2(this%delta%x, this%coef%mult, this%delta%dof%size())
379 end if
380
381 end subroutine les_model_compute_delta
382
383end module les_model
Copy data between host and device (or device and device)
Definition device.F90:66
Compute eddy viscosity.
Common constructor.
Coefficients.
Definition coef.f90:34
Definition comm.F90:1
integer, public pe_rank
MPI rank.
Definition comm.F90:55
subroutine, public device_col2(a_d, b_d, n, strm)
Vector multiplication .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Contains the field_serties_t type.
Defines a field.
Definition field.f90:34
Modular version of the Classic Nek5000 Pn/Pn formulation for fluids.
Defines Gather-scatter operations.
Definition gs_ops.f90:34
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
Implements les_model_t.
Definition les_model.f90:35
character(len=20), dimension(3) delta_known_types
Definition les_model.f90:58
subroutine les_model_compute_delta(this)
Compute the LES lengthscale. For each GLL point, we take the distance between its neighbours in all 3...
subroutine les_model_init_base(this, fluid, nut_name, delta_type, if_ext)
Constructor for the les_model_t (base) class.
subroutine les_model_free_base(this)
Destructor for the les_model_t (base) class.
Definition math.f90:60
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:854
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Routines to generate the right-hand sides for the convection-diffusion equation. Employs the EXT/BDF ...
Definition rhs_maker.f90:38
Compound scheme for the advection and diffusion operators in a transport equation.
Utilities.
Definition utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition utils.f90:346
subroutine, public neko_type_error(base_type, wrong_type, known_types)
Reports an error allocating a type for a particular base pointer class.
Definition utils.f90:313
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
Stores a series (sequence) of fields, logically connected to a base field, and arranged according to ...
Base type of all fluid formulations.
Base abstract type for LES models based on the Boussinesq approximation.
Definition les_model.f90:64
Abstract type to compute extrapolated velocity field for the pressure equation.
Definition rhs_maker.f90:46
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...
#define max(a, b)
Definition tensor.cu:40