165 subroutine scalar_pnpn_init(this, msh, coef, gs, params, user, &
166 ulag, vlag, wlag, time_scheme, rho)
167 class(scalar_pnpn_t),
target,
intent(inout) :: this
168 type(
mesh_t),
target,
intent(in) :: msh
169 type(
coef_t),
target,
intent(in) :: coef
170 type(
gs_t),
target,
intent(inout) :: gs
171 type(json_file),
target,
intent(inout) :: params
172 type(
user_t),
target,
intent(in) :: user
175 real(kind=
rp),
intent(in) :: rho
177 class(
bc_t),
pointer :: bc_i
178 character(len=15),
parameter :: scheme =
'Modular (Pn/Pn)'
184 call this%scheme_init(msh, coef, gs, params, scheme, user, rho)
187 call ax_helm_factory(this%ax, full_formulation = .false.)
190 call scalar_residual_factory(this%res)
193 call rhs_maker_ext_fctry(this%makeext)
196 call rhs_maker_bdf_fctry(this%makebdf)
199 call rhs_maker_oifs_fctry(this%makeoifs)
202 associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
203 dm_xh => this%dm_Xh, nelv => this%msh%nelv)
205 call this%s_res%init(dm_xh,
"s_res")
207 call this%abx1%init(dm_xh,
"abx1")
209 call this%abx2%init(dm_xh,
"abx2")
211 call this%advs%init(dm_xh,
"advs")
213 call this%ds%init(dm_xh,
'ds')
218 call this%setup_bcs_(user)
221 call this%bc_res%init(this%c_Xh, params)
222 do i = 1, this%bcs%size()
223 if (this%bcs%strong(i))
then
224 bc_i => this%bcs%get(i)
225 call this%bc_res%mark_facets(bc_i%marked_facet)
230 call this%bc_res%finalize()
232 call this%bclst_ds%init()
233 call this%bclst_ds%append(this%bc_res)
237 call this%proj_s%init(this%dm_Xh%size(), this%projection_dim, &
238 this%projection_activ_step)
249 call advection_factory(this%adv, params, this%c_Xh, &
250 ulag, vlag, wlag, this%chkp%dtlag, &
256 subroutine scalar_pnpn_restart(this, dtlag, tlag)
257 class(scalar_pnpn_t),
target,
intent(inout) :: this
258 real(kind=rp) :: dtlag(10), tlag(10)
262 n = this%s%dof%size()
264 call col2(this%s%x, this%c_Xh%mult, n)
265 call col2(this%slag%lf(1)%x, this%c_Xh%mult, n)
266 call col2(this%slag%lf(2)%x, this%c_Xh%mult, n)
267 if (neko_bcknd_device .eq. 1)
then
268 call device_memcpy(this%s%x, this%s%x_d, &
269 n, host_to_device, sync = .false.)
270 call device_memcpy(this%slag%lf(1)%x, this%slag%lf(1)%x_d, &
271 n, host_to_device, sync = .false.)
272 call device_memcpy(this%slag%lf(2)%x, this%slag%lf(2)%x_d, &
273 n, host_to_device, sync = .false.)
274 call device_memcpy(this%abx1%x, this%abx1%x_d, &
275 n, host_to_device, sync = .false.)
276 call device_memcpy(this%abx2%x, this%abx2%x_d, &
277 n, host_to_device, sync = .false.)
278 call device_memcpy(this%advs%x, this%advs%x_d, &
279 n, host_to_device, sync = .false.)
282 call this%gs_Xh%op(this%s, gs_op_add)
283 call this%gs_Xh%op(this%slag%lf(1), gs_op_add)
284 call this%gs_Xh%op(this%slag%lf(2), gs_op_add)
328 subroutine scalar_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
329 class(scalar_pnpn_t),
intent(inout) :: this
330 real(kind=rp),
intent(in) :: t
331 integer,
intent(in) :: tstep
332 real(kind=rp),
intent(in) :: dt
333 type(time_scheme_controller_t),
intent(in) :: ext_bdf
334 type(time_step_controller_t),
intent(in) :: dt_controller
338 type(ksp_monitor_t) :: ksp_results(1)
339 character(len=LOG_SIZE) :: log_buf
341 n = this%dm_Xh%size()
343 call profiler_start_region(
'Scalar', 2)
344 associate(u => this%u, v => this%v, w => this%w, s => this%s, &
345 cp => this%cp, rho => this%rho, &
347 s_res => this%s_res, &
348 ax => this%Ax, f_xh => this%f_Xh, xh => this%Xh, &
349 c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
350 slag => this%slag, oifs => this%oifs, &
351 lambda_field => this%lambda_field, &
352 projection_dim => this%projection_dim, &
353 msh => this%msh, res => this%res, makeoifs => this%makeoifs, &
354 makeext => this%makeext, makebdf => this%makebdf, &
355 if_variable_dt => dt_controller%if_variable_dt, &
356 dt_last_change => dt_controller%dt_last_change)
359 call print_debug(this)
361 call this%source_term%compute(t, tstep)
364 if (this%if_gradient_jump_penalty .eqv. .true.)
then
365 call this%gradient_jump_penalty%compute(u, v, w, s)
366 call this%gradient_jump_penalty%perform(f_xh)
370 call this%bcs%apply_scalar(this%f_Xh%x, dm_xh%size(), t, tstep, .false.)
374 call this%adv%compute_scalar(u, v, w, s, this%advs, &
375 xh, this%c_Xh, dm_xh%size())
377 call makeext%compute_scalar(this%abx1, this%abx2, f_xh%x, rho, &
378 ext_bdf%advection_coeffs, n)
380 call makeoifs%compute_scalar(this%advs%x, f_xh%x, rho, dt, n)
383 call this%adv%compute_scalar(u, v, w, s, f_xh, &
384 xh, this%c_Xh, dm_xh%size())
390 call makeext%compute_scalar(this%abx1, this%abx2, f_xh%x, rho, &
391 ext_bdf%advection_coeffs, n)
394 call makebdf%compute_scalar(slag, f_xh%x, s, c_xh%B, rho, dt, &
395 ext_bdf%diffusion_coeffs, ext_bdf%ndiff, n)
401 call this%bcs%apply_scalar(this%s%x, this%dm_Xh%size(), t, tstep, .true.)
404 call this%update_material_properties()
407 call profiler_start_region(
'Scalar_residual', 20)
408 call res%compute(ax, s, s_res, f_xh, c_xh, msh, xh, lambda_field, &
409 rho*cp, ext_bdf%diffusion_coeffs(1), dt, dm_xh%size())
411 call gs_xh%op(s_res, gs_op_add)
415 call this%bclst_ds%apply_scalar(s_res%x, dm_xh%size())
417 call profiler_end_region(
'Scalar_residual', 20)
419 call this%proj_s%pre_solving(s_res%x, tstep, c_xh, n, dt_controller)
421 call this%pc%update()
422 call profiler_start_region(
'Scalar_solve', 21)
423 ksp_results(1) = this%ksp%solve(ax, ds, s_res%x, n, &
424 c_xh, this%bclst_ds, gs_xh)
425 call profiler_end_region(
'Scalar_solve', 21)
427 call this%proj_s%post_solving(ds%x, ax, c_xh, this%bclst_ds, gs_xh, &
428 n, tstep, dt_controller)
431 if (neko_bcknd_device .eq. 1)
then
432 call device_add2s2(s%x_d, ds%x_d, 1.0_rp, n)
434 call add2s2(s%x, ds%x, 1.0_rp, n)
437 call scalar_step_info(tstep, t, dt, ksp_results)
440 call profiler_end_region(
'Scalar', 2)
463 subroutine scalar_pnpn_setup_bcs_(this, user)
464 class(scalar_pnpn_t),
intent(inout) :: this
465 type(user_t),
target,
intent(in) :: user
466 integer :: i, j, n_bcs, zone_size, global_zone_size, ierr
467 type(json_core) :: core
468 type(json_value),
pointer :: bc_object
469 type(json_file) :: bc_subdict
470 class(bc_t),
pointer :: bc_i
473 logical,
allocatable :: marked_zones(:)
474 integer,
allocatable :: zone_indices(:)
477 if (this%params%valid_path(
'case.scalar.boundary_conditions'))
then
478 call this%params%info(
'case.scalar.boundary_conditions', &
480 call this%params%get_core(core)
481 call this%params%get(
'case.scalar.boundary_conditions', bc_object, found)
483 call this%bcs%init(n_bcs)
485 allocate(marked_zones(
size(this%msh%labeled_zones)))
486 marked_zones = .false.
490 call json_extract_item(core, bc_object, i, bc_subdict)
495 call json_get(bc_subdict,
"zone_indices", zone_indices)
497 do j = 1,
size(zone_indices)
498 zone_size = this%msh%labeled_zones(zone_indices(j))%size
499 call mpi_allreduce(zone_size, global_zone_size, 1, &
500 mpi_integer, mpi_max, neko_comm, ierr)
502 if (global_zone_size .eq. 0)
then
503 write(error_unit,
'(A, A, I0, A, A, I0, A)')
"*** ERROR ***: ",&
504 "Zone index ", zone_indices(j), &
505 " is invalid as this zone has 0 size, meaning it ", &
506 "does not in the mesh. Check scalar boundary condition ", &
511 if (marked_zones(zone_indices(j)) .eqv. .true.)
then
512 write(error_unit,
'(A, A, I0, A, A, A, A)')
"*** ERROR ***: ", &
513 "Zone with index ", zone_indices(j), &
514 " has already been assigned a boundary condition. ", &
515 "Please check your boundary_conditions entry for the ", &
516 "scalar and make sure that each zone index appears only ",&
517 "in a single boundary condition."
520 marked_zones(zone_indices(j)) = .true.
526 call bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
527 call this%bcs%append(bc_i)
531 do i = 1,
size(this%msh%labeled_zones)
532 if ((this%msh%labeled_zones(i)%size .gt. 0) .and. &
533 (marked_zones(i) .eqv. .false.))
then
534 write(error_unit,
'(A, A, I0)')
"*** ERROR ***: ", &
535 "No scalar boundary condition assigned to zone ", i