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
fluid_scheme_compressible_euler.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 comm
35 use, intrinsic :: iso_fortran_env, only: error_unit
36 use advection, only : advection_t, advection_factory
38 use dofmap, only : dofmap_t
42 use math, only : col2, copy, col3, addcol3, subcol3
43 use device_math, only : device_col2
44 use field, only : field_t
46 use gs_ops, only : gs_op_add
47 use gather_scatter, only : gs_t
48 use num_types, only : rp
49 use mesh, only : mesh_t
50 use operators, only: div, grad
51 use json_module, only : json_file, json_core, json_value
54 use user_intf, only : user_t
57 use ax_product, only : ax_t, ax_helm_factory
58 use field_list, only : field_list_t
59 use coefs, only: coef_t
60 use space, only : space_t
61 use euler_residual, only: euler_rhs_t, euler_rhs_factory
66 use bc_list, only: bc_list_t
68 use field_math, only : field_copy
69 use bc, only : bc_t
70 use utils, only : neko_error, neko_warning
71 use logger, only : log_size
72 implicit none
73 private
74
75 type, public, extends(fluid_scheme_compressible_t) &
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
79 type(field_t) :: h
80 real(kind=rp) :: c_avisc_low
81 class(advection_t), allocatable :: adv
82 class(ax_t), allocatable :: ax
83 class(euler_rhs_t), allocatable :: euler_rhs
84 type(runge_kutta_time_scheme_t) :: rk_scheme
85
86 ! List of boundary conditions for velocity
87 type(bc_list_t) :: bcs_density
88 contains
89 procedure, pass(this) :: init => fluid_scheme_compressible_euler_init
90 procedure, pass(this) :: free => fluid_scheme_compressible_euler_free
91 procedure, pass(this) :: step => fluid_scheme_compressible_euler_step
92 procedure, pass(this) :: restart => fluid_scheme_compressible_euler_restart
94 procedure, pass(this) :: setup_bcs &
96 procedure, pass(this) :: compute_h
98
99 interface
100
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
114 end interface
115
116 interface
117
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
131 end interface
132
133 interface
134
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
148 end interface
149
150contains
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'
164 integer :: rk_order
165
166 call this%free()
167
168 ! Initialize base class
169 call this%scheme_init(msh, lx, params, scheme, user)
170
171 call euler_rhs_factory(this%euler_rhs)
172
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)
175
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')
182
183 end associate
184
185 if (neko_bcknd_device .eq. 1) then
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)
189 call device_memcpy(p%x, p%x_d, p%dof%size(), &
190 host_to_device, sync = .false.)
191 call device_memcpy(rho_field%x, rho_field%x_d, rho_field%dof%size(), &
192 host_to_device, sync = .false.)
193 call device_memcpy(u%x, u%x_d, u%dof%size(), &
194 host_to_device, sync = .false.)
195 call device_memcpy(v%x, v%x_d, v%dof%size(), &
196 host_to_device, sync = .false.)
197 call device_memcpy(w%x, w%x_d, w%dof%size(), &
198 host_to_device, sync = .false.)
199 call device_memcpy(m_x%x, m_x%x_d, m_x%dof%size(), &
200 host_to_device, sync = .false.)
201 call device_memcpy(m_y%x, m_y%x_d, m_y%dof%size(), &
202 host_to_device, sync = .false.)
203 call device_memcpy(m_z%x, m_z%x_d, m_z%dof%size(), &
204 host_to_device, sync = .false.)
205 end associate
206 end if
207
208 ! Initialize the diffusion operator
209 call ax_helm_factory(this%Ax, full_formulation = .false.)
210
211 ! Compute h
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)
215
216 ! Initialize Runge-Kutta scheme
217 call json_get_or_default(params, 'case.numerics.time_order', rk_order, 4)
218 call this%rk_scheme%init(rk_order)
219
220 ! Set up boundary conditions
221 call this%setup_bcs(user, params)
222
223 end subroutine fluid_scheme_compressible_euler_init
224
227 subroutine fluid_scheme_compressible_euler_free(this)
228 class(fluid_scheme_compressible_euler_t), intent(inout) :: this
229
230 if (allocated(this%Ax)) then
231 deallocate(this%Ax)
232 end if
233
234 call this%drho%free()
235 call this%dm_x%free()
236 call this%dm_y%free()
237 call this%dm_z%free()
238 call this%dE%free()
239
240 ! call this%scheme_free()
241 end subroutine fluid_scheme_compressible_euler_free
242
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)
260 ! number of degrees of freedom
261 integer :: n
262
263 n = this%dm_Xh%size()
264 call this%scratch%request_field(temp, temp_indices(1))
265
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)
279
280 ! Hack: If m_z is always zero, use it to visualize rho
281 ! call field_cfill(m_z, 0.0_rp, n)
282
283 call euler_rhs%step(rho_field, m_x, m_y, m_z, e, &
284 p, u, v, w, ax, &
285 c_xh, gs_xh, h, c_avisc_low, &
286 rk_scheme, dt)
287
289 call this%bcs_density%apply(rho_field, t, tstep)
290
292 ! Update u, v, w
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)
299
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)
309
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)
319
321 call this%bcs_prs%apply(p, t, tstep)
322 ! TODO: Make sure pressure is positive
323 ! E = p / (gamma - 1) + 0.5 * rho * (u^2 + v^2 + w^2)
324 call field_copy(e, p, n)
325 call field_cmult(e, 1.0_rp / (this%gamma - 1.0_rp), n)
326 ! temp = 0.5 * rho * (u^2 + v^2 + w^2)
327 call field_add2(e, temp, n)
328
330
331 ! Hack: If m_z is always zero, use it to visualize rho
332 ! call field_copy(w, rho_field, n)
333
334 end associate
335 call profiler_end_region('Fluid compressible', 1)
336
337 call this%scratch%relinquish_field(temp_indices)
338
339 end subroutine fluid_scheme_compressible_euler_step
340
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
354 logical :: found
355 integer, allocatable :: zone_indices(:)
356 character(len=LOG_SIZE) :: log_buf
357
358 ! Process boundary conditions
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)
363
364 !
365 ! Velocity bcs
366 !
367 call this%bcs_vel%init(n_bcs)
368
369 do i = 1, n_bcs
370 ! Extract BC configuration
371 call json_extract_item(core, bc_object, i, bc_subdict)
372 call json_get(bc_subdict, "zone_indices", zone_indices)
373
374 ! Validate zones
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)
379
380 if (global_zone_size .eq. 0) then
381 write(error_unit,'(A,I0,A)') "Error: Zone ", zone_indices(j), &
382 " has zero size"
383 error stop
384 end if
385 end do
386
387 ! Create BC
388 bc_i => null()
389 call velocity_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
390
391 ! Add to appropriate lists
392 if (associated(bc_i)) then
393 call this%bcs_vel%append(bc_i)
394 end if
395 end do
396
397 !
398 ! Pressure bcs
399 !
400 call this%bcs_prs%init(n_bcs)
401
402 do i = 1, n_bcs
403 ! Create a new json containing just the subdict for this bc
404 call json_extract_item(core, bc_object, i, bc_subdict)
405 bc_i => null()
406 call pressure_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
407
408 ! Not all bcs require an allocation for pressure in particular,
409 ! so we check.
410 if (associated(bc_i)) then
411 call this%bcs_prs%append(bc_i)
412 end if
413 end do
414
415 !
416 ! Density bcs
417 !
418 call this%bcs_density%init(n_bcs)
419
420 do i = 1, n_bcs
421 ! Create a new json containing just the subdict for this bc
422 call json_extract_item(core, bc_object, i, bc_subdict)
423 bc_i => null()
424 call density_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
425
426 ! Not all bcs require an allocation for pressure in particular,
427 ! so we check.
428 if (associated(bc_i)) then
429 call this%bcs_density%append(bc_i)
430 end if
431 end do
432 end if
433 end subroutine fluid_scheme_compressible_euler_setup_bcs
434
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
445
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
449
450 do e = 1, this%c_Xh%msh%nelv
451 do k = 1, this%c_Xh%Xh%lz
452 km = max(1, k-1)
453 kp = min(this%c_Xh%Xh%lz, k+1)
454
455 do j = 1, this%c_Xh%Xh%ly
456 jm = max(1, j-1)
457 jp = min(this%c_Xh%Xh%ly, j+1)
458
459 do i = 1, this%c_Xh%Xh%lx
460 im = max(1, i-1)
461 ip = min(this%c_Xh%Xh%lx, i+1)
462
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
469
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
476
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
483
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)
488
489 end do
490 end do
491 end do
492 end do
493
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())
499 else
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())
502 end if
503
504 end subroutine compute_h
505
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)
513 end subroutine fluid_scheme_compressible_euler_restart
514
Copy data between host and device (or device and device)
Definition device.F90:65
Abstract interface to evaluate rhs.
Definition euler_res.f90:54
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.
Definition advection.f90:34
Defines a Matrix-vector product.
Definition ax.f90:34
Defines a list of bc_t.
Definition bc_list.f90:34
Defines a boundary condition.
Definition bc.f90:34
Coefficients.
Definition coef.f90:34
Definition comm.F90:1
subroutine, public device_col2(a_d, b_d, n)
Vector multiplication .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:46
Defines a mapping of the degrees of freedom.
Definition dofmap.f90:35
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 .
Defines a field.
Definition field.f90:34
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.
Gather-scatter.
Defines Gather-scatter operations.
Definition gs_ops.f90:34
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
Utilities for retrieving parameters from the case files.
Logging routines.
Definition log.f90:34
integer, parameter, public log_size
Definition log.f90:42
Definition math.f90:60
subroutine, public subcol3(a, b, c, n)
Returns .
Definition math.f90:755
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
Defines a mesh.
Definition mesh.f90:34
Build configurations.
integer, parameter neko_bcknd_device
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.
subroutine, public grad(ux, uy, uz, u, coef)
Compute the gradient of a scalar field.
Profiling interface.
Definition profiler.F90:34
subroutine, public profiler_start_region(name, region_id)
Started a named (name) profiler region.
Definition profiler.F90:78
subroutine, public profiler_end_region(name, region_id)
End the most recently started profiler region.
Definition profiler.F90:109
Defines a function space.
Definition space.f90:34
Compound scheme for the advection and diffusion operators in a transport equation.
Implements type time_step_controller.
Interfaces for user interaction with NEKO.
Definition user_intf.f90:34
Utilities.
Definition utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition utils.f90:266
Defines a zero-valued Dirichlet boundary condition.
Base abstract type for computing the advection operator.
Definition advection.f90:46
Base type for a matrix-vector product providing .
Definition ax.f90:43
Base type for a boundary condition.
Definition bc.f90:57
A list of allocatable `bc_t`. Follows the standard interface of lists.
Definition bc_list.f90:47
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
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
The function space for the SEM solution fields.
Definition space.f90:62
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...
Zero-valued Dirichlet boundary condition. Used for no-slip walls, but also for various auxillary cond...
#define max(a, b)
Definition tensor.cu:40