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)
234 real(kind=
rp),
intent(inout) :: t
235 integer,
intent(inout) :: tstep
236 real(kind=
rp),
intent(in) :: dt
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
254 real(kind=
rp) :: dir_value, flux
257 do i = 1,
size(bc_labels)
258 bc_label = trim(bc_labels(i))
259 if (bc_label(1:2) .eq.
'd=')
then
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()
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)
292 if (bc_label(1:4) .eq.
'user')
then
293 call this%user_bc%mark_zone(zones(i))
298 do i = 1, this%n_dir_bcs
299 call bc_list_add(this%bclst_dirichlet, this%dir_bcs(i))
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))
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
329 character(len=LOG_SIZE) :: log_buf
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
342 call json_get(params,
'case.fluid.velocity_solver.type', solver_type)
343 call json_get(params,
'case.fluid.velocity_solver.preconditioner', &
345 call json_get(params,
'case.fluid.velocity_solver.absolute_tolerance', &
349 'case.fluid.velocity_solver.projection_space_size', &
350 this%projection_dim, 20)
352 'case.fluid.velocity_solver.projection_hold_steps', &
353 this%projection_activ_step, 5)
356 write(log_buf,
'(A, A)')
'Type : ', trim(scheme)
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
364 this%dm_Xh => this%u%dof
365 this%params => params
372 call this%slag%init(this%s, 2)
381 call this%set_material_properties(params, user)
383 write(log_buf,
'(A,ES13.6)')
'rho :', this%rho
385 write(log_buf,
'(A,ES13.6)')
'lambda :', this%lambda
387 write(log_buf,
'(A,ES13.6)')
'cp :', this%cp
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.
398 this%nut_field_name =
""
402 call this%lambda_field%init(this%dm_Xh,
"lambda")
405 this%lambda_field%size())
407 call cfill(this%lambda_field%x, this%lambda, this%lambda_field%size())
414 call this%user_bc%init_base(this%c_Xh)
420 this%bc_labels =
"not"
422 if (params%valid_path(
'case.scalar.boundary_types'))
then
423 call json_get(params,
'case.scalar.boundary_types', this%bc_labels, &
431 call this%f_Xh%init(this%dm_Xh, fld_name =
"scalar_rhs")
434 call this%source_term%init(this%f_Xh, this%c_Xh, user)
435 call this%source_term%add(params,
'case.scalar.source_terms')
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, &
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, &
456 if (integer_val .gt. 0)
call this%field_dir_bc%init_field(
'd_s')
458 call bc_list_add(this%bclst_dirichlet, this%field_dir_bc)
463 this%field_dir_bc%update => user%user_dirichlet_update
466 call bc_list_add(this%field_dirichlet_bcs, this%field_dir_bc)
471 'case.fluid.velocity_solver.max_iterations', &
474 'case.fluid.velocity_solver.monitor', &
475 logical_val, .false.)
477 solver_type, integer_val, solver_abstol, logical_val)
479 this%c_Xh, this%dm_Xh, this%gs_Xh, &
480 this%bclst_dirichlet, solver_precon)
484 'case.scalar.gradient_jump_penalty.enabled',&
485 this%if_gradient_jump_penalty, .false.)
487 if (this%if_gradient_jump_penalty .eqv. .true.)
then
488 if ((this%dm_Xh%xh%lx - 1) .eq. 1)
then
490 'case.scalar.gradient_jump_penalty.tau',&
491 gjp_param_a, 0.02_rp)
495 'case.scalar.gradient_jump_penalty.scaling_factor',&
498 'case.scalar.gradient_jump_penalty.scaling_exponent',&
501 call this%gradient_jump_penalty%init(params, this%dm_Xh, this%c_Xh, &
502 gjp_param_a, gjp_param_b)
520 if (
allocated(this%ksp))
then
521 call krylov_solver_destroy(this%ksp)
525 if (
allocated(this%pc))
then
526 call precon_destroy(this%pc)
530 if (
allocated(this%bc_labels))
then
531 deallocate(this%bc_labels)
534 call this%source_term%free()
539 call this%lambda_field%free()
540 call this%slag%free()
544 call this%field_dir_bc%free()
547 if (this%if_gradient_jump_penalty .eqv. .true.)
then
548 call this%gradient_jump_penalty%free()
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
565 if (.not.
allocated(this%ksp))
then
566 call neko_error(
'No Krylov solver for velocity defined')
569 if (.not.
associated(this%Xh))
then
573 if (.not.
associated(this%dm_Xh))
then
577 if (.not.
associated(this%c_Xh))
then
581 if (.not.
associated(this%f_Xh))
then
585 if (.not.
associated(this%params))
then
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
608 call krylov_solver_factory(ksp, n, solver, max_iter, &
609 abstol, monitor = monitor)
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
624 call precon_factory(pc, pctype)
626 select type (pcp => pc)
628 call pcp%init(coef, dof, gs)
630 call pcp%init(coef, dof, gs)
632 call pcp%init(coef, dof, gs)
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, &
642 call pcp%init(dof%msh, dof%Xh, coef, dof, gs, bclst)
656 call this%user_bc%set_eval(usr_eval)
666 real(kind=
rp) :: lambda_factor
668 lambda_factor = this%rho*this%cp/this%pr_turb
670 if (this%variable_material_properties)
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)
677 call cfill(this%lambda_field%x, this%lambda, n)
678 call add2s2(this%lambda_field%x, nut%x, lambda_factor, n)
689 type(json_file),
intent(inout) :: params
690 type(
user_t),
target,
intent(in) :: user
691 character(len=LOG_SIZE) :: log_buf
694 real(kind=
rp) :: dummy_mu, dummy_rho
698 if (.not.
associated(user%material_properties, dummy_mp_ptr))
then
700 write(log_buf,
'(A)')
"Material properties must be set in the user&
703 call user%material_properties(0.0_rp, 0, dummy_rho, dummy_mu, &
704 this%cp, this%lambda, params)
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.")
712 else if (params%valid_path(
'case.scalar.Pe'))
then
713 write(log_buf,
'(A)')
'Non-dimensional scalar material properties &
716 write(log_buf,
'(A)')
'Specific heat capacity will be set to 1,'
718 write(log_buf,
'(A)')
'conductivity to 1/Pe. Assumes density is 1.'
722 call json_get(params,
'case.scalar.Pe', this%lambda)
723 write(log_buf,
'(A,ES13.6)')
'Pe :', this%lambda
730 this%lambda = 1.0_rp/this%lambda
733 call json_get(params,
'case.scalar.lambda', this%lambda)
734 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.