35 use,
intrinsic :: iso_fortran_env, only: error_unit
51 use json_module,
only : json_file, json_core, json_value
77 type(
field_t) :: rho_res, m_x_res, m_y_res, m_z_res, m_e_res
78 type(
field_t) :: drho, dm_x, dm_y, dm_z, de
80 real(kind=
rp) :: c_avisc_low
82 class(
ax_t),
allocatable :: ax
94 procedure, pass(this) :: setup_bcs &
107 module subroutine density_bc_factory(object, scheme, json, coef, user)
108 class(
bc_t),
pointer,
intent(inout) :: object
109 type(fluid_scheme_compressible_euler_t),
intent(in) :: scheme
110 type(json_file),
intent(inout) :: json
111 type(
coef_t),
intent(in) :: coef
112 type(
user_t),
intent(in) :: user
113 end subroutine density_bc_factory
124 module subroutine pressure_bc_factory(object, scheme, json, coef, user)
125 class(
bc_t),
pointer,
intent(inout) :: object
126 type(fluid_scheme_compressible_euler_t),
intent(in) :: scheme
127 type(json_file),
intent(inout) :: json
128 type(
coef_t),
intent(in) :: coef
129 type(
user_t),
intent(in) :: user
130 end subroutine pressure_bc_factory
141 module subroutine velocity_bc_factory(object, scheme, json, coef, user)
142 class(
bc_t),
pointer,
intent(inout) :: object
143 type(fluid_scheme_compressible_euler_t),
intent(in) :: scheme
144 type(json_file),
intent(inout) :: json
145 type(
coef_t),
intent(in) :: coef
146 type(
user_t),
intent(in) :: user
147 end subroutine velocity_bc_factory
157 subroutine fluid_scheme_compressible_euler_init(this, msh, lx, params, user)
158 class(fluid_scheme_compressible_euler_t),
target,
intent(inout) :: this
159 type(
mesh_t),
target,
intent(inout) :: msh
160 integer,
intent(in) :: lx
161 type(json_file),
target,
intent(inout) :: params
162 type(
user_t),
target,
intent(in) :: user
163 character(len=12),
parameter :: scheme =
'compressible'
169 call this%scheme_init(msh, lx, params, scheme, user)
171 call euler_rhs_factory(this%euler_rhs)
173 associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
174 dm_xh => this%dm_Xh, nelv => this%msh%nelv)
176 call this%drho%init(dm_xh,
'drho')
177 call this%dm_x%init(dm_xh,
'dm_x')
178 call this%dm_y%init(dm_xh,
'dm_y')
179 call this%dm_z%init(dm_xh,
'dm_z')
180 call this%dE%init(dm_xh,
'dE')
181 call this%h%init(dm_xh,
'h')
186 associate(p => this%p, rho_field => this%rho_field, &
187 u => this%u, v => this%v, w => this%w, &
188 m_x => this%m_x, m_y => this%m_y, m_z => this%m_z)
191 call device_memcpy(rho_field%x, rho_field%x_d, rho_field%dof%size(), &
209 call ax_helm_factory(this%Ax, full_formulation = .false.)
212 call this%compute_h()
213 call json_get_or_default(params,
'case.numerics.c_avisc_low', &
214 this%c_avisc_low, 0.5_rp)
217 call json_get_or_default(params,
'case.numerics.time_order', rk_order, 4)
218 call this%rk_scheme%init(rk_order)
221 call this%setup_bcs(user, params)
157 subroutine fluid_scheme_compressible_euler_init(this, msh, lx, params, user)
…
223 end subroutine fluid_scheme_compressible_euler_init
227 subroutine fluid_scheme_compressible_euler_free(this)
228 class(fluid_scheme_compressible_euler_t),
intent(inout) :: this
230 if (
allocated(this%Ax))
then
234 call this%drho%free()
235 call this%dm_x%free()
236 call this%dm_y%free()
237 call this%dm_z%free()
227 subroutine fluid_scheme_compressible_euler_free(this)
…
241 end subroutine fluid_scheme_compressible_euler_free
250 subroutine fluid_scheme_compressible_euler_step(this, t, tstep, dt, &
251 ext_bdf, dt_controller)
252 class(fluid_scheme_compressible_euler_t),
target,
intent(inout) :: this
253 real(kind=rp),
intent(in) :: t
254 integer,
intent(in) :: tstep
255 real(kind=rp),
intent(in) :: dt
256 type(time_scheme_controller_t),
intent(in) :: ext_bdf
257 type(time_step_controller_t),
intent(in) :: dt_controller
258 type(field_t),
pointer :: temp
259 integer :: temp_indices(1)
263 n = this%dm_Xh%size()
264 call this%scratch%request_field(temp, temp_indices(1))
266 call profiler_start_region(
'Fluid compressible', 1)
267 associate(u => this%u, v => this%v, w => this%w, p => this%p, &
268 m_x=> this%m_x, m_y => this%m_y, m_z => this%m_z, &
269 xh => this%Xh, msh => this%msh, ax => this%Ax, &
270 c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
271 rho => this%rho, mu => this%mu, e => this%E, &
272 rho_field => this%rho_field, mu_field => this%mu_field, &
273 ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
274 f_x => this%f_x, f_y => this%f_y, f_z => this%f_z, &
275 drho => this%drho, dm_x => this%dm_x, dm_y => this%dm_y, &
276 dm_z => this%dm_z, de => this%dE, &
277 euler_rhs => this%euler_rhs, h => this%h, &
278 c_avisc_low => this%c_avisc_low, rk_scheme => this%rk_scheme)
283 call euler_rhs%step(rho_field, m_x, m_y, m_z, e, &
285 c_xh, gs_xh, h, c_avisc_low, &
289 call this%bcs_density%apply(rho_field, t, tstep)
293 call field_copy(u, m_x, n)
294 call field_invcol2(u, rho_field, n)
295 call field_copy(v, m_y, n)
296 call field_invcol2(v, rho_field, n)
297 call field_copy(w, m_z, n)
298 call field_invcol2(w, rho_field, n)
301 call this%bcs_vel%apply_vector(u%x, v%x, w%x, &
302 dm_xh%size(), t, tstep, strong = .true.)
303 call field_copy(m_x, u, n)
304 call field_col2(m_x, rho_field, n)
305 call field_copy(m_y, v, n)
306 call field_col2(m_y, rho_field, n)
307 call field_copy(m_z, w, n)
308 call field_col2(m_z, rho_field, n)
311 call field_col3(temp, u, u, n)
312 call field_addcol3(temp, v, v, n)
313 call field_addcol3(temp, w, w, n)
314 call field_col2(temp, rho_field, n)
315 call field_cmult(temp, 0.5_rp, n)
316 call field_copy(p, e, n)
317 call field_sub2(p, temp, n)
318 call field_cmult(p, this%gamma - 1.0_rp, n)
321 call this%bcs_prs%apply(p, t, tstep)
324 call field_copy(e, p, n)
325 call field_cmult(e, 1.0_rp / (this%gamma - 1.0_rp), n)
327 call field_add2(e, temp, n)
335 call profiler_end_region(
'Fluid compressible', 1)
337 call this%scratch%relinquish_field(temp_indices)
250 subroutine fluid_scheme_compressible_euler_step(this, t, tstep, dt, &
…
339 end subroutine fluid_scheme_compressible_euler_step
345 subroutine fluid_scheme_compressible_euler_setup_bcs(this, user, params)
346 class(fluid_scheme_compressible_euler_t),
intent(inout) :: this
347 type(user_t),
target,
intent(in) :: user
348 type(json_file),
intent(inout) :: params
349 integer :: i, n_bcs, zone_index, j, zone_size, global_zone_size, ierr
350 class(bc_t),
pointer :: bc_i
351 type(json_core) :: core
352 type(json_value),
pointer :: bc_object
353 type(json_file) :: bc_subdict
355 integer,
allocatable :: zone_indices(:)
356 character(len=LOG_SIZE) :: log_buf
359 if (params%valid_path(
'case.fluid.boundary_conditions'))
then
360 call params%info(
'case.fluid.boundary_conditions', n_children = n_bcs)
361 call params%get_core(core)
362 call params%get(
'case.fluid.boundary_conditions', bc_object, found)
367 call this%bcs_vel%init(n_bcs)
371 call json_extract_item(core, bc_object, i, bc_subdict)
372 call json_get(bc_subdict,
"zone_indices", zone_indices)
375 do j = 1,
size(zone_indices)
376 zone_size = this%msh%labeled_zones(zone_indices(j))%size
377 call mpi_allreduce(zone_size, global_zone_size, 1, &
378 mpi_integer, mpi_max, neko_comm, ierr)
380 if (global_zone_size .eq. 0)
then
381 write(error_unit,
'(A,I0,A)')
"Error: Zone ", zone_indices(j), &
389 call velocity_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
392 if (
associated(bc_i))
then
393 call this%bcs_vel%append(bc_i)
400 call this%bcs_prs%init(n_bcs)
404 call json_extract_item(core, bc_object, i, bc_subdict)
406 call pressure_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
410 if (
associated(bc_i))
then
411 call this%bcs_prs%append(bc_i)
418 call this%bcs_density%init(n_bcs)
422 call json_extract_item(core, bc_object, i, bc_subdict)
424 call density_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
428 if (
associated(bc_i))
then
429 call this%bcs_density%append(bc_i)
345 subroutine fluid_scheme_compressible_euler_setup_bcs(this, user, params)
…
433 end subroutine fluid_scheme_compressible_euler_setup_bcs
439 subroutine compute_h(this)
440 class(fluid_scheme_compressible_euler_t),
intent(inout) :: this
441 integer :: e, i, j, k
442 integer :: im, ip, jm, jp, km, kp
443 real(kind=rp) :: di, dj, dk, ndim_inv
444 integer :: lx_half, ly_half, lz_half
446 lx_half = this%c_Xh%Xh%lx / 2
447 ly_half = this%c_Xh%Xh%ly / 2
448 lz_half = this%c_Xh%Xh%lz / 2
450 do e = 1, this%c_Xh%msh%nelv
451 do k = 1, this%c_Xh%Xh%lz
453 kp = min(this%c_Xh%Xh%lz, k+1)
455 do j = 1, this%c_Xh%Xh%ly
457 jp = min(this%c_Xh%Xh%ly, j+1)
459 do i = 1, this%c_Xh%Xh%lx
461 ip = min(this%c_Xh%Xh%lx, i+1)
463 di = (this%c_Xh%dof%x(ip, j, k, e) - &
464 this%c_Xh%dof%x(im, j, k, e))**2 &
465 + (this%c_Xh%dof%y(ip, j, k, e) - &
466 this%c_Xh%dof%y(im, j, k, e))**2 &
467 + (this%c_Xh%dof%z(ip, j, k, e) - &
468 this%c_Xh%dof%z(im, j, k, e))**2
470 dj = (this%c_Xh%dof%x(i, jp, k, e) - &
471 this%c_Xh%dof%x(i, jm, k, e))**2 &
472 + (this%c_Xh%dof%y(i, jp, k, e) - &
473 this%c_Xh%dof%y(i, jm, k, e))**2 &
474 + (this%c_Xh%dof%z(i, jp, k, e) - &
475 this%c_Xh%dof%z(i, jm, k, e))**2
477 dk = (this%c_Xh%dof%x(i, j, kp, e) - &
478 this%c_Xh%dof%x(i, j, km, e))**2 &
479 + (this%c_Xh%dof%y(i, j, kp, e) - &
480 this%c_Xh%dof%y(i, j, km, e))**2 &
481 + (this%c_Xh%dof%z(i, j, kp, e) - &
482 this%c_Xh%dof%z(i, j, km, e))**2
484 di = sqrt(di) / (ip - im)
485 dj = sqrt(dj) / (jp - jm)
486 dk = sqrt(dk) / (kp - km)
487 this%h%x(i,j,k,e) = (di * dj * dk)**(1.0_rp / 3.0_rp)
494 if (neko_bcknd_device .eq. 1)
then
495 call device_memcpy(this%h%x, this%h%x_d, this%h%dof%size(),&
496 host_to_device, sync = .false.)
497 call this%gs_Xh%op(this%h, gs_op_add)
498 call device_col2(this%h%x_d, this%c_Xh%mult_d, this%h%dof%size())
500 call this%gs_Xh%op(this%h, gs_op_add)
501 call col2(this%h%x, this%c_Xh%mult, this%h%dof%size())
439 subroutine compute_h(this)
…
504 end subroutine compute_h
510 subroutine fluid_scheme_compressible_euler_restart(this, dtlag, tlag)
511 class(fluid_scheme_compressible_euler_t),
target,
intent(inout) :: this
512 real(kind=rp) :: dtlag(10), tlag(10)
510 subroutine fluid_scheme_compressible_euler_restart(this, dtlag, tlag)
…
513 end subroutine fluid_scheme_compressible_euler_restart
Copy data between host and device (or device and device)
Abstract interface to evaluate rhs.
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Retrieves a parameter by name or throws an error.
Subroutines to add advection terms to the RHS of a transport equation.
Defines a Matrix-vector product.
Defines a boundary condition.
subroutine, public device_col2(a_d, b_d, n)
Vector multiplication .
Device abstraction, common interface for various accelerators.
integer, parameter, public host_to_device
Defines a mapping of the degrees of freedom.
Defines inflow dirichlet conditions.
Defines user dirichlet condition for a scalar field.
subroutine, public field_invcol2(a, b, n)
Vector division .
subroutine, public field_cadd(a, s, n)
Add a scalar to vector .
subroutine, public field_col2(a, b, n)
Vector multiplication .
subroutine, public field_sub2(a, b, n)
Vector substraction .
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_copy(a, b, n)
Copy a vector .
subroutine, public field_cmult(a, c, n)
Multiplication by constant c .
subroutine compute_h(this)
Copied from les_model_compute_delta in les_model.f90 TODO: move to a separate module Compute characte...
subroutine fluid_scheme_compressible_euler_init(this, msh, lx, params, user)
Boundary condition factory for density.
subroutine fluid_scheme_compressible_euler_step(this, t, tstep, dt, ext_bdf, dt_controller)
Advance the fluid simulation one timestep.
subroutine fluid_scheme_compressible_euler_restart(this, dtlag, tlag)
Restart the simulation from saved state.
subroutine fluid_scheme_compressible_euler_free(this)
Free allocated memory and cleanup.
subroutine fluid_scheme_compressible_euler_setup_bcs(this, user, params)
Set up boundary conditions for the fluid scheme.
Defines Gather-scatter operations.
integer, parameter, public gs_op_add
Utilities for retrieving parameters from the case files.
integer, parameter, public log_size
subroutine, public subcol3(a, b, c, n)
Returns .
subroutine, public addcol3(a, b, c, n)
Returns .
subroutine, public col2(a, b, n)
Vector multiplication .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
subroutine, public div(res, ux, uy, uz, coef)
Compute the divergence of a vector field.
subroutine, public grad(ux, uy, uz, u, coef)
Compute the gradient of a scalar field.
subroutine, public profiler_start_region(name, region_id)
Started a named (name) profiler region.
subroutine, public profiler_end_region(name, region_id)
End the most recently started profiler region.
Defines a function space.
Compound scheme for the advection and diffusion operators in a transport equation.
Implements type time_step_controller.
Interfaces for user interaction with NEKO.
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Defines a zero-valued Dirichlet boundary condition.
Base abstract type for computing the advection operator.
Base type for a matrix-vector product providing .
Base type for a boundary condition.
A list of allocatable `bc_t`. Follows the standard interface of lists.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Abstract type to compute rhs.
User defined dirichlet condition, for which the user can work with an entire field....
Extension of the user defined dirichlet condition field_dirichlet
field_list_t, To be able to group fields together
Base type of compressible fluid formulations.
The function space for the SEM solution fields.
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...
Provides a tool to set time step dt.
Zero-valued Dirichlet boundary condition. Used for no-slip walls, but also for various auxillary cond...