Neko  0.9.0
A portable framework for high-order spectral element flow simulations
fluid_scheme.f90
Go to the documentation of this file.
1 ! Copyright (c) 2020-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 !
35  use gather_scatter, only : gs_t
36  use mean_sqr_flow, only : mean_sqr_flow_t
37  use neko_config, only : neko_bcknd_device
38  use checkpoint, only : chkp_t
39  use mean_flow, only : mean_flow_t
40  use num_types, only : rp
41  use comm
43  use field, only : field_t
44  use space, only : space_t, gll
45  use dofmap, only : dofmap_t
46  use krylov, only : ksp_t, krylov_solver_factory, krylov_solver_destroy, &
48  use coefs, only: coef_t
49  use wall, only : no_slip_wall_t
50  use inflow, only : inflow_t
52  use blasius, only : blasius_t
53  use dirichlet, only : dirichlet_t
54  use dong_outflow, only : dong_outflow_t
55  use symmetry, only : symmetry_t
58  use jacobi, only : jacobi_t
59  use sx_jacobi, only : sx_jacobi_t
60  use device_jacobi, only : device_jacobi_t
61  use hsmg, only : hsmg_t
62  use precon, only : pc_t, precon_factory, precon_destroy
63  use fluid_stats, only : fluid_stats_t
67  use math, only : cfill, add2s2
70  use operators, only : cfl
74  use json_module, only : json_file
78  use utils, only : neko_error, neko_warning
79  use field_series, only : field_series_t
81  use field_math, only : field_cfill
82  use wall_model_bc, only : wall_model_bc_t
83  use shear_stress, only : shear_stress_t
85  implicit none
86  private
87 
89  type, abstract :: fluid_scheme_t
90  type(field_t), pointer :: u => null()
91  type(field_t), pointer :: v => null()
92  type(field_t), pointer :: w => null()
93  type(field_t), pointer :: p => null()
94  type(field_series_t) :: ulag, vlag, wlag
95  type(space_t) :: xh
96  type(dofmap_t) :: dm_xh
97  type(gs_t) :: gs_xh
98  type(coef_t) :: c_xh
102  type(field_t), pointer :: f_x => null()
104  type(field_t), pointer :: f_y => null()
106  type(field_t), pointer :: f_z => null()
107  class(ksp_t), allocatable :: ksp_vel
108  class(ksp_t), allocatable :: ksp_prs
109  class(pc_t), allocatable :: pc_vel
110  class(pc_t), allocatable :: pc_prs
111  integer :: vel_projection_dim
112  integer :: pr_projection_dim
113  integer :: vel_projection_activ_step
114  integer :: pr_projection_activ_step
115  type(no_slip_wall_t) :: bc_wall
116  class(bc_t), allocatable :: bc_inflow
117  type(wall_model_bc_t) :: bc_wallmodel
119  logical :: if_gradient_jump_penalty
120  type(gradient_jump_penalty_t) :: gradient_jump_penalty_u
121  type(gradient_jump_penalty_t) :: gradient_jump_penalty_v
122  type(gradient_jump_penalty_t) :: gradient_jump_penalty_w
123 
124  ! Attributes for field dirichlet BCs
125  type(field_dirichlet_vector_t) :: user_field_bc_vel
126  type(field_dirichlet_t) :: user_field_bc_prs
127  type(dirichlet_t) :: bc_prs
128  type(dong_outflow_t) :: bc_dong
129  type(symmetry_t) :: bc_sym
130  type(shear_stress_t) :: bc_sh
131  type(bc_list_t) :: bclst_vel
132  type(bc_list_t) :: bclst_vel_neumann
133  type(bc_list_t) :: bclst_prs
134  type(field_t) :: bdry
135  type(json_file), pointer :: params
136  type(mesh_t), pointer :: msh => null()
137  type(chkp_t) :: chkp
138  type(mean_flow_t) :: mean
140  type(mean_sqr_flow_t) :: mean_sqr
141  logical :: forced_flow_rate = .false.
142  logical :: freeze = .false.
144  real(kind=rp) :: mu
146  type(field_t) :: mu_field
148  character(len=:), allocatable :: nut_field_name
150  logical :: variable_material_properties = .false.
152  real(kind=rp) :: rho
154  type(field_t) :: rho_field
155  type(scratch_registry_t) :: scratch
157  character(len=NEKO_MSH_MAX_ZLBL_LEN), allocatable :: bc_labels(:)
158  contains
160  procedure, pass(this) :: fluid_scheme_init_all
161  procedure, pass(this) :: fluid_scheme_init_common
164  procedure, pass(this) :: scheme_free => fluid_scheme_free
166  procedure, pass(this) :: validate => fluid_scheme_validate
168  procedure, pass(this) :: bc_apply_vel => fluid_scheme_bc_apply_vel
170  procedure, pass(this) :: bc_apply_prs => fluid_scheme_bc_apply_prs
172  procedure, pass(this) :: set_usr_inflow => fluid_scheme_set_usr_inflow
174  procedure, pass(this) :: compute_cfl => fluid_compute_cfl
176  procedure, pass(this) :: set_material_properties => &
179  procedure(fluid_scheme_init_intrf), pass(this), deferred :: init
181  procedure(fluid_scheme_free_intrf), pass(this), deferred :: free
183  procedure(fluid_scheme_step_intrf), pass(this), deferred :: step
185  procedure(fluid_scheme_restart_intrf), pass(this), deferred :: restart
186  procedure, private, pass(this) :: set_bc_type_output => &
189  procedure, pass(this) :: update_material_properties => &
191  end type fluid_scheme_t
192 
194  abstract interface
195  subroutine fluid_scheme_init_intrf(this, msh, lx, params, user, &
196  time_scheme)
197  import fluid_scheme_t
198  import json_file
199  import mesh_t
200  import user_t
202  class(fluid_scheme_t), target, intent(inout) :: this
203  type(mesh_t), target, intent(inout) :: msh
204  integer, intent(inout) :: lx
205  type(json_file), target, intent(inout) :: params
206  type(user_t), target, intent(in) :: user
207  type(time_scheme_controller_t), target, intent(in) :: time_scheme
208  end subroutine fluid_scheme_init_intrf
209  end interface
210 
212  abstract interface
213  subroutine fluid_scheme_free_intrf(this)
214  import fluid_scheme_t
215  class(fluid_scheme_t), intent(inout) :: this
216  end subroutine fluid_scheme_free_intrf
217  end interface
218 
220  abstract interface
221  subroutine fluid_scheme_step_intrf(this, t, tstep, dt, ext_bdf, &
222  dt_controller)
223  import fluid_scheme_t
226  import rp
227  class(fluid_scheme_t), target, intent(inout) :: this
228  real(kind=rp), intent(inout) :: t
229  integer, intent(inout) :: tstep
230  real(kind=rp), intent(in) :: dt
231  type(time_scheme_controller_t), intent(inout) :: ext_bdf
232  type(time_step_controller_t), intent(in) :: dt_controller
233  end subroutine fluid_scheme_step_intrf
234  end interface
235 
237  abstract interface
238  subroutine fluid_scheme_restart_intrf(this, dtlag, tlag)
239  import fluid_scheme_t
240  import rp
241  class(fluid_scheme_t), target, intent(inout) :: this
242  real(kind=rp) :: dtlag(10), tlag(10)
243 
244  end subroutine fluid_scheme_restart_intrf
245  end interface
246 
247  interface
248 
249  module subroutine fluid_scheme_factory(object, type_name)
250  class(fluid_scheme_t), intent(inout), allocatable :: object
251  character(len=*) :: type_name
252  end subroutine fluid_scheme_factory
253  end interface
254 
255  public :: fluid_scheme_t, fluid_scheme_factory
256 
257 contains
258 
260  subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, &
261  kspv_init)
262  implicit none
263  class(fluid_scheme_t), target, intent(inout) :: this
264  type(mesh_t), target, intent(inout) :: msh
265  integer, intent(inout) :: lx
266  character(len=*), intent(in) :: scheme
267  type(json_file), target, intent(inout) :: params
268  type(user_t), target, intent(in) :: user
269  logical, intent(in) :: kspv_init
270  type(dirichlet_t) :: bdry_mask
271  character(len=LOG_SIZE) :: log_buf
272  real(kind=rp), allocatable :: real_vec(:)
273  real(kind=rp) :: real_val, kappa, b, z0
274  logical :: logical_val
275  integer :: integer_val, ierr
276  type(json_file) :: wm_json
277  character(len=:), allocatable :: string_val1, string_val2
278 
279  !
280  ! SEM simulation fundamentals
281  !
282 
283  this%msh => msh
284 
285  if (msh%gdim .eq. 2) then
286  call this%Xh%init(gll, lx, lx)
287  else
288  call this%Xh%init(gll, lx, lx, lx)
289  end if
290 
291  call this%dm_Xh%init(msh, this%Xh)
292 
293  call this%gs_Xh%init(this%dm_Xh)
294 
295  call this%c_Xh%init(this%gs_Xh)
296 
297  ! Local scratch registry
298  this%scratch = scratch_registry_t(this%dm_Xh, 10, 2)
299 
300  ! Case parameters
301  this%params => params
302 
303 
304  !
305  ! First section of fluid log
306  !
307 
308  call neko_log%section('Fluid')
309  write(log_buf, '(A, A)') 'Type : ', trim(scheme)
310  call neko_log%message(log_buf)
311 
312  !
313  ! Material properties
314  !
315  call this%set_material_properties(params, user)
316 
317  !
318  ! Turbulence modelling and variable material properties
319  !
320  if (params%valid_path('case.fluid.nut_field')) then
321  call json_get(params, 'case.fluid.nut_field', this%nut_field_name)
322  this%variable_material_properties = .true.
323  else
324  this%nut_field_name = ""
325  end if
326 
327  ! Fill mu and rho field with the physical value
328 
329  call this%mu_field%init(this%dm_Xh, "mu")
330  call this%rho_field%init(this%dm_Xh, "mu")
331  call field_cfill(this%mu_field, this%mu, this%mu_field%size())
332  call field_cfill(this%rho_field, this%rho, this%mu_field%size())
333 
334  ! Since mu, rho is a field, and the none-stress simulation fetches
335  ! data from the host arrays, we need to mirror the constant
336  ! material properties on the host
337  if (neko_bcknd_device .eq. 1) then
338  call cfill(this%mu_field%x, this%mu, this%mu_field%size())
339  call cfill(this%rho_field%x, this%rho, this%rho_field%size())
340  end if
341 
342 
343  ! Projection spaces
344  call json_get_or_default(params, &
345  'case.fluid.velocity_solver.projection_space_size', &
346  this%vel_projection_dim, 20)
347  call json_get_or_default(params, &
348  'case.fluid.pressure_solver.projection_space_size', &
349  this%pr_projection_dim, 20)
350  call json_get_or_default(params, &
351  'case.fluid.velocity_solver.projection_hold_steps', &
352  this%vel_projection_activ_step, 5)
353  call json_get_or_default(params, &
354  'case.fluid.pressure_solver.projection_hold_steps', &
355  this%pr_projection_activ_step, 5)
356 
357 
358  call json_get_or_default(params, 'case.fluid.freeze', this%freeze, .false.)
359 
360  if (params%valid_path("case.fluid.flow_rate_force")) then
361  this%forced_flow_rate = .true.
362  end if
363 
364 
365  if (lx .lt. 10) then
366  write(log_buf, '(A, I1)') 'Poly order : ', lx-1
367  else if (lx .ge. 10) then
368  write(log_buf, '(A, I2)') 'Poly order : ', lx-1
369  else
370  write(log_buf, '(A, I3)') 'Poly order : ', lx-1
371  end if
372  call neko_log%message(log_buf)
373  write(log_buf, '(A, I0)') 'DoFs : ', this%dm_Xh%size()
374  call neko_log%message(log_buf)
375 
376  write(log_buf, '(A,ES13.6)') 'rho :', this%rho
377  call neko_log%message(log_buf)
378  write(log_buf, '(A,ES13.6)') 'mu :', this%mu
379  call neko_log%message(log_buf)
380 
381  call json_get(params, 'case.numerics.dealias', logical_val)
382  write(log_buf, '(A, L1)') 'Dealias : ', logical_val
383  call neko_log%message(log_buf)
384 
385  call json_get_or_default(params, 'case.output_boundary', logical_val, &
386  .false.)
387  write(log_buf, '(A, L1)') 'Save bdry : ', logical_val
388  call neko_log%message(log_buf)
389  if (logical_val) then
390  write(log_buf, '(A)') 'bdry keys: '
391  call neko_log%message(log_buf)
392  write(log_buf, '(A)') 'No-slip wall ''w'' = 1'
393  call neko_log%message(log_buf)
394  write(log_buf, '(A)') 'Inlet/velocity dirichlet ''v'' = 2'
395  call neko_log%message(log_buf)
396  write(log_buf, '(A)') 'Outlet ''o'' = 3'
397  call neko_log%message(log_buf)
398  write(log_buf, '(A)') 'Symmetry ''sym'' = 4'
399  call neko_log%message(log_buf)
400  write(log_buf, '(A)') 'Outlet-normal ''on'' = 5'
401  call neko_log%message(log_buf)
402  write(log_buf, '(A)') 'Periodic = 6'
403  call neko_log%message(log_buf)
404  write(log_buf, '(A)') 'Dirichlet on velocity components: '
405  call neko_log%message(log_buf)
406  write(log_buf, '(A)') ' ''d_vel_u, d_vel_v, d_vel_w'' = 7'
407  call neko_log%message(log_buf)
408  write(log_buf, '(A)') 'Pressure dirichlet ''d_pres'' = 8'
409  call neko_log%message(log_buf)
410  write(log_buf, '(A)') '''d_vel_(u,v,w)'' and ''d_pres'' = 8'
411  call neko_log%message(log_buf)
412  write(log_buf, '(A)') 'Shear stress ''sh'' = 9'
413  call neko_log%message(log_buf)
414  write(log_buf, '(A)') 'Wall modelling ''wm'' = 10'
415  call neko_log%message(log_buf)
416  write(log_buf, '(A)') 'No boundary condition set = 0'
417  call neko_log%message(log_buf)
418  end if
419 
420 
421  !
422  ! Setup velocity boundary conditions
423  !
424  allocate(this%bc_labels(neko_msh_max_zlbls))
425  this%bc_labels = "not"
426 
427  if (params%valid_path('case.fluid.boundary_types')) then
428  call json_get(params, &
429  'case.fluid.boundary_types', &
430  this%bc_labels)
431  end if
432 
433  call bc_list_init(this%bclst_vel)
434 
435  call this%bc_sym%init_base(this%c_Xh)
436  call this%bc_sym%mark_zone(msh%sympln)
437  call this%bc_sym%mark_zones_from_list(msh%labeled_zones, &
438  'sym', this%bc_labels)
439  call this%bc_sym%finalize()
440  call this%bc_sym%init(this%c_Xh)
441  call bc_list_add(this%bclst_vel, this%bc_sym)
442 
443  ! Shear stress conditions
444  call this%bc_sh%init_base(this%c_Xh)
445  call this%bc_sh%mark_zones_from_list(msh%labeled_zones, &
446  'sh', this%bc_labels)
447  call this%bc_sh%finalize()
448  ! This marks the dirichlet and neumann conditions inside, and finalizes
449  ! them.
450  call this%bc_sh%init_shear_stress(this%c_Xh)
451 
452  call bc_list_add(this%bclst_vel, this%bc_sh%symmetry)
453 
454  ! Read stress value, default to [0 0 0]
455  if (this%bc_sh%msk(0) .gt. 0) then
456  call params%get('case.fluid.shear_stress.value', real_vec, logical_val)
457  if (.not. logical_val .and. this%bc_sh%msk(0) .gt. 0) then
458  call neko_warning("No stress values provided for sh boundaries, &
459  & defaulting to 0. Use fluid.shear_stress.value to set.")
460  allocate(real_vec(3))
461  real_vec = 0.0_rp
462  else if (size(real_vec) .ne. 3) then
463  call neko_error ("The shear stress vector provided in &
464  &fluid.shear_stress.value should have 3 components.")
465  end if
466  call this%bc_sh%set_stress(real_vec(1), real_vec(2), real_vec(3))
467  end if
468 
469  call bc_list_init(this%bclst_vel_neumann)
470  call bc_list_add(this%bclst_vel_neumann, this%bc_sh)
471 
472  !
473  ! Inflow
474  !
475  if (params%valid_path('case.fluid.inflow_condition')) then
476  call json_get(params, 'case.fluid.inflow_condition.type', string_val1)
477  if (trim(string_val1) .eq. "uniform") then
478  allocate(inflow_t::this%bc_inflow)
479  else if (trim(string_val1) .eq. "blasius") then
480  allocate(blasius_t::this%bc_inflow)
481  else if (trim(string_val1) .eq. "user") then
482  allocate(usr_inflow_t::this%bc_inflow)
483  else
484  call neko_error('Invalid inflow condition '//string_val1)
485  end if
486 
487  call this%bc_inflow%init_base(this%c_Xh)
488  call this%bc_inflow%mark_zone(msh%inlet)
489  call this%bc_inflow%mark_zones_from_list(msh%labeled_zones, &
490  'v', this%bc_labels)
491  call this%bc_inflow%finalize()
492  call bc_list_add(this%bclst_vel, this%bc_inflow)
493 
494  if (trim(string_val1) .eq. "uniform") then
495  call json_get(params, 'case.fluid.inflow_condition.value', real_vec)
496  select type (bc_if => this%bc_inflow)
497  type is (inflow_t)
498  call bc_if%set_inflow(real_vec)
499  end select
500  else if (trim(string_val1) .eq. "blasius") then
501  select type (bc_if => this%bc_inflow)
502  type is (blasius_t)
503  call json_get(params, 'case.fluid.blasius.delta', real_val)
504  call json_get(params, 'case.fluid.blasius.approximation', &
505  string_val2)
506  call json_get(params, 'case.fluid.blasius.freestream_velocity', &
507  real_vec)
508 
509  call bc_if%set_params(real_vec, real_val, string_val2)
510 
511  end select
512  else if (trim(string_val1) .eq. "user") then
513  end if
514  end if
515 
516  !
517  ! Wall models
518  !
519 
520  call this%bc_wallmodel%init_base(this%c_Xh)
521  call this%bc_wallmodel%mark_zones_from_list(msh%labeled_zones,&
522  'wm', this%bc_labels)
523  call this%bc_wallmodel%finalize()
524 
525  if (this%bc_wallmodel%msk(0) .gt. 0) then
526  call json_extract_object(params, 'case.fluid.wall_modelling', wm_json)
527  call this%bc_wallmodel%init_wall_model_bc(wm_json, this%mu / this%rho)
528  else
529  call this%bc_wallmodel%shear_stress_t%init_shear_stress(this%c_Xh)
530  end if
531 
532  call bc_list_add(this%bclst_vel, this%bc_wallmodel%symmetry)
533  call bc_list_add(this%bclst_vel_neumann, this%bc_wallmodel)
534 
535  call this%bc_wall%init_base(this%c_Xh)
536  call this%bc_wall%mark_zone(msh%wall)
537  call this%bc_wall%mark_zones_from_list(msh%labeled_zones, &
538  'w', this%bc_labels)
539  call this%bc_wall%finalize()
540  call bc_list_add(this%bclst_vel, this%bc_wall)
541 
542  ! Setup field dirichlet bc for u-velocity
543  call this%user_field_bc_vel%bc_u%init_base(this%c_Xh)
544  call this%user_field_bc_vel%bc_u%mark_zones_from_list(msh%labeled_zones, &
545  'd_vel_u', this%bc_labels)
546  call this%user_field_bc_vel%bc_u%finalize()
547 
548  call mpi_allreduce(this%user_field_bc_vel%bc_u%msk(0), integer_val, 1, &
549  mpi_integer, mpi_sum, neko_comm, ierr)
550  if (integer_val .gt. 0) then
551  call this%user_field_bc_vel%bc_u%init_field('d_vel_u')
552  end if
553 
554  ! Setup field dirichlet bc for v-velocity
555  call this%user_field_bc_vel%bc_v%init_base(this%c_Xh)
556  call this%user_field_bc_vel%bc_v%mark_zones_from_list(msh%labeled_zones, &
557  'd_vel_v', this%bc_labels)
558  call this%user_field_bc_vel%bc_v%finalize()
559 
560  call mpi_allreduce(this%user_field_bc_vel%bc_v%msk(0), integer_val, 1, &
561  mpi_integer, mpi_sum, neko_comm, ierr)
562  if (integer_val .gt. 0) then
563  call this%user_field_bc_vel%bc_v%init_field('d_vel_v')
564  end if
565 
566  ! Setup field dirichlet bc for w-velocity
567  call this%user_field_bc_vel%bc_w%init_base(this%c_Xh)
568  call this%user_field_bc_vel%bc_w%mark_zones_from_list(msh%labeled_zones, &
569  'd_vel_w', this%bc_labels)
570  call this%user_field_bc_vel%bc_w%finalize()
571 
572  call mpi_allreduce(this%user_field_bc_vel%bc_w%msk(0), integer_val, 1, &
573  mpi_integer, mpi_sum, neko_comm, ierr)
574  if (integer_val .gt. 0) then
575  call this%user_field_bc_vel%bc_w%init_field('d_vel_w')
576  end if
577 
578  ! Setup our global field dirichlet bc
579  call this%user_field_bc_vel%init_base(this%c_Xh)
580  call this%user_field_bc_vel%mark_zones_from_list(msh%labeled_zones, &
581  'd_vel_u', this%bc_labels)
582  call this%user_field_bc_vel%mark_zones_from_list(msh%labeled_zones, &
583  'd_vel_v', this%bc_labels)
584  call this%user_field_bc_vel%mark_zones_from_list(msh%labeled_zones, &
585  'd_vel_w', this%bc_labels)
586  call this%user_field_bc_vel%finalize()
587 
588  ! Add the field bc to velocity bcs
589  call bc_list_add(this%bclst_vel, this%user_field_bc_vel)
590 
591  !
592  ! Associate our field dirichlet update to the user one.
593  !
594  this%user_field_bc_vel%update => user%user_dirichlet_update
595 
596  !
597  ! Initialize field list and bc list for user_dirichlet_update
598  !
599 
600  ! Note, some of these are potentially not initialized !
601  call this%user_field_bc_vel%field_list%init(4)
602  call this%user_field_bc_vel%field_list%assign_to_field(1, &
603  this%user_field_bc_vel%bc_u%field_bc)
604  call this%user_field_bc_vel%field_list%assign_to_field(2, &
605  this%user_field_bc_vel%bc_v%field_bc)
606  call this%user_field_bc_vel%field_list%assign_to_field(3, &
607  this%user_field_bc_vel%bc_w%field_bc)
608  call this%user_field_bc_vel%field_list%assign_to_field(4, &
609  this%user_field_bc_prs%field_bc)
610 
611  call bc_list_init(this%user_field_bc_vel%bc_list, size = 4)
612  ! Note, bc_list_add only adds if the bc is not empty
613  call bc_list_add(this%user_field_bc_vel%bc_list, &
614  this%user_field_bc_vel%bc_u)
615  call bc_list_add(this%user_field_bc_vel%bc_list, &
616  this%user_field_bc_vel%bc_v)
617  call bc_list_add(this%user_field_bc_vel%bc_list, &
618  this%user_field_bc_vel%bc_w)
619 
620  !
621  ! Check if we need to output boundary types to a separate field
622  call fluid_scheme_set_bc_type_output(this, params)
623 
624  !
625  ! Setup right-hand side fields.
626  !
627  allocate(this%f_x)
628  allocate(this%f_y)
629  allocate(this%f_z)
630  call this%f_x%init(this%dm_Xh, fld_name = "fluid_rhs_x")
631  call this%f_y%init(this%dm_Xh, fld_name = "fluid_rhs_y")
632  call this%f_z%init(this%dm_Xh, fld_name = "fluid_rhs_z")
633 
634  ! Initialize the source term
635  call this%source_term%init(this%f_x, this%f_y, this%f_z, this%c_Xh, user)
636  call this%source_term%add(params, 'case.fluid.source_terms')
637 
638  ! If case.output_boundary is true, set the values for the bc types in the
639  ! output of the field.
640  call this%set_bc_type_output(params)
641 
642  ! Initialize velocity solver
643  if (kspv_init) then
644  call neko_log%section("Velocity solver")
645  call json_get_or_default(params, &
646  'case.fluid.velocity_solver.max_iterations', &
647  integer_val, ksp_max_iter)
648  call json_get(params, 'case.fluid.velocity_solver.type', string_val1)
649  call json_get(params, 'case.fluid.velocity_solver.preconditioner', &
650  string_val2)
651  call json_get(params, 'case.fluid.velocity_solver.absolute_tolerance', &
652  real_val)
653  call json_get_or_default(params, &
654  'case.fluid.velocity_solver.monitor', &
655  logical_val, .false.)
656 
657  call neko_log%message('Type : ('// trim(string_val1) // &
658  ', ' // trim(string_val2) // ')')
659 
660  write(log_buf, '(A,ES13.6)') 'Abs tol :', real_val
661  call neko_log%message(log_buf)
662  call fluid_scheme_solver_factory(this%ksp_vel, this%dm_Xh%size(), &
663  string_val1, integer_val, real_val, logical_val)
664  call fluid_scheme_precon_factory(this%pc_vel, this%ksp_vel, &
665  this%c_Xh, this%dm_Xh, this%gs_Xh, this%bclst_vel, string_val2)
666  call neko_log%end_section()
667  end if
668 
669  ! Assign velocity fields
670  call neko_field_registry%add_field(this%dm_Xh, 'u')
671  call neko_field_registry%add_field(this%dm_Xh, 'v')
672  call neko_field_registry%add_field(this%dm_Xh, 'w')
673  this%u => neko_field_registry%get_field('u')
674  this%v => neko_field_registry%get_field('v')
675  this%w => neko_field_registry%get_field('w')
676 
677  !! Initialize time-lag fields
678  call this%ulag%init(this%u, 2)
679  call this%vlag%init(this%v, 2)
680  call this%wlag%init(this%w, 2)
681 
682 
683  end subroutine fluid_scheme_init_common
684 
686  subroutine fluid_scheme_init_all(this, msh, lx, params, kspv_init, &
687  kspp_init, scheme, user)
688  implicit none
689  class(fluid_scheme_t), target, intent(inout) :: this
690  type(mesh_t), target, intent(inout) :: msh
691  integer, intent(inout) :: lx
692  type(json_file), target, intent(inout) :: params
693  type(user_t), target, intent(in) :: user
694  logical :: kspv_init
695  logical :: kspp_init
696  character(len=*), intent(in) :: scheme
697  real(kind=rp) :: abs_tol
698  integer :: integer_val, ierr
699  logical :: logical_val
700  character(len=:), allocatable :: solver_type, precon_type
701  character(len=LOG_SIZE) :: log_buf
702  real(kind=rp) :: gjp_param_a, gjp_param_b
703 
704  call fluid_scheme_init_common(this, msh, lx, params, scheme, user, &
705  kspv_init)
706 
707  call neko_field_registry%add_field(this%dm_Xh, 'p')
708  this%p => neko_field_registry%get_field('p')
709 
710  !
711  ! Setup pressure boundary conditions
712  !
713  call bc_list_init(this%bclst_prs)
714  call this%bc_prs%init_base(this%c_Xh)
715  call this%bc_prs%mark_zones_from_list(msh%labeled_zones, &
716  'o', this%bc_labels)
717  call this%bc_prs%mark_zones_from_list(msh%labeled_zones, &
718  'on', this%bc_labels)
719 
720  ! Field dirichlet pressure bc
721  call this%user_field_bc_prs%init_base(this%c_Xh)
722  call this%user_field_bc_prs%mark_zones_from_list(msh%labeled_zones, &
723  'd_pres', this%bc_labels)
724  call this%user_field_bc_prs%finalize()
725  call mpi_allreduce(this%user_field_bc_prs%msk(0), integer_val, 1, &
726  mpi_integer, mpi_sum, neko_comm, ierr)
727 
728  if (integer_val .gt. 0) call this%user_field_bc_prs%init_field('d_pres')
729  call bc_list_add(this%bclst_prs, this%user_field_bc_prs)
730  call bc_list_add(this%user_field_bc_vel%bc_list, this%user_field_bc_prs)
731 
732  if (msh%outlet%size .gt. 0) then
733  call this%bc_prs%mark_zone(msh%outlet)
734  end if
735  if (msh%outlet_normal%size .gt. 0) then
736  call this%bc_prs%mark_zone(msh%outlet_normal)
737  end if
738 
739  call this%bc_prs%finalize()
740  call this%bc_prs%set_g(0.0_rp)
741  call bc_list_add(this%bclst_prs, this%bc_prs)
742  call this%bc_dong%init_base(this%c_Xh)
743  call this%bc_dong%mark_zones_from_list(msh%labeled_zones, &
744  'o+dong', this%bc_labels)
745  call this%bc_dong%mark_zones_from_list(msh%labeled_zones, &
746  'on+dong', this%bc_labels)
747  call this%bc_dong%finalize()
748 
749 
750  call this%bc_dong%init(this%c_Xh, params)
751 
752  call bc_list_add(this%bclst_prs, this%bc_dong)
753 
754  ! Pressure solver
755  if (kspp_init) then
756  call neko_log%section("Pressure solver")
757 
758  call json_get_or_default(params, &
759  'case.fluid.pressure_solver.max_iterations', &
760  integer_val, ksp_max_iter)
761  call json_get(params, 'case.fluid.pressure_solver.type', solver_type)
762  call json_get(params, 'case.fluid.pressure_solver.preconditioner', &
763  precon_type)
764  call json_get(params, 'case.fluid.pressure_solver.absolute_tolerance', &
765  abs_tol)
766  call json_get_or_default(params, &
767  'case.fluid.pressure_solver.monitor', &
768  logical_val, .false.)
769  call neko_log%message('Type : ('// trim(solver_type) // &
770  ', ' // trim(precon_type) // ')')
771  write(log_buf, '(A,ES13.6)') 'Abs tol :', abs_tol
772  call neko_log%message(log_buf)
773 
774  call fluid_scheme_solver_factory(this%ksp_prs, this%dm_Xh%size(), &
775  solver_type, integer_val, abs_tol, logical_val)
776  call fluid_scheme_precon_factory(this%pc_prs, this%ksp_prs, &
777  this%c_Xh, this%dm_Xh, this%gs_Xh, this%bclst_prs, precon_type)
778 
779  call neko_log%end_section()
780 
781  end if
782 
783  ! Initiate gradient jump penalty
784  call json_get_or_default(params, &
785  'case.fluid.gradient_jump_penalty.enabled',&
786  this%if_gradient_jump_penalty, .false.)
787 
788  if (this%if_gradient_jump_penalty .eqv. .true.) then
789  if ((this%dm_Xh%xh%lx - 1) .eq. 1) then
790  call json_get_or_default(params, &
791  'case.fluid.gradient_jump_penalty.tau',&
792  gjp_param_a, 0.02_rp)
793  gjp_param_b = 0.0_rp
794  else
795  call json_get_or_default(params, &
796  'case.fluid.gradient_jump_penalty.scaling_factor',&
797  gjp_param_a, 0.8_rp)
798  call json_get_or_default(params, &
799  'case.fluid.gradient_jump_penalty.scaling_exponent',&
800  gjp_param_b, 4.0_rp)
801  end if
802  call this%gradient_jump_penalty_u%init(params, this%dm_Xh, this%c_Xh, &
803  gjp_param_a, gjp_param_b)
804  call this%gradient_jump_penalty_v%init(params, this%dm_Xh, this%c_Xh, &
805  gjp_param_a, gjp_param_b)
806  call this%gradient_jump_penalty_w%init(params, this%dm_Xh, this%c_Xh, &
807  gjp_param_a, gjp_param_b)
808  end if
809 
810  call neko_log%end_section()
811 
812  end subroutine fluid_scheme_init_all
813 
815  subroutine fluid_scheme_free(this)
816  class(fluid_scheme_t), intent(inout) :: this
817 
818  call this%bdry%free()
819 
820  if (allocated(this%bc_inflow)) then
821  call this%bc_inflow%free()
822  end if
823 
824  call this%bc_wall%free()
825  call this%bc_sym%free()
826  call this%bc_sh%free()
827 
828  !
829  ! Free everything related to field_dirichlet BCs
830  !
831  call this%user_field_bc_prs%field_bc%free()
832  call this%user_field_bc_prs%free()
833  call this%user_field_bc_vel%bc_u%field_bc%free()
834  call this%user_field_bc_vel%bc_v%field_bc%free()
835  call this%user_field_bc_vel%bc_w%field_bc%free()
836  call this%user_field_bc_vel%free()
837 
838  call this%Xh%free()
839 
840  if (allocated(this%ksp_vel)) then
841  call krylov_solver_destroy(this%ksp_vel)
842  deallocate(this%ksp_vel)
843  end if
844 
845  if (allocated(this%ksp_prs)) then
846  call krylov_solver_destroy(this%ksp_prs)
847  deallocate(this%ksp_prs)
848  end if
849 
850  if (allocated(this%pc_vel)) then
851  call precon_destroy(this%pc_vel)
852  deallocate(this%pc_vel)
853  end if
854 
855  if (allocated(this%pc_prs)) then
856  call precon_destroy(this%pc_prs)
857  deallocate(this%pc_prs)
858  end if
859 
860  if (allocated(this%bc_labels)) then
861  deallocate(this%bc_labels)
862  end if
863 
864  call this%source_term%free()
865 
866  call this%gs_Xh%free()
867 
868  call this%c_Xh%free()
869 
870  call bc_list_free(this%bclst_vel)
871 
872  call this%scratch%free()
873 
874  nullify(this%params)
875 
876  nullify(this%u)
877  nullify(this%v)
878  nullify(this%w)
879  nullify(this%p)
880 
881  call this%ulag%free()
882  call this%vlag%free()
883  call this%wlag%free()
884 
885 
886  if (associated(this%f_x)) then
887  call this%f_x%free()
888  end if
889 
890  if (associated(this%f_y)) then
891  call this%f_y%free()
892  end if
893 
894  if (associated(this%f_z)) then
895  call this%f_z%free()
896  end if
897 
898  nullify(this%f_x)
899  nullify(this%f_y)
900  nullify(this%f_z)
901 
902  call this%mu_field%free()
903 
904  ! Free gradient jump penalty
905  if (this%if_gradient_jump_penalty .eqv. .true.) then
906  call this%gradient_jump_penalty_u%free()
907  call this%gradient_jump_penalty_v%free()
908  call this%gradient_jump_penalty_w%free()
909  end if
910 
911 
912  end subroutine fluid_scheme_free
913 
916  subroutine fluid_scheme_validate(this)
917  class(fluid_scheme_t), target, intent(inout) :: this
918  ! Variables for retrieving json parameters
919  logical :: logical_val
920 
921  if ( (.not. associated(this%u)) .or. &
922  (.not. associated(this%v)) .or. &
923  (.not. associated(this%w)) .or. &
924  (.not. associated(this%p))) then
925  call neko_error('Fields are not registered')
926  end if
927 
928  if ( (.not. allocated(this%u%x)) .or. &
929  (.not. allocated(this%v%x)) .or. &
930  (.not. allocated(this%w%x)) .or. &
931  (.not. allocated(this%p%x))) then
932  call neko_error('Fields are not allocated')
933  end if
934 
935  if (.not. allocated(this%ksp_vel)) then
936  call neko_error('No Krylov solver for velocity defined')
937  end if
938 
939  if (.not. allocated(this%ksp_prs)) then
940  call neko_error('No Krylov solver for pressure defined')
941  end if
942 
943  if (.not. associated(this%params)) then
944  call neko_error('No parameters defined')
945  end if
946 
947  if (allocated(this%bc_inflow)) then
948  select type (ip => this%bc_inflow)
949  type is (usr_inflow_t)
950  call ip%validate
951  end select
952  end if
953 
954  !
955  ! Setup checkpoint structure (if everything is fine)
956  !
957  call this%chkp%init(this%u, this%v, this%w, this%p)
958 
959  end subroutine fluid_scheme_validate
960 
965  subroutine fluid_scheme_bc_apply_vel(this, t, tstep)
966  class(fluid_scheme_t), intent(inout) :: this
967  real(kind=rp), intent(in) :: t
968  integer, intent(in) :: tstep
969 
970  call bc_list_apply_vector(this%bclst_vel, &
971  this%u%x, this%v%x, this%w%x, this%dm_Xh%size(), t, tstep)
972 
973  end subroutine fluid_scheme_bc_apply_vel
974 
977  subroutine fluid_scheme_bc_apply_prs(this, t, tstep)
978  class(fluid_scheme_t), intent(inout) :: this
979  real(kind=rp), intent(in) :: t
980  integer, intent(in) :: tstep
981 
982  call bc_list_apply_scalar(this%bclst_prs, this%p%x, &
983  this%p%dof%size(), t, tstep)
984 
985  end subroutine fluid_scheme_bc_apply_prs
986 
989  subroutine fluid_scheme_solver_factory(ksp, n, solver, &
990  max_iter, abstol, monitor)
991  class(ksp_t), allocatable, target, intent(inout) :: ksp
992  integer, intent(in), value :: n
993  character(len=*), intent(in) :: solver
994  integer, intent(in) :: max_iter
995  real(kind=rp), intent(in) :: abstol
996  logical, intent(in) :: monitor
997 
998  call krylov_solver_factory(ksp, n, solver, max_iter, abstol, &
999  monitor = monitor)
1000 
1001  end subroutine fluid_scheme_solver_factory
1002 
1004  subroutine fluid_scheme_precon_factory(pc, ksp, coef, dof, gs, bclst, &
1005  pctype)
1006  class(pc_t), allocatable, target, intent(inout) :: pc
1007  class(ksp_t), target, intent(inout) :: ksp
1008  type(coef_t), target, intent(inout) :: coef
1009  type(dofmap_t), target, intent(inout) :: dof
1010  type(gs_t), target, intent(inout) :: gs
1011  type(bc_list_t), target, intent(inout) :: bclst
1012  character(len=*) :: pctype
1013 
1014  call precon_factory(pc, pctype)
1015 
1016  select type (pcp => pc)
1017  type is (jacobi_t)
1018  call pcp%init(coef, dof, gs)
1019  type is (sx_jacobi_t)
1020  call pcp%init(coef, dof, gs)
1021  type is (device_jacobi_t)
1022  call pcp%init(coef, dof, gs)
1023  type is (hsmg_t)
1024  if (len_trim(pctype) .gt. 4) then
1025  if (index(pctype, '+') .eq. 5) then
1026  call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst, &
1027  trim(pctype(6:)))
1028  else
1029  call neko_error('Unknown coarse grid solver')
1030  end if
1031  else
1032  call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst)
1033  end if
1034  end select
1035 
1036  call ksp%set_pc(pc)
1037 
1038  end subroutine fluid_scheme_precon_factory
1039 
1041  subroutine fluid_scheme_set_usr_inflow(this, usr_eval)
1042  class(fluid_scheme_t), intent(inout) :: this
1043  procedure(usr_inflow_eval) :: usr_eval
1044 
1045  select type (bc_if => this%bc_inflow)
1046  type is (usr_inflow_t)
1047  call bc_if%set_eval(usr_eval)
1048  class default
1049  call neko_error("Not a user defined inflow condition")
1050  end select
1051  end subroutine fluid_scheme_set_usr_inflow
1052 
1054  function fluid_compute_cfl(this, dt) result(c)
1055  class(fluid_scheme_t), intent(in) :: this
1056  real(kind=rp), intent(in) :: dt
1057  real(kind=rp) :: c
1058 
1059  c = cfl(dt, this%u%x, this%v%x, this%w%x, &
1060  this%Xh, this%c_Xh, this%msh%nelv, this%msh%gdim)
1061 
1062  end function fluid_compute_cfl
1063 
1066  subroutine fluid_scheme_set_bc_type_output(this, params)
1067  class(fluid_scheme_t), target, intent(inout) :: this
1068  type(json_file), intent(inout) :: params
1069  type(dirichlet_t) :: bdry_mask
1070  logical :: found
1071 
1072  !
1073  ! Check if we need to output boundaries
1074  !
1075  call json_get_or_default(params, 'case.output_boundary', found, .false.)
1076 
1077  if (found) then
1078  call this%bdry%init(this%dm_Xh, 'bdry')
1079  this%bdry = 0.0_rp
1080 
1081  call bdry_mask%init_base(this%c_Xh)
1082  call bdry_mask%mark_zone(this%msh%wall)
1083  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1084  'w', this%bc_labels)
1085  call bdry_mask%finalize()
1086  call bdry_mask%set_g(1.0_rp)
1087  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1088  call bdry_mask%free()
1089 
1090  call bdry_mask%init_base(this%c_Xh)
1091  call bdry_mask%mark_zone(this%msh%inlet)
1092  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1093  'v', this%bc_labels)
1094  call bdry_mask%finalize()
1095  call bdry_mask%set_g(2.0_rp)
1096  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1097  call bdry_mask%free()
1098 
1099  call bdry_mask%init_base(this%c_Xh)
1100  call bdry_mask%mark_zone(this%msh%outlet)
1101  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1102  'o', this%bc_labels)
1103  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1104  'o+dong', this%bc_labels)
1105  call bdry_mask%finalize()
1106  call bdry_mask%set_g(3.0_rp)
1107  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1108  call bdry_mask%free()
1109 
1110  call bdry_mask%init_base(this%c_Xh)
1111  call bdry_mask%mark_zone(this%msh%sympln)
1112  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1113  'sym', this%bc_labels)
1114  call bdry_mask%finalize()
1115  call bdry_mask%set_g(4.0_rp)
1116  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1117  call bdry_mask%free()
1118 
1119  call bdry_mask%init_base(this%c_Xh)
1120  call bdry_mask%mark_zone(this%msh%outlet_normal)
1121  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1122  'on', this%bc_labels)
1123  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1124  'on+dong', this%bc_labels)
1125  call bdry_mask%finalize()
1126  call bdry_mask%set_g(5.0_rp)
1127  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1128  call bdry_mask%free()
1129 
1130  call bdry_mask%init_base(this%c_Xh)
1131  call bdry_mask%mark_zone(this%msh%periodic)
1132  call bdry_mask%finalize()
1133  call bdry_mask%set_g(6.0_rp)
1134  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1135  call bdry_mask%free()
1136 
1137  call bdry_mask%init_base(this%c_Xh)
1138  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1139  'd_vel_u', this%bc_labels)
1140  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1141  'd_vel_v', this%bc_labels)
1142  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1143  'd_vel_w', this%bc_labels)
1144  call bdry_mask%finalize()
1145  call bdry_mask%set_g(7.0_rp)
1146  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1147  call bdry_mask%free()
1148 
1149  call bdry_mask%init_base(this%c_Xh)
1150  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1151  'd_pres', this%bc_labels)
1152  call bdry_mask%finalize()
1153  call bdry_mask%set_g(8.0_rp)
1154  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1155  call bdry_mask%free()
1156 
1157  call bdry_mask%init_base(this%c_Xh)
1158  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1159  'sh', this%bc_labels)
1160  call bdry_mask%finalize()
1161  call bdry_mask%set_g(9.0_rp)
1162  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1163  call bdry_mask%free()
1164 
1165  call bdry_mask%init_base(this%c_Xh)
1166  call bdry_mask%mark_zones_from_list(this%msh%labeled_zones, &
1167  'wm', this%bc_labels)
1168  call bdry_mask%finalize()
1169  call bdry_mask%set_g(10.0_rp)
1170  call bdry_mask%apply_scalar(this%bdry%x, this%dm_Xh%size())
1171  call bdry_mask%free()
1172 
1173  end if
1174 
1175  end subroutine fluid_scheme_set_bc_type_output
1176 
1178  subroutine fluid_scheme_update_material_properties(this)
1179  class(fluid_scheme_t), intent(inout) :: this
1180  type(field_t), pointer :: nut
1181  integer :: n
1182 
1183  if (this%variable_material_properties) then
1184  nut => neko_field_registry%get_field(this%nut_field_name)
1185  n = nut%size()
1186 
1187  if (neko_bcknd_device .eq. 1) then
1188  call device_cfill(this%mu_field%x_d, this%mu, n)
1189  call device_add2s2(this%mu_field%x_d, nut%x_d, this%rho, n)
1190  else
1191  call cfill(this%mu_field%x, this%mu, n)
1192  call add2s2(this%mu_field%x, nut%x, this%rho, n)
1193  end if
1194  end if
1195 
1196  end subroutine fluid_scheme_update_material_properties
1197 
1201  subroutine fluid_scheme_set_material_properties(this, params, user)
1202  class(fluid_scheme_t), intent(inout) :: this
1203  type(json_file), intent(inout) :: params
1204  type(user_t), target, intent(in) :: user
1205  character(len=LOG_SIZE) :: log_buf
1206  ! A local pointer that is needed to make Intel happy
1207  procedure(user_material_properties), pointer :: dummy_mp_ptr
1208  logical :: nondimensional
1209  real(kind=rp) :: dummy_lambda, dummy_cp
1210 
1211  dummy_mp_ptr => dummy_user_material_properties
1212 
1213  if (.not. associated(user%material_properties, dummy_mp_ptr)) then
1214 
1215  write(log_buf, '(A)') "Material properties must be set in the user&
1216  & file!"
1217  call neko_log%message(log_buf)
1218  call user%material_properties(0.0_rp, 0, this%rho, this%mu, &
1219  dummy_cp, dummy_lambda, params)
1220  else
1221  ! Incorrect user input
1222  if (params%valid_path('case.fluid.Re') .and. &
1223  (params%valid_path('case.fluid.mu') .or. &
1224  params%valid_path('case.fluid.rho'))) then
1225  call neko_error("To set the material properties for the fluid,&
1226  & either provide Re OR mu and rho in the case file.")
1227 
1228  ! Non-dimensional case
1229  else if (params%valid_path('case.fluid.Re')) then
1230 
1231  write(log_buf, '(A)') 'Non-dimensional fluid material properties &
1232  & input.'
1233  call neko_log%message(log_buf, lvl = neko_log_verbose)
1234  write(log_buf, '(A)') 'Density will be set to 1, dynamic viscosity to&
1235  & 1/Re.'
1236  call neko_log%message(log_buf, lvl = neko_log_verbose)
1237 
1238  ! Read Re into mu for further manipulation.
1239  call json_get(params, 'case.fluid.Re', this%mu)
1240  write(log_buf, '(A)') 'Read non-dimensional material properties'
1241  call neko_log%message(log_buf)
1242  write(log_buf, '(A,ES13.6)') 'Re :', this%mu
1243  call neko_log%message(log_buf)
1244 
1245  ! Set rho to 1 since the setup is non-dimensional.
1246  this%rho = 1.0_rp
1247  ! Invert the Re to get viscosity.
1248  this%mu = 1.0_rp/this%mu
1249  ! Dimensional case
1250  else
1251  call json_get(params, 'case.fluid.mu', this%mu)
1252  call json_get(params, 'case.fluid.rho', this%rho)
1253  end if
1254 
1255  end if
1256  end subroutine fluid_scheme_set_material_properties
1257 
1258 end module fluid_scheme
Abstract interface to dealocate a fluid formulation.
Abstract interface to initialize a fluid formulation.
Abstract interface to restart a fluid scheme.
Abstract interface to compute a time-step.
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Definition: json_utils.f90:54
Retrieves a parameter by name or throws an error.
Definition: json_utils.f90:45
Abstract interface for setting material properties.
Definition: user_intf.f90:147
Abstract interface defining a user defined inflow condition (pointwise)
Definition: usr_inflow.f90:80
Defines a boundary condition.
Definition: bc.f90:34
subroutine, public bc_list_add(bclst, bc)
Add a condition to a list of boundary conditions.
Definition: bc.f90:505
subroutine, public bc_list_apply_vector(bclst, x, y, z, n, t, tstep)
Apply a list of boundary conditions to a vector field.
Definition: bc.f90:587
subroutine, public bc_list_init(bclst, size)
Constructor for a list of boundary conditions.
Definition: bc.f90:464
subroutine, public bc_list_free(bclst)
Destructor for a list of boundary conditions.
Definition: bc.f90:491
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
Definition: bc.f90:530
Defines a Blasius profile dirichlet condition.
Definition: blasius.f90:34
Defines a checkpoint.
Definition: checkpoint.f90:34
Coefficients.
Definition: coef.f90:34
Definition: comm.F90:1
type(mpi_comm) neko_comm
MPI communicator.
Definition: comm.F90:16
Jacobi preconditioner accelerator backend.
subroutine, public device_add2s2(a_d, b_d, c1, n)
Vector addition with scalar multiplication (multiplication on first argument)
subroutine, public device_cfill(a_d, c, n)
Set all elements to a constant c .
Defines a dirichlet boundary condition.
Definition: dirichlet.f90:34
Defines a mapping of the degrees of freedom.
Definition: dofmap.f90:35
Defines a dong outflow condition.
Defines inflow dirichlet conditions.
Defines inflow dirichlet conditions.
subroutine, public field_cfill(a, c, n)
Set all elements to a constant c .
Definition: field_math.f90:186
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Stores a series fields.
Defines a field.
Definition: field.f90:34
Fluid formulations.
subroutine fluid_scheme_set_material_properties(this, params, user)
Sets rho and mu.
subroutine fluid_scheme_free(this)
Deallocate a fluid formulation.
subroutine fluid_scheme_update_material_properties(this)
Update the values of mu_field if necessary.
subroutine fluid_scheme_set_bc_type_output(this, params)
Set boundary types for the diagnostic output.
real(kind=rp) function fluid_compute_cfl(this, dt)
Compute CFL.
subroutine fluid_scheme_init_all(this, msh, lx, params, kspv_init, kspp_init, scheme, user)
Initialize all components of the current scheme.
subroutine fluid_scheme_set_usr_inflow(this, usr_eval)
Initialize a user defined inflow condition.
subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, kspv_init)
Initialise a fluid scheme.
subroutine fluid_scheme_bc_apply_vel(this, t, tstep)
Apply all boundary conditions defined for velocity Here we perform additional gs operations to take c...
subroutine fluid_scheme_validate(this)
Validate that all fields, solvers etc necessary for performing time-stepping are defined.
subroutine fluid_scheme_bc_apply_prs(this, t, tstep)
Apply all boundary conditions defined for pressure.
Implements the fluid_source_term_t type.
Computes various statistics for the fluid fields. We use the Reynolds decomposition for a field u = ...
Definition: fluid_stats.f90:36
Gather-scatter.
Implements gradient_jump_penalty_t.
Krylov preconditioner.
Definition: pc_hsmg.f90:61
Defines inflow dirichlet conditions.
Definition: inflow.f90:34
Jacobi preconditioner.
Definition: pc_jacobi.f90:34
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
subroutine, public json_extract_object(json, name, object)
Extract object as a separate JSON dictionary.
Definition: json_utils.f90:430
Implements the base abstract type for Krylov solvers plus helper types.
Definition: krylov.f90:34
integer, parameter, public ksp_max_iter
Maximum number of iters.
Definition: krylov.f90:51
Logging routines.
Definition: log.f90:34
integer, parameter, public neko_log_verbose
Verbose log level.
Definition: log.f90:71
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
integer, parameter, public log_size
Definition: log.f90:42
Definition: math.f90:60
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition: math.f90:348
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
Definition: math.f90:673
Defines a mean flow field.
Definition: mean_flow.f90:34
Defines a mean squared flow field.
Defines a mesh.
Definition: mesh.f90:34
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
Definition: mesh.f90:56
integer, parameter, public neko_msh_max_zlbl_len
Max length of a zone label.
Definition: mesh.f90:58
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
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(dt, u, v, w, Xh, coef, nelv, gdim)
Definition: operators.f90:394
Krylov preconditioner.
Definition: precon.f90:34
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
Defines a shear stress boundary condition for a vector field. Maintainer: Timofey Mukha.
Implements the source_term_t type and a wrapper source_term_wrapper_t.
Definition: source_term.f90:34
Defines a function space.
Definition: space.f90:34
integer, parameter, public gll
Definition: space.f90:48
Defines a container for all statistics.
Definition: statistics.f90:34
Jacobi preconditioner SX-Aurora backend.
Mixed Dirichlet-Neumann axis aligned symmetry plane.
Definition: symmetry.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
subroutine, public dummy_user_material_properties(t, tstep, rho, mu, cp, lambda, params)
Definition: user_intf.f90:441
Defines inflow dirichlet conditions.
Definition: usr_inflow.f90:34
Utilities.
Definition: utils.f90:35
subroutine, public neko_warning(warning_msg)
Reports a warning to standard output.
Definition: utils.f90:245
Defines the wall_model_bc_t type. Maintainer: Timofey Mukha.
Defines wall boundary conditions.
Definition: wall.f90:34
A list of boundary conditions.
Definition: bc.f90:104
Base type for a boundary condition.
Definition: bc.f90:51
Blasius profile for inlet (vector valued)
Definition: blasius.f90:48
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
Defines a jacobi preconditioner.
Generic Dirichlet boundary condition on .
Definition: dirichlet.f90:44
Dong outflow condition Follows "A Convective-like Energy-Stable Open Boundary Condition for Simulati...
User defined dirichlet condition, for which the user can work with an entire field....
Extension of the user defined dirichlet condition field_dirichlet
Base type of all fluid formulations.
Wrapper contaning and executing the fluid source terms.
Dirichlet condition for inlet (vector valued)
Definition: inflow.f90:43
Defines a jacobi preconditioner.
Definition: pc_jacobi.f90:45
Base abstract type for a canonical Krylov method, solving .
Definition: krylov.f90:66
Defines a canonical Krylov preconditioner.
Definition: precon.f90:40
A shear stress boundary condition.
The function space for the SEM solution fields.
Definition: space.f90:62
Defines a jacobi preconditioner for SX-Aurora.
Mixed Dirichlet-Neumann symmetry plane condition.
Definition: symmetry.f90:46
Implements the logic to compute the time coefficients for the advection and diffusion operators in a ...
Provides a tool to set time step dt.
User defined dirichlet condition for inlet (vector valued)
Definition: usr_inflow.f90:46
No-slip Wall boundary condition.
Definition: wall.f90:43
A shear stress boundary condition, computing the stress values using a wall model.