Neko  0.9.0
A portable framework for high-order spectral element flow simulations
scalar_scheme.f90
Go to the documentation of this file.
1 ! Copyright (c) 2022-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 !
34 
35 ! todo: module name
37  use gather_scatter, only : gs_t
38  use checkpoint, only : chkp_t
39  use num_types, only: rp
40  use field, only : field_t
41  use field_list, only: field_list_t
42  use space, only : space_t
43  use dofmap, only : dofmap_t
44  use krylov, only : ksp_t, krylov_solver_factory, krylov_solver_destroy, &
46  use coefs, only : coef_t
47  use dirichlet, only : dirichlet_t
48  use neumann, only : neumann_t
49  use jacobi, only : jacobi_t
50  use device_jacobi, only : device_jacobi_t
51  use sx_jacobi, only : sx_jacobi_t
52  use hsmg, only : hsmg_t
53  use precon, only : pc_t, precon_factory, precon_destroy
54  use bc, only : bc_t, bc_list_t, bc_list_free, bc_list_init, &
58  use facet_zone, only : facet_zone_t
64  use json_module, only : json_file
67  use utils, only : neko_error
68  use comm, only: neko_comm, mpi_integer, mpi_sum
70  use math, only : cfill, add2s2
72  use neko_config, only : neko_bcknd_device
73  use field_series, only : field_series_t
76  implicit none
77 
79  type, abstract :: scalar_scheme_t
81  type(field_t), pointer :: u
83  type(field_t), pointer :: v
85  type(field_t), pointer :: w
87  type(field_t), pointer :: s
89  type(field_series_t) :: slag
91  type(space_t), pointer :: xh
93  type(dofmap_t), pointer :: dm_xh
95  type(gs_t), pointer :: gs_xh
97  type(coef_t), pointer :: c_xh
99  type(field_t), pointer :: f_xh => null()
103  class(ksp_t), allocatable :: ksp
105  integer :: ksp_maxiter
107  integer :: projection_dim
108 
109  integer :: projection_activ_step
111  class(pc_t), allocatable :: pc
113  type(dirichlet_t) :: dir_bcs(neko_msh_max_zlbls)
115  type(field_dirichlet_t) :: field_dir_bc
117  type(bc_list_t) :: field_dirichlet_bcs
119  type(neumann_t) :: neumann_bcs(neko_msh_max_zlbls)
121  type(usr_scalar_t) :: user_bc
123  integer :: n_dir_bcs = 0
125  integer :: n_neumann_bcs = 0
127  type(bc_list_t) :: bclst_dirichlet
129  type(bc_list_t) :: bclst_neumann
131  type(json_file), pointer :: params
133  type(mesh_t), pointer :: msh => null()
135  type(chkp_t) :: chkp
137  real(kind=rp) :: lambda
139  type(field_t) :: lambda_field
141  character(len=:), allocatable :: nut_field_name
143  real(kind=rp) :: rho
145  real(kind=rp) :: cp
147  real(kind=rp) :: pr_turb
149  logical :: variable_material_properties = .false.
151  character(len=NEKO_MSH_MAX_ZLBL_LEN), allocatable :: bc_labels(:)
153  logical :: if_gradient_jump_penalty
155  contains
157  procedure, pass(this) :: scheme_init => scalar_scheme_init
159  procedure, pass(this) :: scheme_free => scalar_scheme_free
161  procedure, pass(this) :: validate => scalar_scheme_validate
163  procedure, pass(this) :: set_user_bc => scalar_scheme_set_user_bc
165  procedure, pass(this) :: set_material_properties => &
168  procedure, pass(this) :: update_material_properties => &
171  procedure(scalar_scheme_init_intrf), pass(this), deferred :: init
173  procedure(scalar_scheme_free_intrf), pass(this), deferred :: free
175  procedure(scalar_scheme_step_intrf), pass(this), deferred :: step
177  procedure(scalar_scheme_restart_intrf), pass(this), deferred :: restart
178  end type scalar_scheme_t
179 
181  abstract interface
182  subroutine scalar_scheme_init_intrf(this, msh, coef, gs, params, user, &
183  ulag, vlag, wlag, time_scheme, rho)
184  import scalar_scheme_t
185  import json_file
186  import coef_t
187  import gs_t
188  import mesh_t
189  import user_t
190  import field_series_t
192  import rp
193  class(scalar_scheme_t), target, intent(inout) :: this
194  type(mesh_t), target, intent(inout) :: msh
195  type(coef_t), target, intent(inout) :: coef
196  type(gs_t), target, intent(inout) :: gs
197  type(json_file), target, intent(inout) :: params
198  type(user_t), target, intent(in) :: user
199  type(field_series_t), target, intent(in) :: ulag, vlag, wlag
200  type(time_scheme_controller_t), target, intent(in) :: time_scheme
201  real(kind=rp), intent(in) :: rho
202  end subroutine scalar_scheme_init_intrf
203  end interface
204 
206  abstract interface
207  subroutine scalar_scheme_restart_intrf(this, dtlag, tlag)
208  import scalar_scheme_t
209  import chkp_t
210  import rp
211  class(scalar_scheme_t), target, intent(inout) :: this
212  real(kind=rp) :: dtlag(10), tlag(10)
213  end subroutine scalar_scheme_restart_intrf
214  end interface
215 
216 
218  abstract interface
219  subroutine scalar_scheme_free_intrf(this)
220  import scalar_scheme_t
221  class(scalar_scheme_t), intent(inout) :: this
222  end subroutine scalar_scheme_free_intrf
223  end interface
224 
226  abstract interface
227  subroutine scalar_scheme_step_intrf(this, t, tstep, dt, ext_bdf, &
228  dt_controller)
229  import scalar_scheme_t
232  import rp
233  class(scalar_scheme_t), intent(inout) :: this
234  real(kind=rp), intent(inout) :: t
235  integer, intent(inout) :: tstep
236  real(kind=rp), intent(in) :: dt
237  type(time_scheme_controller_t), intent(inout) :: ext_bdf
238  type(time_step_controller_t), intent(in) :: dt_controller
239  end subroutine scalar_scheme_step_intrf
240  end interface
241 
242 contains
243 
248  subroutine scalar_scheme_add_bcs(this, zones, bc_labels)
249  class(scalar_scheme_t), intent(inout) :: this
250  type(facet_zone_t), intent(inout) :: zones(NEKO_MSH_MAX_ZLBLS)
251  character(len=NEKO_MSH_MAX_ZLBL_LEN), intent(in) :: bc_labels(:)
252  character(len=NEKO_MSH_MAX_ZLBL_LEN) :: bc_label
253  integer :: i
254  real(kind=rp) :: dir_value, flux
255  logical :: bc_exists
256 
257  do i = 1, size(bc_labels)
258  bc_label = trim(bc_labels(i))
259  if (bc_label(1:2) .eq. 'd=') then
260 ! The idea of this commented piece of code is to merge bcs with the same
261 ! Dirichlet value into 1 so that one has less kernel launches. Currently
262 ! segfaults, needs investigation.
263 ! bc_exists = .false.
264 ! bc_idx = 0
265 ! do j = 1, i-1
266 ! if (bc_label .eq. bc_labels(j)) then
267 ! bc_exists = .true.
268 ! bc_idx = j
269 ! end if
270 ! end do
271 
272 ! if (bc_exists) then
273 ! call this%dir_bcs(j)%mark_zone(zones(i))
274 ! else
275  this%n_dir_bcs = this%n_dir_bcs + 1
276  call this%dir_bcs(this%n_dir_bcs)%init_base(this%c_Xh)
277  call this%dir_bcs(this%n_dir_bcs)%mark_zone(zones(i))
278  read(bc_label(3:), *) dir_value
279  call this%dir_bcs(this%n_dir_bcs)%set_g(dir_value)
280  call this%dir_bcs(this%n_dir_bcs)%finalize()
281  end if
282 
283  if (bc_label(1:2) .eq. 'n=') then
284  this%n_neumann_bcs = this%n_neumann_bcs + 1
285  call this%neumann_bcs(this%n_neumann_bcs)%init_base(this%c_Xh)
286  call this%neumann_bcs(this%n_neumann_bcs)%mark_zone(zones(i))
287  read(bc_label(3:), *) flux
288  call this%neumann_bcs(this%n_neumann_bcs)%finalize_neumann(flux)
289  end if
290 
292  if (bc_label(1:4) .eq. 'user') then
293  call this%user_bc%mark_zone(zones(i))
294  end if
295 
296  end do
297 
298  do i = 1, this%n_dir_bcs
299  call bc_list_add(this%bclst_dirichlet, this%dir_bcs(i))
300  end do
301 
302  ! Create list with just Neumann bcs
303  call bc_list_init(this%bclst_neumann, this%n_neumann_bcs)
304  do i = 1, this%n_neumann_bcs
305  call bc_list_add(this%bclst_neumann, this%neumann_bcs(i))
306  end do
307 
308  end subroutine scalar_scheme_add_bcs
309 
318  subroutine scalar_scheme_init(this, msh, c_Xh, gs_Xh, params, scheme, user, &
319  rho)
320  class(scalar_scheme_t), target, intent(inout) :: this
321  type(mesh_t), target, intent(inout) :: msh
322  type(coef_t), target, intent(inout) :: c_Xh
323  type(gs_t), target, intent(inout) :: gs_Xh
324  type(json_file), target, intent(inout) :: params
325  character(len=*), intent(in) :: scheme
326  type(user_t), target, intent(in) :: user
327  real(kind=rp), intent(in) :: rho
328  ! IO buffer for log output
329  character(len=LOG_SIZE) :: log_buf
330  ! Variables for retrieving json parameters
331  logical :: logical_val
332  real(kind=rp) :: real_val, solver_abstol
333  integer :: integer_val, ierr
334  character(len=:), allocatable :: solver_type, solver_precon
335  real(kind=rp) :: gjp_param_a, gjp_param_b
336 
337  this%u => neko_field_registry%get_field('u')
338  this%v => neko_field_registry%get_field('v')
339  this%w => neko_field_registry%get_field('w')
340 
341  call neko_log%section('Scalar')
342  call json_get(params, 'case.fluid.velocity_solver.type', solver_type)
343  call json_get(params, 'case.fluid.velocity_solver.preconditioner', &
344  solver_precon)
345  call json_get(params, 'case.fluid.velocity_solver.absolute_tolerance', &
346  solver_abstol)
347 
348  call json_get_or_default(params, &
349  'case.fluid.velocity_solver.projection_space_size', &
350  this%projection_dim, 20)
351  call json_get_or_default(params, &
352  'case.fluid.velocity_solver.projection_hold_steps', &
353  this%projection_activ_step, 5)
354 
355 
356  write(log_buf, '(A, A)') 'Type : ', trim(scheme)
357  call neko_log%message(log_buf)
358  call neko_log%message('Ksp scalar : ('// trim(solver_type) // &
359  ', ' // trim(solver_precon) // ')')
360  write(log_buf, '(A,ES13.6)') ' `-abs tol :', solver_abstol
361  call neko_log%message(log_buf)
362 
363  this%Xh => this%u%Xh
364  this%dm_Xh => this%u%dof
365  this%params => params
366  this%msh => msh
367  if (.not. neko_field_registry%field_exists('s')) then
368  call neko_field_registry%add_field(this%dm_Xh, 's')
369  end if
370  this%s => neko_field_registry%get_field('s')
371 
372  call this%slag%init(this%s, 2)
373 
374  this%gs_Xh => gs_xh
375  this%c_Xh => c_xh
376 
377  !
378  ! Material properties
379  !
380  this%rho = rho
381  call this%set_material_properties(params, user)
382 
383  write(log_buf, '(A,ES13.6)') 'rho :', this%rho
384  call neko_log%message(log_buf)
385  write(log_buf, '(A,ES13.6)') 'lambda :', this%lambda
386  call neko_log%message(log_buf)
387  write(log_buf, '(A,ES13.6)') 'cp :', this%cp
388  call neko_log%message(log_buf)
389 
390  !
391  ! Turbulence modelling and variable material properties
392  !
393  if (params%valid_path('case.scalar.nut_field')) then
394  call json_get(params, 'case.scalar.Pr_t', this%pr_turb)
395  call json_get(params, 'case.scalar.nut_field', this%nut_field_name)
396  this%variable_material_properties = .true.
397  else
398  this%nut_field_name = ""
399  end if
400 
401  ! Fill lambda field with the physical value
402  call this%lambda_field%init(this%dm_Xh, "lambda")
403  if (neko_bcknd_device .eq. 1) then
404  call device_cfill(this%lambda_field%x_d, this%lambda, &
405  this%lambda_field%size())
406  else
407  call cfill(this%lambda_field%x, this%lambda, this%lambda_field%size())
408  end if
409 
410  !
411  ! Setup scalar boundary conditions
412  !
413  call bc_list_init(this%bclst_dirichlet)
414  call this%user_bc%init_base(this%c_Xh)
415 
416  ! Read boundary types from the case file
417  allocate(this%bc_labels(neko_msh_max_zlbls))
418 
419  ! A filler value
420  this%bc_labels = "not"
421 
422  if (params%valid_path('case.scalar.boundary_types')) then
423  call json_get(params, 'case.scalar.boundary_types', this%bc_labels, &
424  'not')
425  end if
426 
427  !
428  ! Setup right-hand side field.
429  !
430  allocate(this%f_Xh)
431  call this%f_Xh%init(this%dm_Xh, fld_name = "scalar_rhs")
432 
433  ! Initialize the source term
434  call this%source_term%init(this%f_Xh, this%c_Xh, user)
435  call this%source_term%add(params, 'case.scalar.source_terms')
436 
437  call scalar_scheme_add_bcs(this, msh%labeled_zones, this%bc_labels)
438 
439  ! Mark BC zones
440  call this%user_bc%mark_zone(msh%wall)
441  call this%user_bc%mark_zone(msh%inlet)
442  call this%user_bc%mark_zone(msh%outlet)
443  call this%user_bc%mark_zone(msh%outlet_normal)
444  call this%user_bc%mark_zone(msh%sympln)
445  call this%user_bc%finalize()
446  if (this%user_bc%msk(0) .gt. 0) call bc_list_add(this%bclst_dirichlet, &
447  this%user_bc)
448 
449  ! Add field dirichlet BCs
450  call this%field_dir_bc%init_base(this%c_Xh)
451  call this%field_dir_bc%mark_zones_from_list(msh%labeled_zones, &
452  'd_s', this%bc_labels)
453  call this%field_dir_bc%finalize()
454  call mpi_allreduce(this%field_dir_bc%msk(0), integer_val, 1, &
455  mpi_integer, mpi_sum, neko_comm, ierr)
456  if (integer_val .gt. 0) call this%field_dir_bc%init_field('d_s')
457 
458  call bc_list_add(this%bclst_dirichlet, this%field_dir_bc)
459 
460  !
461  ! Associate our field dirichlet update to the user one.
462  !
463  this%field_dir_bc%update => user%user_dirichlet_update
464 
465  call bc_list_init(this%field_dirichlet_bcs, size = 1)
466  call bc_list_add(this%field_dirichlet_bcs, this%field_dir_bc)
467 
468 
469  ! todo parameter file ksp tol should be added
470  call json_get_or_default(params, &
471  'case.fluid.velocity_solver.max_iterations', &
472  integer_val, ksp_max_iter)
473  call json_get_or_default(params, &
474  'case.fluid.velocity_solver.monitor', &
475  logical_val, .false.)
476  call scalar_scheme_solver_factory(this%ksp, this%dm_Xh%size(), &
477  solver_type, integer_val, solver_abstol, logical_val)
478  call scalar_scheme_precon_factory(this%pc, this%ksp, &
479  this%c_Xh, this%dm_Xh, this%gs_Xh, &
480  this%bclst_dirichlet, solver_precon)
481 
482  ! Initiate gradient jump penalty
483  call json_get_or_default(params, &
484  'case.scalar.gradient_jump_penalty.enabled',&
485  this%if_gradient_jump_penalty, .false.)
486 
487  if (this%if_gradient_jump_penalty .eqv. .true.) then
488  if ((this%dm_Xh%xh%lx - 1) .eq. 1) then
489  call json_get_or_default(params, &
490  'case.scalar.gradient_jump_penalty.tau',&
491  gjp_param_a, 0.02_rp)
492  gjp_param_b = 0.0_rp
493  else
494  call json_get_or_default(params, &
495  'case.scalar.gradient_jump_penalty.scaling_factor',&
496  gjp_param_a, 0.8_rp)
497  call json_get_or_default(params, &
498  'case.scalar.gradient_jump_penalty.scaling_exponent',&
499  gjp_param_b, 4.0_rp)
500  end if
501  call this%gradient_jump_penalty%init(params, this%dm_Xh, this%c_Xh, &
502  gjp_param_a, gjp_param_b)
503  end if
504 
505  call neko_log%end_section()
506 
507  end subroutine scalar_scheme_init
508 
509 
511  subroutine scalar_scheme_free(this)
512  class(scalar_scheme_t), intent(inout) :: this
513 
514  nullify(this%Xh)
515  nullify(this%dm_Xh)
516  nullify(this%gs_Xh)
517  nullify(this%c_Xh)
518  nullify(this%params)
519 
520  if (allocated(this%ksp)) then
521  call krylov_solver_destroy(this%ksp)
522  deallocate(this%ksp)
523  end if
524 
525  if (allocated(this%pc)) then
526  call precon_destroy(this%pc)
527  deallocate(this%pc)
528  end if
529 
530  if (allocated(this%bc_labels)) then
531  deallocate(this%bc_labels)
532  end if
533 
534  call this%source_term%free()
535 
536  call bc_list_free(this%bclst_dirichlet)
537  call bc_list_free(this%bclst_neumann)
538 
539  call this%lambda_field%free()
540  call this%slag%free()
541 
542  ! Free everything related to field dirichlet BCs
543  call bc_list_free(this%field_dirichlet_bcs)
544  call this%field_dir_bc%free()
545 
546  ! Free gradient jump penalty
547  if (this%if_gradient_jump_penalty .eqv. .true.) then
548  call this%gradient_jump_penalty%free()
549  end if
550 
551  end subroutine scalar_scheme_free
552 
555  subroutine scalar_scheme_validate(this)
556  class(scalar_scheme_t), target, intent(inout) :: this
557 
558  if ( (.not. allocated(this%u%x)) .or. &
559  (.not. allocated(this%v%x)) .or. &
560  (.not. allocated(this%w%x)) .or. &
561  (.not. allocated(this%s%x))) then
562  call neko_error('Fields are not allocated')
563  end if
564 
565  if (.not. allocated(this%ksp)) then
566  call neko_error('No Krylov solver for velocity defined')
567  end if
568 
569  if (.not. associated(this%Xh)) then
570  call neko_error('No function space defined')
571  end if
572 
573  if (.not. associated(this%dm_Xh)) then
574  call neko_error('No dofmap defined')
575  end if
576 
577  if (.not. associated(this%c_Xh)) then
578  call neko_error('No coefficients defined')
579  end if
580 
581  if (.not. associated(this%f_Xh)) then
582  call neko_error('No rhs allocated')
583  end if
584 
585  if (.not. associated(this%params)) then
586  call neko_error('No parameters defined')
587  end if
588 
589  !
590  ! Setup checkpoint structure (if everything is fine)
591  !
592 ! @todo no io for now
593 ! call this%chkp%init(this%u, this%v, this%w, this%p)
594 
595  end subroutine scalar_scheme_validate
596 
599  subroutine scalar_scheme_solver_factory(ksp, n, solver, max_iter, &
600  abstol, monitor)
601  class(ksp_t), allocatable, target, intent(inout) :: ksp
602  integer, intent(in), value :: n
603  integer, intent(in) :: max_iter
604  character(len=*), intent(in) :: solver
605  real(kind=rp) :: abstol
606  logical, intent(in) :: monitor
607 
608  call krylov_solver_factory(ksp, n, solver, max_iter, &
609  abstol, monitor = monitor)
610 
611  end subroutine scalar_scheme_solver_factory
612 
614  subroutine scalar_scheme_precon_factory(pc, ksp, coef, dof, gs, bclst, &
615  pctype)
616  class(pc_t), allocatable, target, intent(inout) :: pc
617  class(ksp_t), target, intent(inout) :: ksp
618  type(coef_t), target, intent(inout) :: coef
619  type(dofmap_t), target, intent(inout) :: dof
620  type(gs_t), target, intent(inout) :: gs
621  type(bc_list_t), target, intent(inout) :: bclst
622  character(len=*) :: pctype
623 
624  call precon_factory(pc, pctype)
625 
626  select type (pcp => pc)
627  type is (jacobi_t)
628  call pcp%init(coef, dof, gs)
629  type is (sx_jacobi_t)
630  call pcp%init(coef, dof, gs)
631  type is (device_jacobi_t)
632  call pcp%init(coef, dof, gs)
633  type is (hsmg_t)
634  if (len_trim(pctype) .gt. 4) then
635  if (index(pctype, '+') .eq. 5) then
636  call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst, &
637  trim(pctype(6:)))
638  else
639  call neko_error('Unknown coarse grid solver')
640  end if
641  else
642  call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst)
643  end if
644  end select
645 
646  call ksp%set_pc(pc)
647 
648  end subroutine scalar_scheme_precon_factory
649 
652  subroutine scalar_scheme_set_user_bc(this, usr_eval)
653  class(scalar_scheme_t), intent(inout) :: this
654  procedure(usr_scalar_bc_eval) :: usr_eval
655 
656  call this%user_bc%set_eval(usr_eval)
657 
658  end subroutine scalar_scheme_set_user_bc
659 
662  class(scalar_scheme_t), intent(inout) :: this
663  type(field_t), pointer :: nut
664  integer :: n
665  ! Factor to transform nu_t to lambda_t
666  real(kind=rp) :: lambda_factor
667 
668  lambda_factor = this%rho*this%cp/this%pr_turb
669 
670  if (this%variable_material_properties) then
671  nut => neko_field_registry%get_field(this%nut_field_name)
672  n = nut%size()
673  if (neko_bcknd_device .eq. 1) then
674  call device_cfill(this%lambda_field%x_d, this%lambda, n)
675  call device_add2s2(this%lambda_field%x_d, nut%x_d, lambda_factor, n)
676  else
677  call cfill(this%lambda_field%x, this%lambda, n)
678  call add2s2(this%lambda_field%x, nut%x, lambda_factor, n)
679  end if
680  end if
681 
683 
687  subroutine scalar_scheme_set_material_properties(this, params, user)
688  class(scalar_scheme_t), intent(inout) :: this
689  type(json_file), intent(inout) :: params
690  type(user_t), target, intent(in) :: user
691  character(len=LOG_SIZE) :: log_buf
692  ! A local pointer that is needed to make Intel happy
693  procedure(user_material_properties), pointer :: dummy_mp_ptr
694  real(kind=rp) :: dummy_mu, dummy_rho
695 
696  dummy_mp_ptr => dummy_user_material_properties
697 
698  if (.not. associated(user%material_properties, dummy_mp_ptr)) then
699 
700  write(log_buf, '(A)') "Material properties must be set in the user&
701  & file!"
702  call neko_log%message(log_buf)
703  call user%material_properties(0.0_rp, 0, dummy_rho, dummy_mu, &
704  this%cp, this%lambda, params)
705  else
706  if (params%valid_path('case.scalar.Pe') .and. &
707  (params%valid_path('case.scalar.lambda') .or. &
708  params%valid_path('case.scalar.cp'))) then
709  call neko_error("To set the material properties for the scalar,&
710  & either provide Pe OR lambda and cp in the case file.")
711  ! Non-dimensional case
712  else if (params%valid_path('case.scalar.Pe')) then
713  write(log_buf, '(A)') 'Non-dimensional scalar material properties &
714  & input.'
715  call neko_log%message(log_buf, lvl = neko_log_verbose)
716  write(log_buf, '(A)') 'Specific heat capacity will be set to 1,'
717  call neko_log%message(log_buf, lvl = neko_log_verbose)
718  write(log_buf, '(A)') 'conductivity to 1/Pe. Assumes density is 1.'
719  call neko_log%message(log_buf, lvl = neko_log_verbose)
720 
721  ! Read Pe into lambda for further manipulation.
722  call json_get(params, 'case.scalar.Pe', this%lambda)
723  write(log_buf, '(A,ES13.6)') 'Pe :', this%lambda
724  call neko_log%message(log_buf)
725 
726  ! Set cp and rho to 1 since the setup is non-dimensional.
727  this%cp = 1.0_rp
728  this%rho = 1.0_rp
729  ! Invert the Pe to get conductivity
730  this%lambda = 1.0_rp/this%lambda
731  ! Dimensional case
732  else
733  call json_get(params, 'case.scalar.lambda', this%lambda)
734  call json_get(params, 'case.scalar.cp', this%cp)
735  end if
736 
737  end if
739 
740 
741 
742 end module scalar_scheme
Abstract interface defining a dirichlet condition on a list of fields.
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 to dealocate a scalar formulation.
Abstract interface to initialize a scalar formulation.
Abstract interface to restart a scalar formulation.
Abstract interface to compute a time-step.
Abstract interface for setting material properties.
Definition: user_intf.f90:147
Abstract interface defining a user defined scalar boundary condition (pointwise) Just imitating inflo...
Definition: usr_scalar.f90:77
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_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 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 zone as a subset of facets in a mesh.
Definition: facet_zone.f90:34
Defines inflow dirichlet conditions.
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
Gather-scatter.
Implements gradient_jump_penalty_t.
Krylov preconditioner.
Definition: pc_hsmg.f90:61
Jacobi preconditioner.
Definition: pc_jacobi.f90:34
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
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 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
Defines a Neumann boundary condition.
Definition: neumann.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Krylov preconditioner.
Definition: precon.f90:34
Contains the scalar_scheme_t type.
subroutine scalar_scheme_add_bcs(this, zones, bc_labels)
Initialize boundary conditions.
subroutine scalar_scheme_free(this)
Deallocate a scalar formulation.
subroutine scalar_scheme_init(this, msh, c_Xh, gs_Xh, params, scheme, user, rho)
Initialize all related components of the current scheme.
subroutine scalar_scheme_solver_factory(ksp, n, solver, max_iter, abstol, monitor)
Initialize a linear solver.
subroutine scalar_scheme_set_user_bc(this, usr_eval)
Initialize a user defined scalar bc.
subroutine scalar_scheme_validate(this)
Validate that all fields, solvers etc necessary for performing time-stepping are defined.
subroutine scalar_scheme_update_material_properties(this)
Update the values of lambda_field if necessary.
subroutine scalar_scheme_precon_factory(pc, ksp, coef, dof, gs, bclst, pctype)
Initialize a Krylov preconditioner.
subroutine scalar_scheme_set_material_properties(this, params, user)
Set lamdba and cp.
Implements the scalar_source_term_t type.
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
Jacobi preconditioner SX-Aurora backend.
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 dirichlet conditions for scalars.
Definition: usr_scalar.f90:34
Utilities.
Definition: utils.f90:35
A list of boundary conditions.
Definition: bc.f90:104
Base type for a boundary condition.
Definition: bc.f90:51
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
User defined dirichlet condition, for which the user can work with an entire field....
field_list_t, To be able to group fields together
Definition: field_list.f90:13
Defines a jacobi preconditioner.
Definition: pc_jacobi.f90:45
Base abstract type for a canonical Krylov method, solving .
Definition: krylov.f90:66
A Neumann boundary condition for scalar fields. Sets the flux of the field to the chosen value.
Definition: neumann.f90:48
Defines a canonical Krylov preconditioner.
Definition: precon.f90:40
Base type for a scalar advection-diffusion solver.
Wrapper contaning and executing the scalar source terms.
The function space for the SEM solution fields.
Definition: space.f90:62
Defines a jacobi preconditioner for SX-Aurora.
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 scalars.
Definition: usr_scalar.f90:45