Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
adjoint_scalar_pnpn.f90
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
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
61 use advection_adjoint, only : advection_adjoint_t, advection_adjoint_factory
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
71 use bc, only : bc_t
72 use mpi_f08, only: mpi_integer, mpi_sum, mpi_max
73 implicit none
74 private
75
76
78
80 type(field_t) :: s_adj_res
81
83 type(field_t) :: ds_adj
84
86 class(ax_t), allocatable :: ax
87
89 type(projection_t) :: proj_s
90
94 type(zero_dirichlet_t) :: bc_res
95
100 type(bc_list_t) :: bclst_ds
101
103 class(advection_adjoint_t), allocatable :: adv
104
105 ! Time interpolation scheme
106 logical :: oifs
107
108 ! Advection terms for the oifs method
109 type(field_t) :: advs
110
112 class(scalar_residual_t), allocatable :: res
113
115 class(rhs_maker_ext_t), allocatable :: makeext
116
118 class(rhs_maker_bdf_t), allocatable :: makebdf
119
121 class(rhs_maker_oifs_t), allocatable :: makeoifs
122
123 contains
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_
134 end type adjoint_scalar_pnpn_t
135
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
151
152contains
153
170 subroutine adjoint_scalar_pnpn_init(this, msh, coef, gs, params_adjoint, &
171 params_primal, numerics_params, user, chkp, ulag, vlag, wlag, &
172 time_scheme, rho)
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
185 integer :: i
186 class(bc_t), pointer :: bc_i
187 character(len=15), parameter :: scheme = 'Modular (Pn/Pn)'
188 logical :: advection
189
190 call this%free()
191
192 ! Initiliaze base type.
193 call this%scheme_init(msh, coef, gs, params_adjoint, params_primal, &
194 scheme, user, rho)
195
196 ! Setup backend dependent Ax routines
197 call ax_helm_factory(this%ax, full_formulation = .false.)
198
199 ! Setup backend dependent scalar residual routines
200 call scalar_residual_factory(this%res)
201
202 ! Setup backend dependent summation of extrapolation scheme
203 call rhs_maker_ext_fctry(this%makeext)
204
205 ! Setup backend dependent contributions to F from lagged BD terms
206 call rhs_maker_bdf_fctry(this%makebdf)
207
208 ! Setup backend dependent contributions of the OIFS scheme
209 call rhs_maker_oifs_fctry(this%makeoifs)
210
211 ! Initialize variables specific to this plan
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)
214
215 call this%s_adj_res%init(dm_xh, "s_adj_res")
216
217 call this%abx1%init(dm_xh, "abx1")
218
219 call this%abx2%init(dm_xh, "abx2")
220
221 call this%advs%init(dm_xh, "advs")
222
223 call this%ds_adj%init(dm_xh, 'ds_adj')
224
225 end associate
226
227 ! Set up boundary conditions
228 call this%setup_bcs_(user)
229
230 ! Initialize dirichlet bcs for scalar residual
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)
236 end if
237 end do
238
239! call this%bc_res%mark_zones_from_list('d_s', this%bc_labels)
240 call this%bc_res%finalize()
241
242 call this%bclst_ds%init()
243 call this%bclst_ds%append(this%bc_res)
244
245
246 ! Initialize projection space
247 call this%proj_s%init(this%dm_Xh%size(), this%projection_dim, &
248 this%projection_activ_step)
249
250 ! Determine the time-interpolation scheme
251 call json_get_or_default(numerics_params, 'oifs', this%oifs, .false.)
252
253 ! Point to case checkpoint
254 this%chkp => chkp
255
256 ! Initialize advection factory
257 call json_get_or_default(params_adjoint, 'advection', advection, .true.)
258
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, &
262 this%s_adj_lag)
263 ! Add lagged term to checkpoint
264 ! @todo Init chkp object, note, adding 3 slags
265
266 ! Add scalar info to checkpoint
267 ! call this%chkp%add_scalar(this%s)
268 ! this%chkp%abs1 => this%abx1
269 ! this%chkp%abs2 => this%abx2
270 ! this%chkp%slag => this%slag
271
272 end subroutine adjoint_scalar_pnpn_init
273
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)
279 integer :: n
280 dtlag = chkp%dtlag
281 tlag = chkp%tlag
282
283 n = this%s_adj%dof%size()
284
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.)
301 end if
302
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)
306
307 end subroutine adjoint_scalar_pnpn_restart
308
309 subroutine adjoint_scalar_pnpn_free(this)
310 class(adjoint_scalar_pnpn_t), intent(inout) :: this
311
312 !Deallocate scalar field
313 call this%scheme_free()
314
315 call this%bclst_ds%free()
316 call this%proj_s%free()
317
318 call this%s_adj_res%free()
319
320 call this%ds_adj%free()
321
322 call this%abx1%free()
323 call this%abx2%free()
324
325 call this%advs%free()
326
327 if (allocated(this%Ax)) then
328 deallocate(this%Ax)
329 end if
330
331 if (allocated(this%res)) then
332 deallocate(this%res)
333 end if
334
335 if (allocated(this%makeext)) then
336 deallocate(this%makeext)
337 end if
338
339 if (allocated(this%makebdf)) then
340 deallocate(this%makebdf)
341 end if
342
343 if (allocated(this%makeoifs)) then
344 deallocate(this%makeoifs)
345 end if
346
347 end subroutine adjoint_scalar_pnpn_free
348
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
354 ! Number of degrees of freedom
355 integer :: n
356 ! Linear solver results monitor
357 type(ksp_monitor_t) :: ksp_results(1)
358
359 n = this%dm_Xh%size()
360
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)
373
374 ! Logs extra information the log level is NEKO_LOG_DEBUG or above.
375 call print_debug(this)
376 ! Compute the source terms
377 call this%source_term%compute(time)
378
379 ! Apply weak boundary conditions, that contribute to the source terms.
380 call this%bcs%apply_scalar(this%f_Xh%x, dm_xh%size(), time, .false.)
381
382 ! if (oifs) then
383 ! call neko_error("oifs not implemented for adjoint scalar")
384 ! ! Add the advection operators to the right-hans-side.
385 ! call this%adv%compute_scalar(u, v, w, s_adj, this%advs, &
386 ! Xh, this%c_Xh, dm_Xh%size())
387
388 ! call makeext%compute_scalar(this%abx1, this%abx2, f_Xh%x, rho, &
389 ! ext_bdf%advection_coeffs, n)
390
391 ! call makeoifs%compute_scalar(this%advs%x, f_Xh%x, rho, dt, n)
392 ! else
393 ! Add the advection operators to the right-hans-side.
394 call this%adv%compute_adjoint_scalar(u, v, w, s_adj, f_xh, &
395 xh, this%c_Xh, dm_xh%size())
396
397 ! At this point the RHS contains the sum of the advection operator,
398 ! Neumann boundary sources and additional source terms, evaluated using
399 ! the scalar field from the previous time-step.
400 ! Now, this value is used in the explicit time scheme to advance these
401 ! terms in time.
402 call makeext%compute_scalar(this%abx1, this%abx2, f_xh%x, &
403 rho%x(1,1,1,1), ext_bdf%advection_coeffs, n)
404
405 ! Add the RHS contributions coming from the BDF scheme.
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)
408 ! end if
409
410 call s_adj_lag%update()
411
413 call this%bcs%apply_scalar(this%s_adj%x, this%dm_Xh%size(), time, &
414 .true.)
415
416 ! Update material properties if necessary
417 call this%update_material_properties(time)
418
419 ! Compute scalar residual.
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), &
423 dt, dm_xh%size())
424
425 call gs_xh%op(s_adj_res, gs_op_add)
426
427
428 ! Apply a 0-valued Dirichlet boundary conditions on the ds_adj.
429 call this%bclst_ds%apply_scalar(s_adj_res%x, dm_xh%size())
430
431 call profiler_end_region('Adjoint_scalar_residual', 20)
432
433 call this%proj_s%pre_solving(s_adj_res%x, tstep, c_xh, n, dt_controller)
434
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)
440
441 ksp_results(1)%name = 'Adjoint Scalar'
442
443 call this%proj_s%post_solving(ds_adj%x, ax, c_xh, this%bclst_ds, gs_xh, &
444 n, tstep, dt_controller)
445
446 ! Update the solution
447 if (neko_bcknd_device .eq. 1) then
448 call device_add2s2(s_adj%x_d, ds_adj%x_d, 1.0_rp, n)
449 else
450 call add2s2(s_adj%x, ds_adj%x, 1.0_rp, n)
451 end if
452
453 call scalar_step_info(time, ksp_results)
454
455 end associate
456 call profiler_end_region('Scalar', 2)
457 end subroutine adjoint_scalar_pnpn_step
458
459 subroutine print_debug(this)
460 class(adjoint_scalar_pnpn_t), intent(inout) :: this
461 ! character(len=LOG_SIZE) :: log_buf
462 integer :: n
463
464 n = this%dm_Xh%size()
465 ! TODO come back to this
466 !write(log_buf,'(A,A,E15.7,A,E15.7,A,E15.7)') 'Adjoint scalar debug', &
467 ! ' l2norm s_adj', glsc2(this%s_adj%x, this%s_adj%x, n), &
468 ! ' slag1', glsc2(this%s_adj_lag%lf(1)%x, this%s_adj_lag%lf(1)%x, n), &
469 ! ' slag2', glsc2(this%s_adj_lag%lf(2)%x, this%s_adj_lag%lf(2)%x, n)
470 !call neko_log%message(log_buf, lvl=NEKO_LOG_DEBUG)
471 !write(log_buf,'(A,A,E15.7,A,E15.7)') 'Adjoint scalar debug2', &
472 ! ' l2norm abx1', glsc2(this%abx1%x, this%abx1%x, n), &
473 ! ' abx2', glsc2(this%abx2%x, this%abx2%x, n)
474 !call neko_log%message(log_buf, lvl=NEKO_LOG_DEBUG)
475 end subroutine print_debug
476
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
488 logical :: found
489 ! Monitor which boundary zones have been marked
490 logical, allocatable :: marked_zones(:)
491 integer, allocatable :: zone_indices(:)
492
493 if (this%params%valid_path('boundary_conditions')) then
494 call this%params%info('boundary_conditions', &
495 n_children = n_bcs)
496 call this%params%get_core(core)
497 call this%params%get('boundary_conditions', bc_object, found)
498
499 call this%bcs%init(n_bcs)
500
501 allocate(marked_zones(size(this%msh%labeled_zones)))
502 marked_zones = .false.
503
504 do i = 1, n_bcs
505 ! Create a new json containing just the subdict for this bc
506 call json_extract_item(core, bc_object, i, bc_subdict)
507
508 ! Check that we are not trying to assing a bc to zone, for which one
509 ! has already been assigned and that the zone has more than 0 size
510 ! in the mesh.
511 call json_get(bc_subdict, "zone_indices", zone_indices)
512
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)
517
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 ", &
523 i, "."
524 error stop
525 end if
526
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."
534 error stop
535 else
536 marked_zones(zone_indices(j)) = .true.
537 end if
538 end do
539
540 bc_i => null()
541
542 call adjoint_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
543 call this%bcs%append(bc_i)
544 end do
545
546 ! Make sure all labeled zones with non-zero size have been marked
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
552 error stop
553 end if
554 end do
555 else
556 ! Check that there are no labeled zones, i.e. all are periodic.
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 ", &
562 this%s%name
563 error stop
564 end if
565 end do
566 end if
567 end subroutine adjoint_scalar_pnpn_setup_bcs_
568
569end module adjoint_scalar_pnpn
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.