36 use comm,
only: neko_comm
37 use utils,
only: neko_error
38 use num_types,
only: rp
39 use,
intrinsic :: iso_fortran_env, only: error_unit
40 use rhs_maker,
only : rhs_maker_bdf_t, rhs_maker_ext_t, rhs_maker_oifs_t, &
41 rhs_maker_ext_fctry, rhs_maker_bdf_fctry, rhs_maker_oifs_fctry
43 use checkpoint,
only : chkp_t
44 use field,
only : field_t
45 use bc_list,
only : bc_list_t
46 use mesh,
only : mesh_t
47 use coefs,
only : coef_t
48 use device,
only : host_to_device, device_memcpy
49 use gather_scatter,
only : gs_t, gs_op_add
50 use scalar_residual,
only : scalar_residual_t, scalar_residual_factory
51 use ax_product,
only : ax_t, ax_helm_factory
52 use field_series,
only: field_series_t
53 use facet_normal,
only : facet_normal_t
54 use krylov,
only : ksp_monitor_t
55 use device_math,
only : device_add2s2, device_col2
56 use scalar_aux,
only : scalar_step_info
57 use time_scheme_controller,
only : time_scheme_controller_t
58 use projection,
only : projection_t
59 use math,
only : glsc2, col2, add2s2
60 use logger,
only : neko_log, log_size, neko_log_debug
62 use profiler,
only : profiler_start_region, profiler_end_region
63 use json_utils,
only : json_get, json_get_or_default, json_extract_item
64 use json_module,
only : json_file, json_core, json_value
65 use user_intf,
only : user_t
66 use neko_config,
only : neko_bcknd_device
67 use zero_dirichlet,
only : zero_dirichlet_t
68 use time_step_controller,
only : time_step_controller_t
69 use scratch_registry,
only : neko_scratch_registry
70 use time_state,
only : time_state_t
72 use mpi_f08,
only: mpi_integer, mpi_sum, mpi_max
80 type(field_t) :: s_adj_res
83 type(field_t) :: ds_adj
86 class(ax_t),
allocatable :: ax
89 type(projection_t) :: proj_s
94 type(zero_dirichlet_t) :: bc_res
100 type(bc_list_t) :: bclst_ds
109 type(field_t) :: advs
112 class(scalar_residual_t),
allocatable :: res
115 class(rhs_maker_ext_t),
allocatable :: makeext
118 class(rhs_maker_bdf_t),
allocatable :: makebdf
121 class(rhs_maker_oifs_t),
allocatable :: makeoifs
125 procedure, pass(this) :: init => adjoint_scalar_pnpn_init
127 procedure, pass(this) :: restart => adjoint_scalar_pnpn_restart
129 procedure, pass(this) :: free => adjoint_scalar_pnpn_free
131 procedure, pass(this) :: step => adjoint_scalar_pnpn_step
133 procedure, pass(this) :: setup_bcs_ => adjoint_scalar_pnpn_setup_bcs_
143 module subroutine adjoint_bc_factory(object, scheme, json, coef, user)
144 class(bc_t),
pointer,
intent(inout) :: object
145 type(adjoint_scalar_pnpn_t),
intent(in) :: scheme
146 type(json_file),
intent(inout) :: json
147 type(coef_t),
intent(in) :: coef
148 type(user_t),
intent(in) :: user
149 end subroutine adjoint_bc_factory
150 end interface adjoint_bc_factory
170 subroutine adjoint_scalar_pnpn_init(this, msh, coef, gs, params_adjoint, &
171 params_primal, numerics_params, user, chkp, ulag, vlag, wlag, &
173 class(adjoint_scalar_pnpn_t),
target,
intent(inout) :: this
174 type(mesh_t),
target,
intent(in) :: msh
175 type(coef_t),
target,
intent(in) :: coef
176 type(gs_t),
target,
intent(inout) :: gs
177 type(json_file),
target,
intent(inout) :: params_adjoint
178 type(json_file),
target,
intent(inout) :: params_primal
179 type(json_file),
target,
intent(inout) :: numerics_params
180 type(user_t),
target,
intent(in) :: user
181 type(chkp_t),
target,
intent(inout) :: chkp
182 type(field_series_t),
target,
intent(in) :: ulag, vlag, wlag
183 type(time_scheme_controller_t),
target,
intent(in) :: time_scheme
184 type(field_t),
target,
intent(in) :: rho
186 class(bc_t),
pointer :: bc_i
187 character(len=15),
parameter :: scheme =
'Modular (Pn/Pn)'
193 call this%scheme_init(msh, coef, gs, params_adjoint, params_primal, &
197 call ax_helm_factory(this%ax, full_formulation = .false.)
200 call scalar_residual_factory(this%res)
203 call rhs_maker_ext_fctry(this%makeext)
206 call rhs_maker_bdf_fctry(this%makebdf)
209 call rhs_maker_oifs_fctry(this%makeoifs)
212 associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
213 dm_xh => this%dm_Xh, nelv => this%msh%nelv)
215 call this%s_adj_res%init(dm_xh,
"s_adj_res")
217 call this%abx1%init(dm_xh,
"abx1")
219 call this%abx2%init(dm_xh,
"abx2")
221 call this%advs%init(dm_xh,
"advs")
223 call this%ds_adj%init(dm_xh,
'ds_adj')
228 call this%setup_bcs_(user)
231 call this%bc_res%init(this%c_Xh, params_adjoint)
232 do i = 1, this%bcs%size()
233 if (this%bcs%strong(i))
then
234 bc_i => this%bcs%get(i)
235 call this%bc_res%mark_facets(bc_i%marked_facet)
240 call this%bc_res%finalize()
242 call this%bclst_ds%init()
243 call this%bclst_ds%append(this%bc_res)
247 call this%proj_s%init(this%dm_Xh%size(), this%projection_dim, &
248 this%projection_activ_step)
251 call json_get_or_default(numerics_params,
'oifs', this%oifs, .false.)
257 call json_get_or_default(params_adjoint,
'advection', advection, .true.)
259 call advection_adjoint_factory(this%adv, numerics_params, this%c_Xh, &
260 ulag, vlag, wlag, this%chkp%dtlag, &
261 this%chkp%tlag, time_scheme, .not. advection, &
272 end subroutine adjoint_scalar_pnpn_init
275 subroutine adjoint_scalar_pnpn_restart(this, chkp)
276 class(adjoint_scalar_pnpn_t),
target,
intent(inout) :: this
277 type(chkp_t),
intent(inout) :: chkp
278 real(kind=rp) :: dtlag(10), tlag(10)
283 n = this%s_adj%dof%size()
285 call col2(this%s_adj%x, this%c_Xh%mult, n)
286 call col2(this%s_adj_lag%lf(1)%x, this%c_Xh%mult, n)
287 call col2(this%s_adj_lag%lf(2)%x, this%c_Xh%mult, n)
288 if (neko_bcknd_device .eq. 1)
then
289 call device_memcpy(this%s_adj%x, this%s_adj%x_d, &
290 n, host_to_device, sync = .false.)
291 call device_memcpy(this%s_adj_lag%lf(1)%x, this%s_adj_lag%lf(1)%x_d, &
292 n, host_to_device, sync = .false.)
293 call device_memcpy(this%s_adj_lag%lf(2)%x, this%s_adj_lag%lf(2)%x_d, &
294 n, host_to_device, sync = .false.)
295 call device_memcpy(this%abx1%x, this%abx1%x_d, &
296 n, host_to_device, sync = .false.)
297 call device_memcpy(this%abx2%x, this%abx2%x_d, &
298 n, host_to_device, sync = .false.)
299 call device_memcpy(this%advs%x, this%advs%x_d, &
300 n, host_to_device, sync = .false.)
303 call this%gs_Xh%op(this%s_adj, gs_op_add)
304 call this%gs_Xh%op(this%s_adj_lag%lf(1), gs_op_add)
305 call this%gs_Xh%op(this%s_adj_lag%lf(2), gs_op_add)
307 end subroutine adjoint_scalar_pnpn_restart
309 subroutine adjoint_scalar_pnpn_free(this)
310 class(adjoint_scalar_pnpn_t),
intent(inout) :: this
313 call this%scheme_free()
315 call this%bclst_ds%free()
316 call this%proj_s%free()
318 call this%s_adj_res%free()
320 call this%ds_adj%free()
322 call this%abx1%free()
323 call this%abx2%free()
325 call this%advs%free()
327 if (
allocated(this%Ax))
then
331 if (
allocated(this%res))
then
335 if (
allocated(this%makeext))
then
336 deallocate(this%makeext)
339 if (
allocated(this%makebdf))
then
340 deallocate(this%makebdf)
343 if (
allocated(this%makeoifs))
then
344 deallocate(this%makeoifs)
347 end subroutine adjoint_scalar_pnpn_free
349 subroutine adjoint_scalar_pnpn_step(this, time, ext_bdf, dt_controller)
350 class(adjoint_scalar_pnpn_t),
intent(inout) :: this
351 type(time_state_t),
intent(in) :: time
352 type(time_scheme_controller_t),
intent(in) :: ext_bdf
353 type(time_step_controller_t),
intent(in) :: dt_controller
357 type(ksp_monitor_t) :: ksp_results(1)
359 n = this%dm_Xh%size()
361 call profiler_start_region(
'Scalar', 2)
362 associate(u => this%u, v => this%v, w => this%w, s_adj => this%s_adj, &
363 cp => this%cp, rho => this%rho, lambda => this%lambda, &
364 ds_adj => this%ds_adj, &
365 s_adj_res => this%s_adj_res, &
366 ax => this%Ax, f_xh => this%f_Xh, xh => this%Xh, &
367 c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
368 s_adj_lag => this%s_adj_lag, oifs => this%oifs, &
369 projection_dim => this%projection_dim, &
370 msh => this%msh, res => this%res, makeoifs => this%makeoifs, &
371 makeext => this%makeext, makebdf => this%makebdf, &
372 t => time%t, tstep => time%tstep, dt => time%dt)
375 call print_debug(this)
377 call this%source_term%compute(time)
380 call this%bcs%apply_scalar(this%f_Xh%x, dm_xh%size(), time, .false.)
394 call this%adv%compute_adjoint_scalar(u, v, w, s_adj, f_xh, &
395 xh, this%c_Xh, dm_xh%size())
402 call makeext%compute_scalar(this%abx1, this%abx2, f_xh%x, &
403 rho%x(1,1,1,1), ext_bdf%advection_coeffs, n)
406 call makebdf%compute_scalar(s_adj_lag, f_xh%x, s_adj, c_xh%B, &
407 rho%x(1,1,1,1), dt, ext_bdf%diffusion_coeffs, ext_bdf%ndiff, n)
410 call s_adj_lag%update()
413 call this%bcs%apply_scalar(this%s_adj%x, this%dm_Xh%size(), time, &
417 call this%update_material_properties(time)
420 call profiler_start_region(
'Adjoint_scalar_residual', 20)
421 call res%compute(ax, s_adj, s_adj_res, f_xh, c_xh, msh, xh, &
422 lambda, rho%x(1,1,1,1)*cp%x(1,1,1,1), ext_bdf%diffusion_coeffs(1), &
425 call gs_xh%op(s_adj_res, gs_op_add)
429 call this%bclst_ds%apply_scalar(s_adj_res%x, dm_xh%size())
431 call profiler_end_region(
'Adjoint_scalar_residual', 20)
433 call this%proj_s%pre_solving(s_adj_res%x, tstep, c_xh, n, dt_controller)
435 call this%pc%update()
436 call profiler_start_region(
'Adjoint_scalar_solve', 21)
437 ksp_results(1) = this%ksp%solve(ax, ds_adj, s_adj_res%x, n, &
438 c_xh, this%bclst_ds, gs_xh)
439 call profiler_end_region(
'Adjoint_scalar_solve', 21)
441 ksp_results(1)%name =
'Adjoint Scalar'
443 call this%proj_s%post_solving(ds_adj%x, ax, c_xh, this%bclst_ds, gs_xh, &
444 n, tstep, dt_controller)
447 if (neko_bcknd_device .eq. 1)
then
448 call device_add2s2(s_adj%x_d, ds_adj%x_d, 1.0_rp, n)
450 call add2s2(s_adj%x, ds_adj%x, 1.0_rp, n)
453 call scalar_step_info(time, ksp_results)
456 call profiler_end_region(
'Scalar', 2)
457 end subroutine adjoint_scalar_pnpn_step
459 subroutine print_debug(this)
460 class(adjoint_scalar_pnpn_t),
intent(inout) :: this
464 n = this%dm_Xh%size()
475 end subroutine print_debug
480 subroutine adjoint_scalar_pnpn_setup_bcs_(this, user)
481 class(adjoint_scalar_pnpn_t),
intent(inout) :: this
482 type(user_t),
target,
intent(in) :: user
483 integer :: i, j, n_bcs, zone_size, global_zone_size, ierr
484 type(json_core) :: core
485 type(json_value),
pointer :: bc_object
486 type(json_file) :: bc_subdict
487 class(bc_t),
pointer :: bc_i
490 logical,
allocatable :: marked_zones(:)
491 integer,
allocatable :: zone_indices(:)
493 if (this%params%valid_path(
'boundary_conditions'))
then
494 call this%params%info(
'boundary_conditions', &
496 call this%params%get_core(core)
497 call this%params%get(
'boundary_conditions', bc_object, found)
499 call this%bcs%init(n_bcs)
501 allocate(marked_zones(
size(this%msh%labeled_zones)))
502 marked_zones = .false.
506 call json_extract_item(core, bc_object, i, bc_subdict)
511 call json_get(bc_subdict,
"zone_indices", zone_indices)
513 do j = 1,
size(zone_indices)
514 zone_size = this%msh%labeled_zones(zone_indices(j))%size
515 call mpi_allreduce(zone_size, global_zone_size, 1, &
516 mpi_integer, mpi_max, neko_comm, ierr)
518 if (global_zone_size .eq. 0)
then
519 write(error_unit,
'(A, A, I0, A, A, I0, A)')
"*** ERROR ***: ",&
520 "Zone index ", zone_indices(j), &
521 " is invalid as this zone has 0 size, meaning it ", &
522 "does not exist in the mesh. Check adjoint scalar BC ", &
527 if (marked_zones(zone_indices(j)) .eqv. .true.)
then
528 write(error_unit,
'(A, A, I0, A, A, A, A)')
"*** ERROR ***: ", &
529 "Zone with index ", zone_indices(j), &
530 " has already been assigned a boundary condition. ", &
531 "Please check your boundary_conditions entry for the ", &
532 "adjoint scalar and make sure that each zone index ", &
533 "appears only in a single boundary condition."
536 marked_zones(zone_indices(j)) = .true.
542 call adjoint_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
543 call this%bcs%append(bc_i)
547 do i = 1,
size(this%msh%labeled_zones)
548 if ((this%msh%labeled_zones(i)%size .gt. 0) .and. &
549 (marked_zones(i) .eqv. .false.))
then
550 write(error_unit,
'(A, A, I0)')
"*** ERROR ***: ", &
551 "No adjoint scalar boundary condition assigned to zone ", i
557 do i = 1,
size(this%msh%labeled_zones)
558 if (this%msh%labeled_zones(i)%size .gt. 0)
then
559 write(error_unit,
'(A, A, A)')
"*** ERROR ***: ", &
560 "No boundary_conditions entry in the case file for " // &
561 " adjoint scalar ", &
567 end subroutine adjoint_scalar_pnpn_setup_bcs_
Boundary condition factory. Both constructs and initializes the object.
Contains the adjoint_scalar_pnpn_t type.
Contains the adjoint_scalar_scheme_t type.
Subroutines to add advection terms to the RHS of a transport equation.
Base type for a scalar advection-diffusion solver.
Base abstract type for computing the advection operator.