Neko  0.8.99
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 
217  abstract interface
218  subroutine scalar_scheme_free_intrf(this)
219  import scalar_scheme_t
220  class(scalar_scheme_t), intent(inout) :: this
221  end subroutine scalar_scheme_free_intrf
222  end interface
223 
225  abstract interface
226  subroutine scalar_scheme_step_intrf(this, t, tstep, dt, ext_bdf, &
227  dt_controller)
228  import scalar_scheme_t
231  import rp
232  class(scalar_scheme_t), intent(inout) :: this
233  real(kind=rp), intent(inout) :: t
234  integer, intent(inout) :: tstep
235  real(kind=rp), intent(in) :: dt
236  type(time_scheme_controller_t), intent(inout) :: ext_bdf
237  type(time_step_controller_t), intent(in) :: dt_controller
238  end subroutine scalar_scheme_step_intrf
239  end interface
240 
241 contains
242 
247  subroutine scalar_scheme_add_bcs(this, zones, bc_labels)
248  class(scalar_scheme_t), intent(inout) :: this
249  type(facet_zone_t), intent(inout) :: zones(NEKO_MSH_MAX_ZLBLS)
250  character(len=NEKO_MSH_MAX_ZLBL_LEN), intent(in) :: bc_labels(:)
251  character(len=NEKO_MSH_MAX_ZLBL_LEN) :: bc_label
252  integer :: i
253  real(kind=rp) :: dir_value, flux
254  logical :: bc_exists
255 
256  do i = 1, size(bc_labels)
257  bc_label = trim(bc_labels(i))
258  if (bc_label(1:2) .eq. 'd=') then
259 ! The idea of this commented piece of code is to merge bcs with the same
260 ! Dirichlet value into 1 so that one has less kernel launches. Currently
261 ! segfaults, needs investigation.
262 ! bc_exists = .false.
263 ! bc_idx = 0
264 ! do j = 1, i-1
265 ! if (bc_label .eq. bc_labels(j)) then
266 ! bc_exists = .true.
267 ! bc_idx = j
268 ! end if
269 ! end do
270 
271 ! if (bc_exists) then
272 ! call this%dir_bcs(j)%mark_zone(zones(i))
273 ! else
274  this%n_dir_bcs = this%n_dir_bcs + 1
275  call this%dir_bcs(this%n_dir_bcs)%init_base(this%c_Xh)
276  call this%dir_bcs(this%n_dir_bcs)%mark_zone(zones(i))
277  read(bc_label(3:), *) dir_value
278  call this%dir_bcs(this%n_dir_bcs)%set_g(dir_value)
279  call this%dir_bcs(this%n_dir_bcs)%finalize()
280  end if
281 
282  if (bc_label(1:2) .eq. 'n=') then
283  this%n_neumann_bcs = this%n_neumann_bcs + 1
284  call this%neumann_bcs(this%n_neumann_bcs)%init_base(this%c_Xh)
285  call this%neumann_bcs(this%n_neumann_bcs)%mark_zone(zones(i))
286  read(bc_label(3:), *) flux
287  call this%neumann_bcs(this%n_neumann_bcs)%finalize_neumann(flux)
288  end if
289 
291  if (bc_label(1:4) .eq. 'user') then
292  call this%user_bc%mark_zone(zones(i))
293  end if
294 
295  end do
296 
297  do i = 1, this%n_dir_bcs
298  call bc_list_add(this%bclst_dirichlet, this%dir_bcs(i))
299  end do
300 
301  ! Create list with just Neumann bcs
302  call bc_list_init(this%bclst_neumann, this%n_neumann_bcs)
303  do i = 1, this%n_neumann_bcs
304  call bc_list_add(this%bclst_neumann, this%neumann_bcs(i))
305  end do
306 
307  end subroutine scalar_scheme_add_bcs
308 
317  subroutine scalar_scheme_init(this, msh, c_Xh, gs_Xh, params, scheme, user, &
318  rho)
319  class(scalar_scheme_t), target, intent(inout) :: this
320  type(mesh_t), target, intent(inout) :: msh
321  type(coef_t), target, intent(inout) :: c_Xh
322  type(gs_t), target, intent(inout) :: gs_Xh
323  type(json_file), target, intent(inout) :: params
324  character(len=*), intent(in) :: scheme
325  type(user_t), target, intent(in) :: user
326  real(kind=rp), intent(in) :: rho
327  ! IO buffer for log output
328  character(len=LOG_SIZE) :: log_buf
329  ! Variables for retrieving json parameters
330  logical :: logical_val
331  real(kind=rp) :: real_val, solver_abstol
332  integer :: integer_val, ierr
333  character(len=:), allocatable :: solver_type, solver_precon
334  real(kind=rp) :: gjp_param_a, gjp_param_b
335 
336  this%u => neko_field_registry%get_field('u')
337  this%v => neko_field_registry%get_field('v')
338  this%w => neko_field_registry%get_field('w')
339 
340  call neko_log%section('Scalar')
341  call json_get(params, 'case.fluid.velocity_solver.type', solver_type)
342  call json_get(params, 'case.fluid.velocity_solver.preconditioner', &
343  solver_precon)
344  call json_get(params, 'case.fluid.velocity_solver.absolute_tolerance', &
345  solver_abstol)
346 
347  call json_get_or_default(params, &
348  'case.fluid.velocity_solver.projection_space_size', &
349  this%projection_dim, 20)
350  call json_get_or_default(params, &
351  'case.fluid.velocity_solver.projection_hold_steps', &
352  this%projection_activ_step, 5)
353 
354 
355  write(log_buf, '(A, A)') 'Type : ', trim(scheme)
356  call neko_log%message(log_buf)
357  call neko_log%message('Ksp scalar : ('// trim(solver_type) // &
358  ', ' // trim(solver_precon) // ')')
359  write(log_buf, '(A,ES13.6)') ' `-abs tol :', solver_abstol
360  call neko_log%message(log_buf)
361 
362  this%Xh => this%u%Xh
363  this%dm_Xh => this%u%dof
364  this%params => params
365  this%msh => msh
366  if (.not. neko_field_registry%field_exists('s')) then
367  call neko_field_registry%add_field(this%dm_Xh, 's')
368  end if
369  this%s => neko_field_registry%get_field('s')
370 
371  call this%slag%init(this%s, 2)
372 
373  this%gs_Xh => gs_xh
374  this%c_Xh => c_xh
375 
376  !
377  ! Material properties
378  !
379  this%rho = rho
380  call this%set_material_properties(params, user)
381 
382  write(log_buf, '(A,ES13.6)') 'rho :', this%rho
383  call neko_log%message(log_buf)
384  write(log_buf, '(A,ES13.6)') 'lambda :', this%lambda
385  call neko_log%message(log_buf)
386  write(log_buf, '(A,ES13.6)') 'cp :', this%cp
387  call neko_log%message(log_buf)
388 
389  !
390  ! Turbulence modelling and variable material properties
391  !
392  if (params%valid_path('case.scalar.nut_field')) then
393  call json_get(params, 'case.scalar.Pr_t', this%pr_turb)
394  call json_get(params, 'case.scalar.nut_field', this%nut_field_name)
395  this%variable_material_properties = .true.
396  else
397  this%nut_field_name = ""
398  end if
399 
400  ! Fill lambda field with the physical value
401  call this%lambda_field%init(this%dm_Xh, "lambda")
402  if (neko_bcknd_device .eq. 1) then
403  call device_cfill(this%lambda_field%x_d, this%lambda, &
404  this%lambda_field%size())
405  else
406  call cfill(this%lambda_field%x, this%lambda, this%lambda_field%size())
407  end if
408 
409  !
410  ! Setup scalar boundary conditions
411  !
412  call bc_list_init(this%bclst_dirichlet)
413  call this%user_bc%init_base(this%c_Xh)
414 
415  ! Read boundary types from the case file
416  allocate(this%bc_labels(neko_msh_max_zlbls))
417 
418  ! A filler value
419  this%bc_labels = "not"
420 
421  if (params%valid_path('case.scalar.boundary_types')) then
422  call json_get(params, 'case.scalar.boundary_types', this%bc_labels, &
423  'not')
424  end if
425 
426  !
427  ! Setup right-hand side field.
428  !
429  allocate(this%f_Xh)
430  call this%f_Xh%init(this%dm_Xh, fld_name = "scalar_rhs")
431 
432  ! Initialize the source term
433  call this%source_term%init(this%f_Xh, this%c_Xh, user)
434  call this%source_term%add(params, 'case.scalar.source_terms')
435 
436  call scalar_scheme_add_bcs(this, msh%labeled_zones, this%bc_labels)
437 
438  ! Mark BC zones
439  call this%user_bc%mark_zone(msh%wall)
440  call this%user_bc%mark_zone(msh%inlet)
441  call this%user_bc%mark_zone(msh%outlet)
442  call this%user_bc%mark_zone(msh%outlet_normal)
443  call this%user_bc%mark_zone(msh%sympln)
444  call this%user_bc%finalize()
445  if (this%user_bc%msk(0) .gt. 0) call bc_list_add(this%bclst_dirichlet, &
446  this%user_bc)
447 
448  ! Add field dirichlet BCs
449  call this%field_dir_bc%init_base(this%c_Xh)
450  call this%field_dir_bc%mark_zones_from_list(msh%labeled_zones, &
451  'd_s', this%bc_labels)
452  call this%field_dir_bc%finalize()
453  call mpi_allreduce(this%field_dir_bc%msk(0), integer_val, 1, &
454  mpi_integer, mpi_sum, neko_comm, ierr)
455  if (integer_val .gt. 0) call this%field_dir_bc%init_field('d_s')
456 
457  call bc_list_add(this%bclst_dirichlet, this%field_dir_bc)
458 
459  !
460  ! Associate our field dirichlet update to the user one.
461  !
462  this%field_dir_bc%update => user%user_dirichlet_update
463 
464  call bc_list_init(this%field_dirichlet_bcs, size = 1)
465  call bc_list_add(this%field_dirichlet_bcs, this%field_dir_bc)
466 
467 
468  ! todo parameter file ksp tol should be added
469  call json_get_or_default(params, &
470  'case.fluid.velocity_solver.max_iterations', &
471  integer_val, ksp_max_iter)
472  call json_get_or_default(params, &
473  'case.fluid.velocity_solver.monitor', &
474  logical_val, .false.)
475  call scalar_scheme_solver_factory(this%ksp, this%dm_Xh%size(), &
476  solver_type, integer_val, solver_abstol, logical_val)
477  call scalar_scheme_precon_factory(this%pc, this%ksp, &
478  this%c_Xh, this%dm_Xh, this%gs_Xh, &
479  this%bclst_dirichlet, solver_precon)
480 
481  ! Initiate gradient jump penalty
482  call json_get_or_default(params, &
483  'case.scalar.gradient_jump_penalty.enabled',&
484  this%if_gradient_jump_penalty, .false.)
485 
486  if (this%if_gradient_jump_penalty .eqv. .true.) then
487  if ((this%dm_Xh%xh%lx - 1) .eq. 1) then
488  call json_get_or_default(params, &
489  'case.scalar.gradient_jump_penalty.tau',&
490  gjp_param_a, 0.02_rp)
491  gjp_param_b = 0.0_rp
492  else
493  call json_get_or_default(params, &
494  'case.scalar.gradient_jump_penalty.scaling_factor',&
495  gjp_param_a, 0.8_rp)
496  call json_get_or_default(params, &
497  'case.scalar.gradient_jump_penalty.scaling_exponent',&
498  gjp_param_b, 4.0_rp)
499  end if
500  call this%gradient_jump_penalty%init(params, this%dm_Xh, this%c_Xh, &
501  gjp_param_a, gjp_param_b)
502  end if
503 
504  call neko_log%end_section()
505 
506  end subroutine scalar_scheme_init
507 
508 
510  subroutine scalar_scheme_free(this)
511  class(scalar_scheme_t), intent(inout) :: this
512 
513  nullify(this%Xh)
514  nullify(this%dm_Xh)
515  nullify(this%gs_Xh)
516  nullify(this%c_Xh)
517  nullify(this%params)
518 
519  if (allocated(this%ksp)) then
520  call krylov_solver_destroy(this%ksp)
521  deallocate(this%ksp)
522  end if
523 
524  if (allocated(this%pc)) then
525  call precon_destroy(this%pc)
526  deallocate(this%pc)
527  end if
528 
529  if (allocated(this%bc_labels)) then
530  deallocate(this%bc_labels)
531  end if
532 
533  call this%source_term%free()
534 
535  call bc_list_free(this%bclst_dirichlet)
536  call bc_list_free(this%bclst_neumann)
537 
538  call this%lambda_field%free()
539  call this%slag%free()
540 
541  ! Free everything related to field dirichlet BCs
542  call bc_list_free(this%field_dirichlet_bcs)
543  call this%field_dir_bc%free()
544 
545  ! Free gradient jump penalty
546  if (this%if_gradient_jump_penalty .eqv. .true.) then
547  call this%gradient_jump_penalty%free()
548  end if
549 
550  end subroutine scalar_scheme_free
551 
554  subroutine scalar_scheme_validate(this)
555  class(scalar_scheme_t), target, intent(inout) :: this
556 
557  if ( (.not. allocated(this%u%x)) .or. &
558  (.not. allocated(this%v%x)) .or. &
559  (.not. allocated(this%w%x)) .or. &
560  (.not. allocated(this%s%x))) then
561  call neko_error('Fields are not allocated')
562  end if
563 
564  if (.not. allocated(this%ksp)) then
565  call neko_error('No Krylov solver for velocity defined')
566  end if
567 
568  if (.not. associated(this%Xh)) then
569  call neko_error('No function space defined')
570  end if
571 
572  if (.not. associated(this%dm_Xh)) then
573  call neko_error('No dofmap defined')
574  end if
575 
576  if (.not. associated(this%c_Xh)) then
577  call neko_error('No coefficients defined')
578  end if
579 
580  if (.not. associated(this%f_Xh)) then
581  call neko_error('No rhs allocated')
582  end if
583 
584  if (.not. associated(this%params)) then
585  call neko_error('No parameters defined')
586  end if
587 
588  !
589  ! Setup checkpoint structure (if everything is fine)
590  !
591 ! @todo no io for now
592 ! call this%chkp%init(this%u, this%v, this%w, this%p)
593 
594  end subroutine scalar_scheme_validate
595 
598  subroutine scalar_scheme_solver_factory(ksp, n, solver, max_iter, &
599  abstol, monitor)
600  class(ksp_t), allocatable, target, intent(inout) :: ksp
601  integer, intent(in), value :: n
602  integer, intent(in) :: max_iter
603  character(len=*), intent(in) :: solver
604  real(kind=rp) :: abstol
605  logical, intent(in) :: monitor
606 
607  call krylov_solver_factory(ksp, n, solver, max_iter, &
608  abstol, monitor = monitor)
609 
610  end subroutine scalar_scheme_solver_factory
611 
613  subroutine scalar_scheme_precon_factory(pc, ksp, coef, dof, gs, bclst, &
614  pctype)
615  class(pc_t), allocatable, target, intent(inout) :: pc
616  class(ksp_t), target, intent(inout) :: ksp
617  type(coef_t), target, intent(inout) :: coef
618  type(dofmap_t), target, intent(inout) :: dof
619  type(gs_t), target, intent(inout) :: gs
620  type(bc_list_t), target, intent(inout) :: bclst
621  character(len=*) :: pctype
622 
623  call precon_factory(pc, pctype)
624 
625  select type (pcp => pc)
626  type is (jacobi_t)
627  call pcp%init(coef, dof, gs)
628  type is (sx_jacobi_t)
629  call pcp%init(coef, dof, gs)
630  type is (device_jacobi_t)
631  call pcp%init(coef, dof, gs)
632  type is (hsmg_t)
633  if (len_trim(pctype) .gt. 4) then
634  if (index(pctype, '+') .eq. 5) then
635  call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst, &
636  trim(pctype(6:)))
637  else
638  call neko_error('Unknown coarse grid solver')
639  end if
640  else
641  call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst)
642  end if
643  end select
644 
645  call ksp%set_pc(pc)
646 
647  end subroutine scalar_scheme_precon_factory
648 
651  subroutine scalar_scheme_set_user_bc(this, usr_eval)
652  class(scalar_scheme_t), intent(inout) :: this
653  procedure(usr_scalar_bc_eval) :: usr_eval
654 
655  call this%user_bc%set_eval(usr_eval)
656 
657  end subroutine scalar_scheme_set_user_bc
658 
661  class(scalar_scheme_t), intent(inout) :: this
662  type(field_t), pointer :: nut
663  integer :: n
664  ! Factor to transform nu_t to lambda_t
665  real(kind=rp) :: lambda_factor
666 
667  lambda_factor = this%rho*this%cp/this%pr_turb
668 
669  if (this%variable_material_properties) then
670  nut => neko_field_registry%get_field(this%nut_field_name)
671  n = nut%size()
672  if (neko_bcknd_device .eq. 1) then
673  call device_cfill(this%lambda_field%x_d, this%lambda, n)
674  call device_add2s2(this%lambda_field%x_d, nut%x_d, lambda_factor, n)
675  else
676  call cfill(this%lambda_field%x, this%lambda, n)
677  call add2s2(this%lambda_field%x, nut%x, lambda_factor, n)
678  end if
679  end if
680 
682 
686  subroutine scalar_scheme_set_material_properties(this, params, user)
687  class(scalar_scheme_t), intent(inout) :: this
688  type(json_file), intent(inout) :: params
689  type(user_t), target, intent(in) :: user
690  character(len=LOG_SIZE) :: log_buf
691  ! A local pointer that is needed to make Intel happy
692  procedure(user_material_properties), pointer :: dummy_mp_ptr
693  real(kind=rp) :: dummy_mu, dummy_rho
694 
695  dummy_mp_ptr => dummy_user_material_properties
696 
697  if (.not. associated(user%material_properties, dummy_mp_ptr)) then
698 
699  write(log_buf, '(A)') "Material properties must be set in the user&
700  & file!"
701  call neko_log%message(log_buf)
702  call user%material_properties(0.0_rp, 0, dummy_rho, dummy_mu, &
703  this%cp, this%lambda, params)
704  else
705  if (params%valid_path('case.scalar.Pe') .and. &
706  (params%valid_path('case.scalar.lambda') .or. &
707  params%valid_path('case.scalar.cp'))) then
708  call neko_error("To set the material properties for the scalar,&
709  & either provide Pe OR lambda and cp in the case file.")
710  ! Non-dimensional case
711  else if (params%valid_path('case.scalar.Pe')) then
712  write(log_buf, '(A)') 'Non-dimensional scalar material properties &
713  & input.'
714  call neko_log%message(log_buf, lvl = neko_log_verbose)
715  write(log_buf, '(A)') 'Specific heat capacity will be set to 1,'
716  call neko_log%message(log_buf, lvl = neko_log_verbose)
717  write(log_buf, '(A)') 'conductivity to 1/Pe. Assumes density is 1.'
718  call neko_log%message(log_buf, lvl = neko_log_verbose)
719 
720  ! Read Pe into lambda for further manipulation.
721  call json_get(params, 'case.scalar.Pe', this%lambda)
722  write(log_buf, '(A,ES13.6)') 'Pe :', this%lambda
723  call neko_log%message(log_buf)
724 
725  ! Set cp and rho to 1 since the setup is non-dimensional.
726  this%cp = 1.0_rp
727  this%rho = 1.0_rp
728  ! Invert the Pe to get conductivity
729  this%lambda = 1.0_rp/this%lambda
730  ! Dimensional case
731  else
732  call json_get(params, 'case.scalar.lambda', this%lambda)
733  call json_get(params, 'case.scalar.cp', this%cp)
734  end if
735 
736  end if
738 
739 
740 
741 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