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 / this%coef%Xh%lx &
317 / this%coef%Xh%ly / this%coef%Xh%lz)**(1.0_rp / 3.0_rp)
318 end do
319 else if (this%delta_type .eq. "pointwise") then
320 do e = 1, this%coef%msh%nelv
321 do k = 1, this%coef%Xh%lz
322 km = max(1, k-1)
323 kp = min(this%coef%Xh%lz, k+1)
324
325 do j = 1, this%coef%Xh%ly
326 jm = max(1, j-1)
327 jp = min(this%coef%Xh%ly, j+1)
328
329 do i = 1, this%coef%Xh%lx
330 im = max(1, i-1)
331 ip = min(this%coef%Xh%lx, i+1)
332
333 di = (this%coef%dof%x(ip, j, k, e) - &
334 this%coef%dof%x(im, j, k, e))**2 &
335 + (this%coef%dof%y(ip, j, k, e) - &
336 this%coef%dof%y(im, j, k, e))**2 &
337 + (this%coef%dof%z(ip, j, k, e) - &
338 this%coef%dof%z(im, j, k, e))**2
339
340 dj = (this%coef%dof%x(i, jp, k, e) - &
341 this%coef%dof%x(i, jm, k, e))**2 &
342 + (this%coef%dof%y(i, jp, k, e) - &
343 this%coef%dof%y(i, jm, k, e))**2 &
344 + (this%coef%dof%z(i, jp, k, e) - &
345 this%coef%dof%z(i, jm, k, e))**2
346
347 dk = (this%coef%dof%x(i, j, kp, e) - &
348 this%coef%dof%x(i, j, km, e))**2 &
349 + (this%coef%dof%y(i, j, kp, e) - &
350 this%coef%dof%y(i, j, km, e))**2 &
351 + (this%coef%dof%z(i, j, kp, e) - &
352 this%coef%dof%z(i, j, km, e))**2
353
354 di = sqrt(di) / (ip - im)
355 dj = sqrt(dj) / (jp - jm)
356 dk = sqrt(dk) / (kp - km)
357 this%delta%x(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
358
359 end do
360 end do
361 end do
362 end do
363 else
364 call neko_type_error("delta_type for LES model", &
365 this%delta_type, delta_known_types)
366 stop
367 end if
368
369 if (neko_bcknd_device .eq. 1) then
370 call device_memcpy(this%delta%x, this%delta%x_d, this%delta%dof%size(),&
371 host_to_device, sync = .false.)
372 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
373 call device_col2(this%delta%x_d, this%coef%mult_d, this%delta%dof%size())
374 else
375 call this%coef%gs_h%op(this%delta%x, this%delta%dof%size(), gs_op_add)
376 call col2(this%delta%x, this%coef%mult, this%delta%dof%size())
377 end if
378
379 end subroutine les_model_compute_delta
380
381end 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:860
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:284
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:251
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