44 use krylov,
only :
ksp_t, krylov_solver_factory, krylov_solver_destroy, &
53 use precon,
only :
pc_t, precon_factory, precon_destroy
64 use json_module,
only : json_file
95 type(
gs_t),
pointer :: gs_xh
105 integer :: ksp_maxiter
107 integer :: projection_dim
109 integer :: projection_activ_step
123 integer :: n_dir_bcs = 0
125 integer :: n_neumann_bcs = 0
131 type(json_file),
pointer :: params
137 real(kind=
rp) :: lambda
141 character(len=:),
allocatable :: nut_field_name
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
165 procedure, pass(this) :: set_material_properties => &
168 procedure, pass(this) :: update_material_properties => &
183 ulag, vlag, wlag, time_scheme, rho)
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
201 real(kind=
rp),
intent(in) :: rho
212 real(kind=
rp) :: dtlag(10), tlag(10)
233 real(kind=
rp),
intent(inout) :: t
234 integer,
intent(inout) :: tstep
235 real(kind=
rp),
intent(in) :: dt
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
253 real(kind=
rp) :: dir_value, flux
256 do i = 1,
size(bc_labels)
257 bc_label = trim(bc_labels(i))
258 if (bc_label(1:2) .eq.
'd=')
then
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()
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)
291 if (bc_label(1:4) .eq.
'user')
then
292 call this%user_bc%mark_zone(zones(i))
297 do i = 1, this%n_dir_bcs
298 call bc_list_add(this%bclst_dirichlet, this%dir_bcs(i))
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))
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
328 character(len=LOG_SIZE) :: log_buf
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
341 call json_get(params,
'case.fluid.velocity_solver.type', solver_type)
342 call json_get(params,
'case.fluid.velocity_solver.preconditioner', &
344 call json_get(params,
'case.fluid.velocity_solver.absolute_tolerance', &
348 'case.fluid.velocity_solver.projection_space_size', &
349 this%projection_dim, 20)
351 'case.fluid.velocity_solver.projection_hold_steps', &
352 this%projection_activ_step, 5)
355 write(log_buf,
'(A, A)')
'Type : ', trim(scheme)
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
363 this%dm_Xh => this%u%dof
364 this%params => params
371 call this%slag%init(this%s, 2)
380 call this%set_material_properties(params, user)
382 write(log_buf,
'(A,ES13.6)')
'rho :', this%rho
384 write(log_buf,
'(A,ES13.6)')
'lambda :', this%lambda
386 write(log_buf,
'(A,ES13.6)')
'cp :', this%cp
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.
397 this%nut_field_name =
""
401 call this%lambda_field%init(this%dm_Xh,
"lambda")
404 this%lambda_field%size())
406 call cfill(this%lambda_field%x, this%lambda, this%lambda_field%size())
413 call this%user_bc%init_base(this%c_Xh)
419 this%bc_labels =
"not"
421 if (params%valid_path(
'case.scalar.boundary_types'))
then
422 call json_get(params,
'case.scalar.boundary_types', this%bc_labels, &
430 call this%f_Xh%init(this%dm_Xh, fld_name =
"scalar_rhs")
433 call this%source_term%init(this%f_Xh, this%c_Xh, user)
434 call this%source_term%add(params,
'case.scalar.source_terms')
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, &
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, &
455 if (integer_val .gt. 0)
call this%field_dir_bc%init_field(
'd_s')
457 call bc_list_add(this%bclst_dirichlet, this%field_dir_bc)
462 this%field_dir_bc%update => user%user_dirichlet_update
465 call bc_list_add(this%field_dirichlet_bcs, this%field_dir_bc)
470 'case.fluid.velocity_solver.max_iterations', &
473 'case.fluid.velocity_solver.monitor', &
474 logical_val, .false.)
476 solver_type, integer_val, solver_abstol, logical_val)
478 this%c_Xh, this%dm_Xh, this%gs_Xh, &
479 this%bclst_dirichlet, solver_precon)
483 'case.scalar.gradient_jump_penalty.enabled',&
484 this%if_gradient_jump_penalty, .false.)
486 if (this%if_gradient_jump_penalty .eqv. .true.)
then
487 if ((this%dm_Xh%xh%lx - 1) .eq. 1)
then
489 'case.scalar.gradient_jump_penalty.tau',&
490 gjp_param_a, 0.02_rp)
494 'case.scalar.gradient_jump_penalty.scaling_factor',&
497 'case.scalar.gradient_jump_penalty.scaling_exponent',&
500 call this%gradient_jump_penalty%init(params, this%dm_Xh, this%c_Xh, &
501 gjp_param_a, gjp_param_b)
519 if (
allocated(this%ksp))
then
520 call krylov_solver_destroy(this%ksp)
524 if (
allocated(this%pc))
then
525 call precon_destroy(this%pc)
529 if (
allocated(this%bc_labels))
then
530 deallocate(this%bc_labels)
533 call this%source_term%free()
538 call this%lambda_field%free()
539 call this%slag%free()
543 call this%field_dir_bc%free()
546 if (this%if_gradient_jump_penalty .eqv. .true.)
then
547 call this%gradient_jump_penalty%free()
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
564 if (.not.
allocated(this%ksp))
then
565 call neko_error(
'No Krylov solver for velocity defined')
568 if (.not.
associated(this%Xh))
then
572 if (.not.
associated(this%dm_Xh))
then
576 if (.not.
associated(this%c_Xh))
then
580 if (.not.
associated(this%f_Xh))
then
584 if (.not.
associated(this%params))
then
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
607 call krylov_solver_factory(ksp, n, solver, max_iter, &
608 abstol, monitor = monitor)
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
623 call precon_factory(pc, pctype)
625 select type (pcp => pc)
627 call pcp%init(coef, dof, gs)
629 call pcp%init(coef, dof, gs)
631 call pcp%init(coef, dof, gs)
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, &
641 call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst)
655 call this%user_bc%set_eval(usr_eval)
665 real(kind=
rp) :: lambda_factor
667 lambda_factor = this%rho*this%cp/this%pr_turb
669 if (this%variable_material_properties)
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)
676 call cfill(this%lambda_field%x, this%lambda, n)
677 call add2s2(this%lambda_field%x, nut%x, lambda_factor, n)
688 type(json_file),
intent(inout) :: params
689 type(
user_t),
target,
intent(in) :: user
690 character(len=LOG_SIZE) :: log_buf
693 real(kind=
rp) :: dummy_mu, dummy_rho
697 if (.not.
associated(user%material_properties, dummy_mp_ptr))
then
699 write(log_buf,
'(A)')
"Material properties must be set in the user&
702 call user%material_properties(0.0_rp, 0, dummy_rho, dummy_mu, &
703 this%cp, this%lambda, params)
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.")
711 else if (params%valid_path(
'case.scalar.Pe'))
then
712 write(log_buf,
'(A)')
'Non-dimensional scalar material properties &
715 write(log_buf,
'(A)')
'Specific heat capacity will be set to 1,'
717 write(log_buf,
'(A)')
'conductivity to 1/Pe. Assumes density is 1.'
721 call json_get(params,
'case.scalar.Pe', this%lambda)
722 write(log_buf,
'(A,ES13.6)')
'Pe :', this%lambda
729 this%lambda = 1.0_rp/this%lambda
732 call json_get(params,
'case.scalar.lambda', this%lambda)
733 call json_get(params,
'case.scalar.cp', this%cp)
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...
Retrieves a parameter by name or throws an error.
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.
Abstract interface defining a user defined scalar boundary condition (pointwise) Just imitating inflo...
Defines a boundary condition.
subroutine, public bc_list_add(bclst, bc)
Add a condition to a list of boundary conditions.
subroutine, public bc_list_init(bclst, size)
Constructor for a list of boundary conditions.
subroutine, public bc_list_free(bclst)
Destructor for a list of boundary conditions.
subroutine, public bc_list_apply_scalar(bclst, x, n, t, tstep)
Apply a list of boundary conditions to a scalar field.
type(mpi_comm) neko_comm
MPI communicator.
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.
Defines a mapping of the degrees of freedom.
Defines a zone as a subset of facets in a mesh.
Defines inflow dirichlet conditions.
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Implements gradient_jump_penalty_t.
Utilities for retrieving parameters from the case files.
Implements the base abstract type for Krylov solvers plus helper types.
integer, parameter, public ksp_max_iter
Maximum number of iters.
integer, parameter, public neko_log_verbose
Verbose log level.
type(log_t), public neko_log
Global log stream.
integer, parameter, public log_size
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
subroutine, public add2s2(a, b, c1, n)
Vector addition with scalar multiplication (multiplication on second argument)
integer, parameter, public neko_msh_max_zlbls
Max num. zone labels.
integer, parameter, public neko_msh_max_zlbl_len
Max length of a zone label.
integer, parameter neko_bcknd_device
Defines a Neumann boundary condition.
integer, parameter, public rp
Global precision used in computations.
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.
Defines a function space.
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.
subroutine, public dummy_user_material_properties(t, tstep, rho, mu, cp, lambda, params)
Defines dirichlet conditions for scalars.
A list of boundary conditions.
Base type for a boundary condition.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Defines a jacobi preconditioner.
Generic Dirichlet boundary condition on .
User defined dirichlet condition, for which the user can work with an entire field....
field_list_t, To be able to group fields together
Implements the gradient jump penalty.
Defines a jacobi preconditioner.
Base abstract type for a canonical Krylov method, solving .
A Neumann boundary condition for scalar fields. Sets the flux of the field to the chosen value.
Defines a canonical Krylov preconditioner.
Base type for a scalar advection-diffusion solver.
Wrapper contaning and executing the scalar source terms.
The function space for the SEM solution fields.
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.