Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
adjoint_fluid_pnpn.f90
Go to the documentation of this file.
1! Copyright (c) 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!
34module adjoint_fluid_pnpn
35 use, intrinsic :: iso_fortran_env, only: error_unit
36 use coefs, only: coef_t
37 use symmetry, only: symmetry_t
38 use field_registry, only: neko_field_registry
39 use logger, only: neko_log, log_size, neko_log_debug
40 use num_types, only: rp
41 use krylov, only: ksp_monitor_t
42 use pnpn_residual, only: pnpn_prs_res_t, pnpn_vel_res_t, &
43 pnpn_prs_res_factory, pnpn_vel_res_factory, &
44 pnpn_prs_res_stress_factory, pnpn_vel_res_stress_factory
45 use rhs_maker, only: rhs_maker_sumab_t, rhs_maker_bdf_t, rhs_maker_ext_t, &
46 rhs_maker_oifs_t, rhs_maker_sumab_fctry, rhs_maker_bdf_fctry, &
47 rhs_maker_ext_fctry, rhs_maker_oifs_fctry
48 ! use adjoint_volflow, only: adjoint_volflow_t
49 use adjoint_fluid_scheme, only: adjoint_fluid_scheme_t
50 use device_mathops, only: device_opcolv, device_opadd2cm
51 use fluid_aux, only: fluid_step_info
52 use time_scheme_controller, only: time_scheme_controller_t
53 use projection, only: projection_t
54 use device, only : device_memcpy, host_to_device, device_event_sync,&
55 device_event_create, device_event_destroy, device_stream_wait_event, &
56 glb_cmd_queue
57 use advection_adjoint, only: advection_adjoint_t
58 use advection_adjoint_fctry, only: advection_adjoint_factory
59 use profiler, only: profiler_start_region, profiler_end_region
60 use json_module, only: json_file, json_core, json_value
61 use json_utils, only: json_get, json_get_or_default, json_extract_item
62 use json_module, only: json_file
63 use ax_product, only: ax_t, ax_helm_factory
64 use field, only: field_t
65 use dirichlet, only: dirichlet_t
66 use shear_stress, only: shear_stress_t
67 use wall_model_bc, only: wall_model_bc_t
68 use facet_normal, only: facet_normal_t
69 use non_normal, only: non_normal_t
70 use mesh, only: mesh_t
71 use user_intf, only: user_t
72 use time_step_controller, only: time_step_controller_t
73 use gs_ops, only: gs_op_add
74 use neko_config, only: neko_bcknd_device
75 use mathops, only: opadd2cm, opcolv
76 use bc_list, only: bc_list_t
77 use zero_dirichlet, only: zero_dirichlet_t
78 use utils, only: neko_error, neko_type_error
79 use field_math, only: field_add2, field_copy
80 use bc, only: bc_t
81 use file, only: file_t
82 use operators, only: ortho
83 use vector, only: vector_t
84 use device_math, only: device_vlsc3, device_cmult
85 use math, only: vlsc3, cmult
86 use json_utils_ext, only: json_key_fallback
87 use, intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr, c_associated
88 use comm, only: neko_comm, mpi_real_precision
89 use mpi_f08, only: mpi_sum, mpi_max, mpi_allreduce, mpi_comm_world, &
90 mpi_integer, mpi_logical, mpi_lor
91
92 implicit none
93 private
94
95
96 type, public, extends(adjoint_fluid_scheme_t) :: adjoint_fluid_pnpn_t
97
99 type(field_t) :: p_res, u_res, v_res, w_res
100
103 type(field_t) :: dp, du, dv, dw
104
105 type(field_t), pointer :: u_b => null()
106 type(field_t), pointer :: v_b => null()
107 type(field_t), pointer :: w_b => null()
108 type(field_t), pointer :: p_b => null()
109
110 !
111 ! Implicit operators, i.e. the left-hand-side of the Helmholz problem.
112 !
113
114 ! Coupled Helmholz operator for velocity
115 class(ax_t), allocatable :: ax_vel
116 ! Helmholz operator for pressure
117 class(ax_t), allocatable :: ax_prs
118
119 !
120 ! Projections for solver speed-up
121 !
122
124 type(projection_t) :: proj_prs
126 type(projection_t) :: proj_u
128 type(projection_t) :: proj_v
130 type(projection_t) :: proj_w
131
132 !
133 ! Special Karniadakis scheme boundary conditions in the pressure equation
134 !
135
137 type(facet_normal_t) :: bc_prs_surface
138
140 type(facet_normal_t) :: bc_sym_surface
141
142 !
143 ! Boundary conditions and lists for residuals and solution increments
144 !
145
147 type(zero_dirichlet_t) :: bc_vel_res
149 type(zero_dirichlet_t) :: bc_du
151 type(zero_dirichlet_t) :: bc_dv
153 type(zero_dirichlet_t) :: bc_dw
155 type(zero_dirichlet_t) :: bc_dp
156
158 type(bc_list_t) :: bclst_vel_res
159 type(bc_list_t) :: bclst_du
160 type(bc_list_t) :: bclst_dv
161 type(bc_list_t) :: bclst_dw
162 type(bc_list_t) :: bclst_dp
163
164
165 ! Checker for wether we have a strong pressure bc. If not, the pressure
166 ! is demeaned at every time step.
167 logical :: prs_dirichlet = .false.
168
169
170 ! The advection operator.
171 class(advection_adjoint_t), allocatable :: adv
172
173 ! Time OIFS interpolation scheme for advection.
174 logical :: oifs
175
176 ! Time variables
177 type(field_t) :: abx1, aby1, abz1
178 type(field_t) :: abx2, aby2, abz2
179
180 ! Advection terms for the oifs method
181 type(field_t) :: advx, advy, advz
182
184 class(pnpn_prs_res_t), allocatable :: prs_res
185
187 class(pnpn_vel_res_t), allocatable :: vel_res
188
190 class(rhs_maker_sumab_t), allocatable :: sumab
191
193 class(rhs_maker_ext_t), allocatable :: makeabf
194
196 class(rhs_maker_bdf_t), allocatable :: makebdf
197
199 class(rhs_maker_oifs_t), allocatable :: makeoifs
200
202 ! type(fluid_volflow_t) :: vol_flow
204 type(c_ptr) :: event = c_null_ptr
205
206 ! ======================================================================= !
207 ! Addressable attributes
208
209 real(kind=rp) :: norm_scaling
210 real(kind=rp) :: norm_target
211 real(kind=rp) :: norm_tolerance
212
213 ! ======================================================================= !
214 ! Definition of shorthands and local variables
215
217 real(kind=rp) :: norm_l2_base
218
220 real(kind=rp) :: norm_l2_upper
222 real(kind=rp) :: norm_l2_lower
223
225 type(file_t) :: file_output
226
227 contains
229 procedure, pass(this) :: init => adjoint_fluid_pnpn_init
231 procedure, pass(this) :: free => adjoint_fluid_pnpn_free
233 procedure, pass(this) :: step => adjoint_fluid_pnpn_step
235 procedure, pass(this) :: restart => adjoint_fluid_pnpn_restart
237 procedure, pass(this) :: setup_bcs => adjoint_fluid_pnpn_setup_bcs
239 procedure, pass(this) :: write_boundary_conditions => &
240 adjoint_fluid_pnpn_write_boundary_conditions
241
243 procedure, public, pass(this) :: pw_compute_ => power_iterations_compute
244
245 end type adjoint_fluid_pnpn_t
246
247 interface
248
255 module subroutine pressure_bc_factory(object, scheme, json, coef, user)
256 class(bc_t), pointer, intent(inout) :: object
257 type(adjoint_fluid_pnpn_t), intent(in) :: scheme
258 type(json_file), intent(inout) :: json
259 type(coef_t), intent(in) :: coef
260 type(user_t), intent(in) :: user
261 end subroutine pressure_bc_factory
262 end interface
263
264 interface
265
272 module subroutine velocity_bc_factory(object, scheme, json, coef, user)
273 class(bc_t), pointer, intent(inout) :: object
274 type(adjoint_fluid_pnpn_t), intent(in) :: scheme
275 type(json_file), intent(inout) :: json
276 type(coef_t), intent(in) :: coef
277 type(user_t), intent(in) :: user
278 end subroutine velocity_bc_factory
279 end interface
280
281contains
282
283 subroutine adjoint_fluid_pnpn_init(this, msh, lx, params, user, time_scheme)
284 class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
285 type(mesh_t), target, intent(inout) :: msh
286 integer, intent(in) :: lx
287 type(json_file), target, intent(inout) :: params
288 type(user_t), target, intent(in) :: user
289 type(time_scheme_controller_t), target, intent(in) :: time_scheme
290 character(len=15), parameter :: scheme = 'Adjoint (Pn/Pn)'
291 real(kind=rp) :: abs_tol
292 character(len=LOG_SIZE) :: log_buf
293 integer :: solver_maxiter
294 character(len=:), allocatable :: solver_type, precon_type
295 logical :: monitor, found
296 logical :: advection
297
298 ! Temporary field pointers
299 character(len=:), allocatable :: file_name
300 character(len=256) :: header_line
301
302 call this%free()
303
304 ! Initialize base class
305 call this%init_base(msh, lx, params, scheme, user, .true.)
306
307 ! Add pressure field to the registry. For this scheme it is in the same
308 ! Xh as the velocity
309 call neko_field_registry%add_field(this%dm_Xh, 'p_adj')
310 this%p_adj => neko_field_registry%get_field('p_adj')
311
312 !
313 ! Select governing equations via associated residual and Ax types
314 !
315
316 if (this%variable_material_properties .eqv. .true.) then
317 ! Setup backend dependent Ax routines
318 call ax_helm_factory(this%Ax_vel, full_formulation = .true.)
319
320 ! Setup backend dependent prs residual routines
321 call pnpn_prs_res_stress_factory(this%prs_res)
322
323 ! Setup backend dependent vel residual routines
324 call pnpn_vel_res_stress_factory(this%vel_res)
325 else
326 ! Setup backend dependent Ax routines
327 call ax_helm_factory(this%Ax_vel, full_formulation = .false.)
328
329 ! Setup backend dependent prs residual routines
330 call pnpn_prs_res_factory(this%prs_res)
331
332 ! Setup backend dependent vel residual routines
333 call pnpn_vel_res_factory(this%vel_res)
334 end if
335
336 ! Setup Ax for the pressure
337 call ax_helm_factory(this%Ax_prs, full_formulation = .false.)
338
339
340 ! Setup backend dependent summation of AB/BDF
341 call rhs_maker_sumab_fctry(this%sumab)
342
343 ! Setup backend dependent summation of extrapolation scheme
344 call rhs_maker_ext_fctry(this%makeabf)
345
346 ! Setup backend depenent contributions to F from lagged BD terms
347 call rhs_maker_bdf_fctry(this%makebdf)
348
349 ! Setup backend dependent summations of the OIFS method
350 call rhs_maker_oifs_fctry(this%makeoifs)
351
352 ! Initialize variables specific to this plan
353 associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
354 dm_xh => this%dm_Xh, nelv => this%msh%nelv)
355
356 call this%p_res%init(dm_xh, "p_res")
357 call this%u_res%init(dm_xh, "u_res")
358 call this%v_res%init(dm_xh, "v_res")
359 call this%w_res%init(dm_xh, "w_res")
360 call this%abx1%init(dm_xh, "abx1")
361 call this%aby1%init(dm_xh, "aby1")
362 call this%abz1%init(dm_xh, "abz1")
363 call this%abx2%init(dm_xh, "abx2")
364 call this%aby2%init(dm_xh, "aby2")
365 call this%abz2%init(dm_xh, "abz2")
366 call this%advx%init(dm_xh, "advx")
367 call this%advy%init(dm_xh, "advy")
368 call this%advz%init(dm_xh, "advz")
369 end associate
370
371 call this%du%init(this%dm_Xh, 'du')
372 call this%dv%init(this%dm_Xh, 'dv')
373 call this%dw%init(this%dm_Xh, 'dw')
374 call this%dp%init(this%dm_Xh, 'dp')
375
376 ! Set up boundary conditions
377 call this%setup_bcs(user, params)
378
379 ! Check if we need to output boundaries
380 call json_get_or_default(params, 'case.output_boundary', found, .false.)
381 if (found) call this%write_boundary_conditions()
382
383 ! Intialize projection space
384 if (this%variable_material_properties .and. &
385 this%vel_projection_dim .gt. 0) then
386 call neko_error("Velocity projection not available for full stress &
387 &formulation")
388 end if
389
390
391 call this%proj_prs%init(this%dm_Xh%size(), this%pr_projection_dim, &
392 this%pr_projection_activ_step)
393
394 call this%proj_u%init(this%dm_Xh%size(), this%vel_projection_dim, &
395 this%vel_projection_activ_step)
396 call this%proj_v%init(this%dm_Xh%size(), this%vel_projection_dim, &
397 this%vel_projection_activ_step)
398 call this%proj_w%init(this%dm_Xh%size(), this%vel_projection_dim, &
399 this%vel_projection_activ_step)
400
401
402 ! Add lagged term to checkpoint
403 call this%chkp%add_lag(this%ulag, this%vlag, this%wlag)
404
405 ! Determine the time-interpolation scheme
406 call json_get_or_default(params, 'case.numerics.oifs', this%oifs, .false.)
407
408 ! Initialize the advection factory
409 call json_get_or_default(params, 'case.fluid.advection', advection, .true.)
410
411 ! Todo: make sure this actually follow the forward advection factory
412 call advection_adjoint_factory(this%adv, params, this%c_Xh)
413
414 if (params%valid_path('case.fluid.flow_rate_force')) then
415 call neko_error("Flow rate forcing not available for adjoint_fluid_pnpn")
416 ! call this%vol_flow%init(this%dm_Xh, params)
417 end if
418
419 ! Setup pressure solver
420 call neko_log%section("Pressure solver")
421
422 call json_get_or_default(params, &
423 'case.fluid.pressure_solver.max_iterations', &
424 solver_maxiter, 800)
425 call json_get(params, 'case.fluid.pressure_solver.type', solver_type)
426 call json_get(params, 'case.fluid.pressure_solver.preconditioner', &
427 precon_type)
428 call json_get(params, 'case.fluid.pressure_solver.absolute_tolerance', &
429 abs_tol)
430 call json_get_or_default(params, 'case.fluid.velocity_solver.monitor', &
431 monitor, .false.)
432 call neko_log%message('Type : ('// trim(solver_type) // &
433 ', ' // trim(precon_type) // ')')
434 write(log_buf, '(A,ES13.6)') 'Abs tol :', abs_tol
435 call neko_log%message(log_buf)
436
437 call this%solver_factory(this%ksp_prs, this%dm_Xh%size(), &
438 solver_type, solver_maxiter, abs_tol, monitor)
439 call this%precon_factory_(this%pc_prs, this%ksp_prs, &
440 this%c_Xh, this%dm_Xh, this%gs_Xh, this%bcs_prs, precon_type)
441
442 call neko_log%end_section()
443
444 ! Setup gather-scatter event (only relevant on devices)
445 if (neko_bcknd_device .eq. 1) then
446 call device_event_create(this%event, 2)
447 end if
448
449 ! ------------------------------------------------------------------------ !
450 ! Handling the rescaling and baseflow
451
452 ! Read the norm scaling from the json file
453 call json_get_or_default(params, 'norm_scaling', &
454 this%norm_scaling, 0.5_rp)
455
456 ! The baseflow is the solution to the forward.
457 ! Userdefined baseflows can be invoked via setting initial conditions
458 ! call neko_field_registry%add_field(this%dm_Xh, 'u')
459 ! call neko_field_registry%add_field(this%dm_Xh, 'v')
460 ! call neko_field_registry%add_field(this%dm_Xh, 'w')
461 ! call neko_field_registry%add_field(this%dm_Xh, 'p')
462 this%u_b => neko_field_registry%get_field('u')
463 this%v_b => neko_field_registry%get_field('v')
464 this%w_b => neko_field_registry%get_field('w')
465 this%p_b => neko_field_registry%get_field('p')
466
467 ! Read the json file
468 call json_get_or_default(params, 'norm_target', &
469 this%norm_target, -1.0_rp)
470 call json_get_or_default(params, 'norm_tolerance', &
471 this%norm_tolerance, 10.0_rp)
472
473 ! Build the header
474 call json_get_or_default(params, 'output_file', &
475 file_name, 'power_iterations.csv')
476 this%file_output = file_t(trim(file_name))
477 write(header_line, '(A)') 'Time, Norm, Scaling'
478 call this%file_output%set_header(header_line)
479
480 end subroutine adjoint_fluid_pnpn_init
481
482 subroutine adjoint_fluid_pnpn_restart(this, dtlag, tlag)
483 class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
484 real(kind=rp) :: dtlag(10), tlag(10)
485 integer :: i, j, n
486
487 n = this%u_adj%dof%size()
488 if (allocated(this%chkp%previous_mesh%elements) .or. &
489 this%chkp%previous_Xh%lx .ne. this%Xh%lx) then
490 associate(u => this%u_adj, v => this%v_adj, w => this%w_adj, &
491 p => this%p_adj, c_xh => this%c_Xh, ulag => this%ulag, &
492 vlag => this%vlag, wlag => this%wlag)
493 do concurrent(j = 1:n)
494 u%x(j,1,1,1) = u%x(j,1,1,1) * c_xh%mult(j,1,1,1)
495 v%x(j,1,1,1) = v%x(j,1,1,1) * c_xh%mult(j,1,1,1)
496 w%x(j,1,1,1) = w%x(j,1,1,1) * c_xh%mult(j,1,1,1)
497 p%x(j,1,1,1) = p%x(j,1,1,1) * c_xh%mult(j,1,1,1)
498 end do
499 do i = 1, this%ulag%size()
500 do concurrent(j = 1:n)
501 ulag%lf(i)%x(j,1,1,1) = ulag%lf(i)%x(j,1,1,1) &
502 * c_xh%mult(j,1,1,1)
503 vlag%lf(i)%x(j,1,1,1) = vlag%lf(i)%x(j,1,1,1) &
504 * c_xh%mult(j,1,1,1)
505 wlag%lf(i)%x(j,1,1,1) = wlag%lf(i)%x(j,1,1,1) &
506 * c_xh%mult(j,1,1,1)
507 end do
508 end do
509 end associate
510 end if
511
512 if (neko_bcknd_device .eq. 1) then
513 associate(u => this%u_adj, v => this%v_adj, w => this%w_adj, &
514 ulag => this%ulag, vlag => this%vlag, wlag => this%wlag,&
515 p => this%p_adj)
516 call device_memcpy(u%x, u%x_d, u%dof%size(), &
517 host_to_device, sync = .false.)
518 call device_memcpy(v%x, v%x_d, v%dof%size(), &
519 host_to_device, sync = .false.)
520 call device_memcpy(w%x, w%x_d, w%dof%size(), &
521 host_to_device, sync = .false.)
522 call device_memcpy(p%x, p%x_d, p%dof%size(), &
523 host_to_device, sync = .false.)
524 call device_memcpy(ulag%lf(1)%x, ulag%lf(1)%x_d, &
525 u%dof%size(), host_to_device, sync = .false.)
526 call device_memcpy(ulag%lf(2)%x, ulag%lf(2)%x_d, &
527 u%dof%size(), host_to_device, sync = .false.)
528
529 call device_memcpy(vlag%lf(1)%x, vlag%lf(1)%x_d, &
530 v%dof%size(), host_to_device, sync = .false.)
531 call device_memcpy(vlag%lf(2)%x, vlag%lf(2)%x_d, &
532 v%dof%size(), host_to_device, sync = .false.)
533
534 call device_memcpy(wlag%lf(1)%x, wlag%lf(1)%x_d, &
535 w%dof%size(), host_to_device, sync = .false.)
536 call device_memcpy(wlag%lf(2)%x, wlag%lf(2)%x_d, &
537 w%dof%size(), host_to_device, sync = .false.)
538 call device_memcpy(this%abx1%x, this%abx1%x_d, &
539 w%dof%size(), host_to_device, sync = .false.)
540 call device_memcpy(this%abx2%x, this%abx2%x_d, &
541 w%dof%size(), host_to_device, sync = .false.)
542 call device_memcpy(this%aby1%x, this%aby1%x_d, &
543 w%dof%size(), host_to_device, sync = .false.)
544 call device_memcpy(this%aby2%x, this%aby2%x_d, &
545 w%dof%size(), host_to_device, sync = .false.)
546 call device_memcpy(this%abz1%x, this%abz1%x_d, &
547 w%dof%size(), host_to_device, sync = .false.)
548 call device_memcpy(this%abz2%x, this%abz2%x_d, &
549 w%dof%size(), host_to_device, sync = .false.)
550 call device_memcpy(this%advx%x, this%advx%x_d, &
551 w%dof%size(), host_to_device, sync = .false.)
552 call device_memcpy(this%advy%x, this%advy%x_d, &
553 w%dof%size(), host_to_device, sync = .false.)
554 call device_memcpy(this%advz%x, this%advz%x_d, &
555 w%dof%size(), host_to_device, sync = .false.)
556 end associate
557 end if
558
559 ! Make sure that continuity is maintained (important for interpolation)
560 ! Do not do this for lagged rhs
561 ! (derivatives are not necessairly coninous across elements)
562
563 if (allocated(this%chkp%previous_mesh%elements) &
564 .or. this%chkp%previous_Xh%lx .ne. this%Xh%lx) then
565 call this%gs_Xh%op(this%u_adj, gs_op_add)
566 call this%gs_Xh%op(this%v_adj, gs_op_add)
567 call this%gs_Xh%op(this%w_adj, gs_op_add)
568 call this%gs_Xh%op(this%p_adj, gs_op_add)
569
570 do i = 1, this%ulag%size()
571 call this%gs_Xh%op(this%ulag%lf(i), gs_op_add)
572 call this%gs_Xh%op(this%vlag%lf(i), gs_op_add)
573 call this%gs_Xh%op(this%wlag%lf(i), gs_op_add)
574 end do
575 end if
576
577 !! If we would decide to only restart from lagged fields instead of asving
578 !! abx1, aby1 etc.
579 !! Observe that one also needs to recompute the focing at the old time steps
580 !u_temp = this%ulag%lf(2)
581 !v_temp = this%vlag%lf(2)
582 !w_temp = this%wlag%lf(2)
583 !! Compute the source terms
584 !call this%source_term%compute(tlag(2), -1)
585 !
586 !! Pre-multiply the source terms with the mass matrix.
587 !if (NEKO_BCKND_DEVICE .eq. 1) then
588 ! call device_opcolv(this%f_adj_x%x_d, this%f_adj_y%x_d, this%f_adj_z%x_d,
589 ! this%c_Xh%B_d, this%msh%gdim, n)
590 !else
591 ! call opcolv(this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, this%c_Xh%B,
592 ! this%msh%gdim, n)
593 !end if
594
595 !! Add the advection operators to the right-hand-side.
596 !call this%adv%compute(u_temp, v_temp, w_temp, &
597 ! this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, &
598 ! this%Xh, this%c_Xh, this%dm_Xh%size())
599 !this%abx2 = this%f_x
600 !this%aby2 = this%f_y
601 !this%abz2 = this%f_z
602 !
603 !u_temp = this%ulag%lf(1)
604 !v_temp = this%vlag%lf(1)
605 !w_temp = this%wlag%lf(1)
606 !call this%source_term%compute(tlag(1), 0)
607
608 !! Pre-multiply the source terms with the mass matrix.
609 !if (NEKO_BCKND_DEVICE .eq. 1) then
610 ! call device_opcolv(this%f_adj_x%x_d, this%f_adj_y%x_d, this%f_adj_z%x_d,
611 ! this%c_Xh%B_d, this%msh%gdim, n)
612 !else
613 ! call opcolv(this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, this%c_Xh%B,
614 ! this%msh%gdim, n)
615 !end if
616
617 !! Pre-multiply the source terms with the mass matrix.
618 !if (NEKO_BCKND_DEVICE .eq. 1) then
619 ! call device_opcolv(this%f_adj_x%x_d, this%f_adj_y%x_d, this%f_adj_z%x_d,
620 ! this%c_Xh%B_d, this%msh%gdim, n)
621 !else
622 ! call opcolv(this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, this%c_Xh%B,
623 ! this%msh%gdim, n)
624 !end if
625
626 !call this%adv%compute(u_temp, v_temp, w_temp, &
627 ! this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, &
628 ! this%Xh, this%c_Xh, this%dm_Xh%size())
629 !this%abx1 = this%f_x
630 !this%aby1 = this%f_y
631 !this%abz1 = this%f_z
632
633 end subroutine adjoint_fluid_pnpn_restart
634
635 subroutine adjoint_fluid_pnpn_free(this)
636 class(adjoint_fluid_pnpn_t), intent(inout) :: this
637
638 !Deallocate velocity and pressure fields
639 call this%scheme_free()
640
641 call this%bc_prs_surface%free()
642 call this%bc_sym_surface%free()
643 call this%bclst_vel_res%free()
644 call this%bclst_dp%free()
645 call this%proj_prs%free()
646 call this%proj_u%free()
647 call this%proj_v%free()
648 call this%proj_w%free()
649
650 call this%p_res%free()
651 call this%u_res%free()
652 call this%v_res%free()
653 call this%w_res%free()
654
655 call this%du%free()
656 call this%dv%free()
657 call this%dw%free()
658 call this%dp%free()
659
660 call this%abx1%free()
661 call this%aby1%free()
662 call this%abz1%free()
663
664 call this%abx2%free()
665 call this%aby2%free()
666 call this%abz2%free()
667
668 call this%advx%free()
669 call this%advy%free()
670 call this%advz%free()
671
672 if (allocated(this%Ax_vel)) then
673 deallocate(this%Ax_vel)
674 end if
675
676 if (allocated(this%Ax_prs)) then
677 deallocate(this%Ax_prs)
678 end if
679
680 if (allocated(this%prs_res)) then
681 deallocate(this%prs_res)
682 end if
683
684 if (allocated(this%vel_res)) then
685 deallocate(this%vel_res)
686 end if
687
688 if (allocated(this%sumab)) then
689 deallocate(this%sumab)
690 end if
691
692 if (allocated(this%makeabf)) then
693 deallocate(this%makeabf)
694 end if
695
696 if (allocated(this%makebdf)) then
697 deallocate(this%makebdf)
698 end if
699
700 if (allocated(this%makeoifs)) then
701 deallocate(this%makeoifs)
702 end if
703
704 ! call this%vol_flow%free()
705 if (neko_bcknd_device .eq. 1) then
706 if (c_associated(this%event)) then
707 call device_event_destroy(this%event)
708 end if
709 end if
710
711 end subroutine adjoint_fluid_pnpn_free
712
719 subroutine adjoint_fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
720 class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
721 real(kind=rp), intent(in) :: t
722 integer, intent(in) :: tstep
723 real(kind=rp), intent(in) :: dt
724 type(time_scheme_controller_t), intent(in) :: ext_bdf
725 type(time_step_controller_t), intent(in) :: dt_controller
726 ! number of degrees of freedom
727 integer :: n
728 ! Solver results monitors (pressure + 3 velocity)
729 type(ksp_monitor_t) :: ksp_results(4)
730 ! Extrapolated velocity for the pressure residual
731 type(field_t), pointer :: u_e, v_e, w_e
732 ! Indices for tracking temporary fields
733 integer :: temp_indices(3)
734
735 if (this%freeze) return
736
737 n = this%dm_Xh%size()
738
739 call profiler_start_region('Adjoint', 1)
740 associate(u => this%u_adj, v => this%v_adj, w => this%w_adj, &
741 p => this%p_adj, &
742 du => this%du, dv => this%dv, dw => this%dw, dp => this%dp, &
743 u_b => this%u_b, v_b => this%v_b, w_b => this%w_b, &
744 u_res => this%u_res, v_res => this%v_res, w_res => this%w_res, &
745 p_res => this%p_res, ax_vel => this%Ax_vel, ax_prs => this%Ax_prs, &
746 xh => this%Xh, &
747 c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
748 ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
749 msh => this%msh, prs_res => this%prs_res, &
750 source_term => this%source_term, &
751 vel_res => this%vel_res, sumab => this%sumab, &
752 makeabf => this%makeabf, makebdf => this%makebdf, &
753 vel_projection_dim => this%vel_projection_dim, &
754 pr_projection_dim => this%pr_projection_dim, &
755 rho => this%rho, mu => this%mu, oifs => this%oifs, &
756 rho_field => this%rho_field, mu_field => this%mu_field, &
757 f_x => this%f_adj_x, f_y => this%f_adj_y, f_z => this%f_adj_z, &
758 if_variable_dt => dt_controller%if_variable_dt, &
759 dt_last_change => dt_controller%dt_last_change, &
760 event => this%event)
761
762 ! Get temporary arrays
763 call this%scratch%request_field(u_e, temp_indices(1))
764 call this%scratch%request_field(v_e, temp_indices(2))
765 call this%scratch%request_field(w_e, temp_indices(3))
766 call sumab%compute_fluid(u_e, v_e, w_e, u, v, w, &
767 ulag, vlag, wlag, ext_bdf%advection_coeffs, ext_bdf%nadv)
768
769 ! Compute the source terms
770 call this%source_term%compute(t, tstep)
771
772 ! Add Neumann bc contributions to the RHS
773 call this%bcs_vel%apply_vector(f_x%x, f_y%x, f_z%x, &
774 this%dm_Xh%size(), t, tstep, strong = .false.)
775
776 ! Compute the grandient jump penalty term
777 if (this%if_gradient_jump_penalty .eqv. .true.) then
778 call neko_error("Gradient jump penalty not implemented for adjoint")
779 ! call this%gradient_jump_penalty_u_adj%compute(u, v, w, u)
780 ! call this%gradient_jump_penalty_v_adj%compute(u, v, w, v)
781 ! call this%gradient_jump_penalty_w_adj%compute(u, v, w, w)
782 ! call this%gradient_jump_penalty_u_adj%perform(f_x)
783 ! call this%gradient_jump_penalty_v_adj%perform(f_y)
784 ! call this%gradient_jump_penalty_w_adj%perform(f_z)
785 end if
786
787 if (oifs) then
788 call neko_error("OIFS not implemented for adjoint")
789 ! ! Add the advection operators to the right-hand-side.
790 ! call this%adv%compute_adjoint(u, v, w, &
791 ! this%advx, this%advy, this%advz, &
792 ! Xh, this%c_Xh, dm_Xh%size(), dt)
793
794 ! At this point the RHS contains the sum of the advection operator and
795 ! additional source terms, evaluated using the velocity field from the
796 ! previous time-step. Now, this value is used in the explicit time
797 ! scheme to advance both terms in time.
798 ! call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
799 ! this%abx2, this%aby2, this%abz2, &
800 ! f_x%x, f_y%x, f_z%x, &
801 ! rho, ext_bdf%advection_coeffs, n)
802
803 ! Now, the source terms from the previous time step are added to the
804 ! RHS.
805 ! call makeoifs%compute_fluid(this%advx%x, this%advy%x, this%advz%x, &
806 ! f_x%x, f_y%x, f_z%x, &
807 ! rho, dt, n)
808 else
809 call this%adv%compute_adjoint(u, v, w, u_b, v_b, w_b, &
810 f_x, f_y, f_z, &
811 xh, this%c_Xh, dm_xh%size())
812
813 ! At this point the RHS contains the sum of the advection operator and
814 ! additional source terms, evaluated using the velocity field from the
815 ! previous time-step. Now, this value is used in the explicit time
816 ! scheme to advance both terms in time.
817 call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
818 this%abx2, this%aby2, this%abz2, &
819 f_x%x, f_y%x, f_z%x, &
820 rho, ext_bdf%advection_coeffs, n)
821
822 ! Add the RHS contributions coming from the BDF scheme.
823 call makebdf%compute_fluid(ulag, vlag, wlag, f_x%x, f_y%x, f_z%x, &
824 u, v, w, c_xh%B, rho, dt, &
825 ext_bdf%diffusion_coeffs, ext_bdf%ndiff, n)
826 end if
827
828 call ulag%update()
829 call vlag%update()
830 call wlag%update()
831
832 call this%bc_apply_vel(t, tstep, strong = .true.)
833 call this%bc_apply_prs(t, tstep)
834
835 ! Update material properties if necessary
836 call this%update_material_properties()
837
838 ! Compute pressure residual.
839 call profiler_start_region('Pressure_residual', 18)
840
841 call prs_res%compute(p, p_res,&
842 u, v, w, &
843 u_e, v_e, w_e, &
844 f_x, f_y, f_z, &
845 c_xh, gs_xh, &
846 this%bc_prs_surface, this%bc_sym_surface,&
847 ax_prs, ext_bdf%diffusion_coeffs(1), dt, &
848 mu_field, rho_field, event)
849
850 ! De-mean the pressure residual when no strong pressure boundaries present
851 if (.not. this%prs_dirichlet) call ortho(p_res%x, this%glb_n_points, n)
852
853 call gs_xh%op(p_res, gs_op_add, event)
854 if (neko_bcknd_device .eq. 1) then
855 call device_stream_wait_event(glb_cmd_queue, event, 0)
856 end if
857
858 ! Set the residual to zero at strong pressure boundaries.
859 call this%bclst_dp%apply_scalar(p_res%x, p%dof%size(), t, tstep)
860
861
862 call profiler_end_region('Pressure_residual', 18)
863
864
865 call this%proj_prs%pre_solving(p_res%x, tstep, c_xh, n, dt_controller, &
866 'Pressure')
867
868 call this%pc_prs%update()
869
870 call profiler_start_region('Pressure_solve', 3)
871
872 ! Solve for the pressure increment.
873 ksp_results(1) = this%ksp_prs%solve(ax_prs, dp, p_res%x, n, c_xh, &
874 this%bclst_dp, gs_xh)
875
876
877 call profiler_end_region('Pressure_solve', 3)
878
879 call this%proj_prs%post_solving(dp%x, ax_prs, c_xh, &
880 this%bclst_dp, gs_xh, n, tstep, dt_controller)
881
882 ! Update the pressure with the increment. Demean if necessary.
883 call field_add2(p, dp, n)
884 if (.not. this%prs_dirichlet) call ortho(p%x, this%glb_n_points, n)
885
886 ! Compute velocity residual.
887 call profiler_start_region('Velocity_residual', 19)
888 call vel_res%compute(ax_vel, u, v, w, &
889 u_res, v_res, w_res, &
890 p, &
891 f_x, f_y, f_z, &
892 c_xh, msh, xh, &
893 mu_field, rho_field, ext_bdf%diffusion_coeffs(1), &
894 dt, dm_xh%size())
895
896 call gs_xh%op(u_res, gs_op_add, event)
897 call gs_xh%op(v_res, gs_op_add, event)
898 call gs_xh%op(w_res, gs_op_add, event)
899 if (neko_bcknd_device .eq. 1) then
900 call device_stream_wait_event(glb_cmd_queue, event, 0)
901 end if
902
903 ! Set residual to zero at strong velocity boundaries.
904 call this%bclst_vel_res%apply(u_res, v_res, w_res, t, tstep)
905
906
907 call profiler_end_region('Velocity_residual', 19)
908
909 call this%proj_u%pre_solving(u_res%x, tstep, c_xh, n, dt_controller)
910 call this%proj_v%pre_solving(v_res%x, tstep, c_xh, n, dt_controller)
911 call this%proj_w%pre_solving(w_res%x, tstep, c_xh, n, dt_controller)
912
913 call this%pc_vel%update()
914
915 call profiler_start_region("Velocity_solve", 4)
916 ksp_results(2:4) = this%ksp_vel%solve_coupled(ax_vel, du, dv, dw, &
917 u_res%x, v_res%x, w_res%x, n, c_xh, &
918 this%bclst_du, this%bclst_dv, this%bclst_dw, gs_xh, &
919 this%ksp_vel%max_iter)
920 call profiler_end_region("Velocity_solve", 4)
921
922 call this%proj_u%post_solving(du%x, ax_vel, c_xh, &
923 this%bclst_du, gs_xh, n, tstep, dt_controller)
924 call this%proj_v%post_solving(dv%x, ax_vel, c_xh, &
925 this%bclst_dv, gs_xh, n, tstep, dt_controller)
926 call this%proj_w%post_solving(dw%x, ax_vel, c_xh, &
927 this%bclst_dw, gs_xh, n, tstep, dt_controller)
928
929 if (neko_bcknd_device .eq. 1) then
930 call device_opadd2cm(u%x_d, v%x_d, w%x_d, &
931 du%x_d, dv%x_d, dw%x_d, 1.0_rp, n, msh%gdim)
932 else
933 call opadd2cm(u%x, v%x, w%x, du%x, dv%x, dw%x, 1.0_rp, n, msh%gdim)
934 end if
935
936 if (this%forced_flow_rate) then
937 call neko_error('Forced flow rate is not implemented for the adjoint')
938 !call this%vol_flow%adjust( u, v, w, p, u_res, v_res, w_res, p_res, &
939 ! c_Xh, gs_Xh, ext_bdf, rho, mu, dt, &
940 ! this%bclst_dp, this%bclst_du, this%bclst_dv, &
941 ! this%bclst_dw, this%bclst_vel_res, Ax_vel, Ax_prs, this%ksp_prs,&
942 ! this%ksp_vel, this%pc_prs, this%pc_vel, this%ksp_prs%max_iter, &
943 ! this%ksp_vel%max_iter)
944 end if
945
946 call fluid_step_info(tstep, t, dt, ksp_results, this%strict_convergence)
947
948 call this%scratch%relinquish_field(temp_indices)
949
950 end associate
951 call profiler_end_region('Adjoint', 1)
952
953 ! Compute the norm of the field and determine if we should do a rescale.
954 ! TODO: HARRY
955 ! we need to discuss this rescaling in the context of topology optimization
956 ! vs stability analysis,
957 ! for now I'm commenting this out
958 !call this%PW_compute_(t, tstep)
959
960 end subroutine adjoint_fluid_pnpn_step
961
962 ! ========================================================================== !
963
966 subroutine adjoint_fluid_pnpn_setup_bcs(this, user, params)
967 use mpi_f08, only: mpi_in_place
968
969 class(adjoint_fluid_pnpn_t), intent(inout) :: this
970 type(user_t), target, intent(in) :: user
971 type(json_file), intent(inout) :: params
972 integer :: i, n_bcs, j, zone_size, global_zone_size, ierr
973 class(bc_t), pointer :: bc_i
974 type(json_core) :: core
975 type(json_value), pointer :: bc_object
976 type(json_file) :: bc_subdict
977 logical :: found
978 ! Monitor which boundary zones have been marked
979 logical, allocatable :: marked_zones(:)
980 integer, allocatable :: zone_indices(:)
981 character(len=:), allocatable :: json_key
982
983 ! Lists for the residuals and solution increments
984 call this%bclst_vel_res%init()
985 call this%bclst_du%init()
986 call this%bclst_dv%init()
987 call this%bclst_dw%init()
988 call this%bclst_dp%init()
989
990 call this%bc_vel_res%init_from_components(this%c_Xh)
991 call this%bc_du%init_from_components(this%c_Xh)
992 call this%bc_dv%init_from_components(this%c_Xh)
993 call this%bc_dw%init_from_components(this%c_Xh)
994 call this%bc_dp%init_from_components(this%c_Xh)
995
996 ! Special PnPn boundary conditions for pressure
997 call this%bc_prs_surface%init_from_components(this%c_Xh)
998 call this%bc_sym_surface%init_from_components(this%c_Xh)
999
1000 json_key = 'case.adjoint_fluid.boundary_conditions'
1001
1002 ! Populate bcs_vel and bcs_prs based on the case file
1003 if (params%valid_path(json_key)) then
1004 call params%info(json_key, n_children = n_bcs)
1005 call params%get_core(core)
1006 call params%get(json_key, bc_object, found)
1007
1008 !
1009 ! Velocity bcs
1010 !
1011 call this%bcs_vel%init(n_bcs)
1012
1013 allocate(marked_zones(size(this%msh%labeled_zones)))
1014 marked_zones = .false.
1015
1016 do i = 1, n_bcs
1017 ! Create a new json containing just the subdict for this bc
1018 call json_extract_item(core, bc_object, i, bc_subdict)
1019
1020 call json_get(bc_subdict, "zone_indices", zone_indices)
1021
1022 ! Check that we are not trying to assing a bc to zone, for which one
1023 ! has already been assigned and that the zone has more than 0 size
1024 ! in the mesh.
1025 do j = 1, size(zone_indices)
1026 zone_size = this%msh%labeled_zones(zone_indices(j))%size
1027 call mpi_allreduce(zone_size, global_zone_size, 1, &
1028 mpi_integer, mpi_max, neko_comm, ierr)
1029
1030 if (global_zone_size .eq. 0) then
1031 write(error_unit, '(A, A, I0, A, A, I0, A)') "*** ERROR ***: ",&
1032 "Zone index ", zone_indices(j), &
1033 " is invalid as this zone has 0 size, meaning it does ", &
1034 "not in the mesh. Check adjoint boundary condition ", &
1035 i, "."
1036 error stop
1037 end if
1038
1039 if (marked_zones(zone_indices(j)) .eqv. .true.) then
1040 write(error_unit, '(A, A, I0, A, A, A, A)') "*** ERROR ***: ", &
1041 "Zone with index ", zone_indices(j), &
1042 " has already been assigned a boundary condition. ", &
1043 "Please check your boundary_conditions entry for the ", &
1044 "adjoint and make sure that each zone index appears ", &
1045 "only in a single boundary condition."
1046 error stop
1047 else
1048 marked_zones(zone_indices(j)) = .true.
1049 end if
1050 end do
1051
1052 bc_i => null()
1053 call velocity_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
1054
1055 ! Not all bcs require an allocation for velocity in particular,
1056 ! so we check.
1057 if (associated(bc_i)) then
1058
1059 ! We need to treat mixed bcs separately because they are by
1060 ! convention marked weak and currently contain nested
1061 ! bcs, some of which are strong.
1062 select type (bc_i)
1063 type is (symmetry_t)
1064 ! Symmetry has 3 internal bcs, but only one actually contains
1065 ! markings.
1066 ! Symmetry's apply_scalar doesn't do anything, so we need to
1067 ! mark individual nested bcs to the du,dv,dw, whereas the
1068 ! vel_res can just get symmetry as a whole, because on this
1069 ! list we call apply_vector.
1070 ! Additionally we have to mark the special surface bc for p.
1071 call this%bclst_vel_res%append(bc_i)
1072 call this%bc_du%mark_facets(bc_i%bc_x%marked_facet)
1073 call this%bc_dv%mark_facets(bc_i%bc_y%marked_facet)
1074 call this%bc_dw%mark_facets(bc_i%bc_z%marked_facet)
1075
1076 call this%bcs_vel%append(bc_i)
1077
1078 call this%bc_sym_surface%mark_facets(bc_i%marked_facet)
1079 type is (non_normal_t)
1080 ! This is a bc for the residuals and increments, not the
1081 ! velocity itself. So, don't append to bcs_vel
1082 call this%bclst_vel_res%append(bc_i)
1083 call this%bc_du%mark_facets(bc_i%bc_x%marked_facet)
1084 call this%bc_dv%mark_facets(bc_i%bc_y%marked_facet)
1085 call this%bc_dw%mark_facets(bc_i%bc_z%marked_facet)
1086 type is (shear_stress_t)
1087 ! Same as symmetry
1088 call this%bclst_vel_res%append(bc_i%symmetry)
1089 call this%bclst_du%append(bc_i%symmetry%bc_x)
1090 call this%bclst_dv%append(bc_i%symmetry%bc_y)
1091 call this%bclst_dw%append(bc_i%symmetry%bc_z)
1092
1093 call this%bcs_vel%append(bc_i)
1094 type is (wall_model_bc_t)
1095 ! Same as symmetry
1096 call this%bclst_vel_res%append(bc_i%symmetry)
1097 call this%bclst_du%append(bc_i%symmetry%bc_x)
1098 call this%bclst_dv%append(bc_i%symmetry%bc_y)
1099 call this%bclst_dw%append(bc_i%symmetry%bc_z)
1100
1101 call this%bcs_vel%append(bc_i)
1102 class default
1103
1104 ! For the default case we use our dummy zero_dirichlet bcs to
1105 ! mark the same faces as in ordinary velocity dirichlet
1106 ! conditions.
1107 ! Additionally we mark the special PnPn pressure bc.
1108 if (bc_i%strong .eqv. .true.) then
1109 call this%bc_vel_res%mark_facets(bc_i%marked_facet)
1110 call this%bc_du%mark_facets(bc_i%marked_facet)
1111 call this%bc_dv%mark_facets(bc_i%marked_facet)
1112 call this%bc_dw%mark_facets(bc_i%marked_facet)
1113
1114 call this%bc_prs_surface%mark_facets(bc_i%marked_facet)
1115 end if
1116
1117 call this%bcs_vel%append(bc_i)
1118 end select
1119 end if
1120 end do
1121
1122 ! Make sure all labeled zones with non-zero size have been marked
1123 do i = 1, size(this%msh%labeled_zones)
1124 if ((this%msh%labeled_zones(i)%size .gt. 0) .and. &
1125 (marked_zones(i) .eqv. .false.)) then
1126 write(error_unit, '(A, A, I0)') "*** ERROR ***: ", &
1127 "No adjoint boundary condition assigned to zone ", i
1128 error stop
1129 end if
1130 end do
1131
1132 !
1133 ! Pressure bcs
1134 !
1135 call this%bcs_prs%init(n_bcs)
1136
1137 do i = 1, n_bcs
1138 ! Create a new json containing just the subdict for this bc
1139 call json_extract_item(core, bc_object, i, bc_subdict)
1140 bc_i => null()
1141 call pressure_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
1142
1143 ! Not all bcs require an allocation for pressure in particular,
1144 ! so we check.
1145 if (associated(bc_i)) then
1146 call this%bcs_prs%append(bc_i)
1147
1148 ! Mark strong bcs in the dummy dp bc to force zero change.
1149 if (bc_i%strong .eqv. .true.) then
1150 call this%bc_dp%mark_facets(bc_i%marked_facet)
1151 end if
1152
1153 end if
1154
1155 end do
1156 end if
1157
1158 call this%bc_prs_surface%finalize()
1159 call this%bc_sym_surface%finalize()
1160
1161 call this%bc_vel_res%finalize()
1162 call this%bc_du%finalize()
1163 call this%bc_dv%finalize()
1164 call this%bc_dw%finalize()
1165 call this%bc_dp%finalize()
1166
1167 call this%bclst_vel_res%append(this%bc_vel_res)
1168 call this%bclst_du%append(this%bc_du)
1169 call this%bclst_dv%append(this%bc_dv)
1170 call this%bclst_dw%append(this%bc_dw)
1171 call this%bclst_dp%append(this%bc_dp)
1172
1173 ! If we have no strong pressure bcs, we will demean the pressure
1174 this%prs_dirichlet = .not. this%bclst_dp%is_empty()
1175 call mpi_allreduce(mpi_in_place, this%prs_dirichlet, 1, &
1176 mpi_logical, mpi_lor, neko_comm)
1177
1178 end subroutine adjoint_fluid_pnpn_setup_bcs
1179
1181 subroutine adjoint_fluid_pnpn_write_boundary_conditions(this)
1182 use inflow, only: inflow_t
1183 use field_dirichlet, only: field_dirichlet_t
1184 use blasius, only: blasius_t
1185 use field_dirichlet_vector, only: field_dirichlet_vector_t
1186 use usr_inflow, only: usr_inflow_t
1187 use dong_outflow, only: dong_outflow_t
1188 class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
1189 type(dirichlet_t) :: bdry_mask
1190 type(field_t), pointer :: bdry_field
1191 type(file_t) :: bdry_file
1192 integer :: temp_index, i
1193 class(bc_t), pointer :: bci
1194 character(len=LOG_SIZE) :: log_buf
1195
1196 call neko_log%section("Adjoint boundary conditions")
1197 write(log_buf, '(A)') &
1198 'Marking using integer keys in boundary_adjoint0.f00000'
1199 call neko_log%message(log_buf)
1200 write(log_buf, '(A)') 'Condition-value pairs: '
1201 call neko_log%message(log_buf)
1202 write(log_buf, '(A)') ' no_slip = 1'
1203 call neko_log%message(log_buf)
1204 write(log_buf, '(A)') ' velocity_value = 2'
1205 call neko_log%message(log_buf)
1206 write(log_buf, '(A)') ' outflow, normal_outflow (+dong) = 3'
1207 call neko_log%message(log_buf)
1208 write(log_buf, '(A)') ' symmetry = 4'
1209 call neko_log%message(log_buf)
1210 write(log_buf, '(A)') ' user_velocity_pointwise = 5'
1211 call neko_log%message(log_buf)
1212 write(log_buf, '(A)') ' periodic = 6'
1213 call neko_log%message(log_buf)
1214 write(log_buf, '(A)') ' user_velocity = 7'
1215 call neko_log%message(log_buf)
1216 write(log_buf, '(A)') ' user_pressure = 8'
1217 call neko_log%message(log_buf)
1218 write(log_buf, '(A)') ' shear_stress = 9'
1219 call neko_log%message(log_buf)
1220 write(log_buf, '(A)') ' wall_modelling = 10'
1221 call neko_log%message(log_buf)
1222 write(log_buf, '(A)') ' blasius_profile = 11'
1223 call neko_log%message(log_buf)
1224 call neko_log%end_section()
1225
1226 call this%scratch%request_field(bdry_field, temp_index)
1227 bdry_field = 0.0_rp
1228
1229
1230
1231 call bdry_mask%init_from_components(this%c_Xh, 6.0_rp)
1232 call bdry_mask%mark_zone(this%msh%periodic)
1233 call bdry_mask%finalize()
1234 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1235 call bdry_mask%free()
1236
1237 do i = 1, this%bcs_prs%size()
1238 bci => this%bcs_prs%get(i)
1239 select type (bc => bci)
1240 type is (zero_dirichlet_t)
1241 call bdry_mask%init_from_components(this%c_Xh, 3.0_rp)
1242 call bdry_mask%mark_facets(bci%marked_facet)
1243 call bdry_mask%finalize()
1244 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1245 call bdry_mask%free()
1246 type is (dong_outflow_t)
1247 call bdry_mask%init_from_components(this%c_Xh, 3.0_rp)
1248 call bdry_mask%mark_facets(bci%marked_facet)
1249 call bdry_mask%finalize()
1250 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1251 call bdry_mask%free()
1252 type is (field_dirichlet_t)
1253 call bdry_mask%init_from_components(this%c_Xh, 8.0_rp)
1254 call bdry_mask%mark_facets(bci%marked_facet)
1255 call bdry_mask%finalize()
1256 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1257 call bdry_mask%free()
1258 end select
1259 end do
1260
1261 do i = 1, this%bcs_vel%size()
1262 bci => this%bcs_vel%get(i)
1263 select type (bc => bci)
1264 type is (zero_dirichlet_t)
1265 call bdry_mask%init_from_components(this%c_Xh, 1.0_rp)
1266 call bdry_mask%mark_facets(bci%marked_facet)
1267 call bdry_mask%finalize()
1268 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1269 call bdry_mask%free()
1270 type is (inflow_t)
1271 call bdry_mask%init_from_components(this%c_Xh, 2.0_rp)
1272 call bdry_mask%mark_facets(bci%marked_facet)
1273 call bdry_mask%finalize()
1274 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1275 call bdry_mask%free()
1276 type is (symmetry_t)
1277 call bdry_mask%init_from_components(this%c_Xh, 4.0_rp)
1278 call bdry_mask%mark_facets(bci%marked_facet)
1279 call bdry_mask%finalize()
1280 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1281 call bdry_mask%free()
1282 type is (usr_inflow_t)
1283 call bdry_mask%init_from_components(this%c_Xh, 5.0_rp)
1284 call bdry_mask%mark_facets(bci%marked_facet)
1285 call bdry_mask%finalize()
1286 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1287 call bdry_mask%free()
1288 type is (field_dirichlet_vector_t)
1289 call bdry_mask%init_from_components(this%c_Xh, 7.0_rp)
1290 call bdry_mask%mark_facets(bci%marked_facet)
1291 call bdry_mask%finalize()
1292 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1293 call bdry_mask%free()
1294 type is (shear_stress_t)
1295 call bdry_mask%init_from_components(this%c_Xh, 9.0_rp)
1296 call bdry_mask%mark_facets(bci%marked_facet)
1297 call bdry_mask%finalize()
1298 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1299 call bdry_mask%free()
1300 type is (wall_model_bc_t)
1301 call bdry_mask%init_from_components(this%c_Xh, 10.0_rp)
1302 call bdry_mask%mark_facets(bci%marked_facet)
1303 call bdry_mask%finalize()
1304 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1305 call bdry_mask%free()
1306 type is (blasius_t)
1307 call bdry_mask%init_from_components(this%c_Xh, 11.0_rp)
1308 call bdry_mask%mark_facets(bci%marked_facet)
1309 call bdry_mask%finalize()
1310 call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
1311 call bdry_mask%free()
1312 end select
1313 end do
1314
1315
1316 bdry_file = file_t('boundary_adjoint.fld')
1317 call bdry_file%write(bdry_field)
1318
1319 call this%scratch%relinquish_field(temp_index)
1320 end subroutine adjoint_fluid_pnpn_write_boundary_conditions
1321
1322 ! End of section to verify
1323 ! ========================================================================== !
1324
1325 subroutine rescale_fluid(fluid_data, scale)
1326
1327 class(adjoint_fluid_pnpn_t), intent(inout) :: fluid_data
1329 real(kind=rp), intent(in) :: scale
1330
1331 ! Local variables
1332 integer :: i
1333
1334 ! Scale the velocity fields
1335 if (neko_bcknd_device .eq. 1) then
1336 call device_cmult(fluid_data%u_adj%x_d, scale, fluid_data%u_adj%size())
1337 call device_cmult(fluid_data%v_adj%x_d, scale, fluid_data%v_adj%size())
1338 call device_cmult(fluid_data%w_adj%x_d, scale, fluid_data%w_adj%size())
1339 else
1340 call cmult(fluid_data%u_adj%x, scale, fluid_data%u_adj%size())
1341 call cmult(fluid_data%v_adj%x, scale, fluid_data%v_adj%size())
1342 call cmult(fluid_data%w_adj%x, scale, fluid_data%w_adj%size())
1343 end if
1344
1345 ! Scale the right hand sides
1346 if (neko_bcknd_device .eq. 1) then
1347 call device_cmult(fluid_data%f_adj_x%x_d, scale, &
1348 fluid_data%f_adj_x%size())
1349 call device_cmult(fluid_data%f_adj_y%x_d, scale, &
1350 fluid_data%f_adj_y%size())
1351 call device_cmult(fluid_data%f_adj_z%x_d, scale, &
1352 fluid_data%f_adj_z%size())
1353 ! HARRY
1354 ! maybe the abx's too
1355 call device_cmult(fluid_data%abx1%x_d, scale, fluid_data%abx1%size())
1356 call device_cmult(fluid_data%aby1%x_d, scale, fluid_data%aby1%size())
1357 call device_cmult(fluid_data%abz1%x_d, scale, fluid_data%abz1%size())
1358 call device_cmult(fluid_data%abx2%x_d, scale, fluid_data%abx2%size())
1359 call device_cmult(fluid_data%aby2%x_d, scale, fluid_data%aby2%size())
1360 call device_cmult(fluid_data%abz2%x_d, scale, fluid_data%abz2%size())
1361
1362 else
1363 call cmult(fluid_data%f_adj_x%x, scale, fluid_data%f_adj_x%size())
1364 call cmult(fluid_data%f_adj_y%x, scale, fluid_data%f_adj_y%size())
1365 call cmult(fluid_data%f_adj_z%x, scale, fluid_data%f_adj_z%size())
1366
1367 call cmult(fluid_data%abx1%x, scale, fluid_data%abx1%size())
1368 call cmult(fluid_data%aby1%x, scale, fluid_data%aby1%size())
1369 call cmult(fluid_data%abz1%x, scale, fluid_data%abz1%size())
1370
1371 call cmult(fluid_data%abx2%x, scale, fluid_data%abx2%size())
1372 call cmult(fluid_data%aby2%x, scale, fluid_data%aby2%size())
1373 call cmult(fluid_data%abz2%x, scale, fluid_data%abz2%size())
1374 end if
1375
1376 ! Scale the lag terms
1377 if (neko_bcknd_device .eq. 1) then
1378 do i = 1, fluid_data%ulag%size()
1379 call device_cmult(fluid_data%ulag%lf(i)%x_d, &
1380 scale, fluid_data%ulag%lf(i)%size())
1381 end do
1382
1383 do i = 1, fluid_data%vlag%size()
1384 call device_cmult(fluid_data%vlag%lf(i)%x_d, &
1385 scale, fluid_data%vlag%lf(i)%size())
1386 end do
1387
1388 do i = 1, fluid_data%wlag%size()
1389 call device_cmult(fluid_data%wlag%lf(i)%x_d, &
1390 scale, fluid_data%wlag%lf(i)%size())
1391 end do
1392 else
1393 do i = 1, fluid_data%ulag%size()
1394 call cmult(fluid_data%ulag%lf(i)%x, &
1395 scale, fluid_data%ulag%lf(i)%size())
1396 end do
1397
1398 do i = 1, fluid_data%vlag%size()
1399 call cmult(fluid_data%vlag%lf(i)%x, &
1400 scale, fluid_data%vlag%lf(i)%size())
1401 end do
1402
1403 do i = 1, fluid_data%wlag%size()
1404 call cmult(fluid_data%wlag%lf(i)%x, &
1405 scale, fluid_data%wlag%lf(i)%size())
1406 end do
1407 end if
1408
1409 end subroutine rescale_fluid
1410
1411 function norm(x, y, z, B, volume, n)
1412 use mpi_f08, only: mpi_in_place
1413
1414 integer, intent(in) :: n
1415 real(kind=rp), dimension(n), intent(in) :: x, y, z
1416 real(kind=rp), dimension(n), intent(in) :: b
1417 real(kind=rp), intent(in) :: volume
1418
1419 real(kind=rp) :: norm
1420
1421 norm = vlsc3(x, x, b, n) + vlsc3(y, y, b, n) + vlsc3(z, z, b, n)
1422
1423 call mpi_allreduce(mpi_in_place, norm, 1, &
1424 mpi_real_precision, mpi_sum, mpi_comm_world)
1425
1426 norm = sqrt(norm / volume)
1427 end function norm
1428
1429 function device_norm(x_d, y_d, z_d, B_d, volume, n)
1430 use mpi_f08, only: mpi_in_place
1431
1432 type(c_ptr), intent(in) :: x_d, y_d, z_d
1433 type(c_ptr), intent(in) :: b_d
1434 real(kind=rp), intent(in) :: volume
1435 integer, intent(in) :: n
1436
1437 real(kind=rp) :: device_norm
1438
1439 device_norm = device_vlsc3(x_d, x_d, b_d, n) + &
1440 device_vlsc3(y_d, y_d, b_d, n) + &
1441 device_vlsc3(z_d, z_d, b_d, n)
1442
1443 call mpi_allreduce(mpi_in_place, device_norm, 1, &
1444 mpi_real_precision, mpi_sum, mpi_comm_world)
1445
1446 device_norm = sqrt(device_norm / volume)
1447
1448 end function device_norm
1449
1453 subroutine power_iterations_compute(this, t, tstep)
1454 class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
1455 real(kind=rp), intent(in) :: t
1456 integer, intent(in) :: tstep
1457
1458 ! Local variables
1459 real(kind=rp) :: scaling_factor
1460 real(kind=rp) :: norm_l2, norm_l2_base
1461 character(len=256) :: log_message
1462 type(vector_t) :: data_line
1463 integer :: n
1464
1465 n = this%c_Xh%dof%size()
1466 if (tstep .eq. 1) then
1467 if (neko_bcknd_device .eq. 1) then
1468 norm_l2_base = device_norm(this%u_adj%x_d, this%v_adj%x_d, &
1469 this%w_adj%x_d, &
1470 this%c_Xh%B_d, this%c_Xh%volume, n)
1471 else
1472 norm_l2_base = this%norm_scaling * norm(this%u_adj%x, this%v_adj%x, &
1473 this%w_adj%x, &
1474 this%c_Xh%B, this%c_Xh%volume, n)
1475 end if
1476 if (this%norm_target .lt. 0.0_rp) then
1477 this%norm_target = norm_l2_base
1478 end if
1479
1480 this%norm_l2_upper = this%norm_tolerance * this%norm_target
1481 this%norm_l2_lower = this%norm_target / this%norm_tolerance
1482
1483 end if
1484
1485 ! Compute the norm of the velocity field and eigenvalue estimate
1486 if (neko_bcknd_device .eq. 1) then
1487 norm_l2 = device_norm(this%u_adj%x_d, this%v_adj%x_d, this%w_adj%x_d, &
1488 this%c_Xh%B_d, this%c_Xh%volume, n)
1489 else
1490 norm_l2 = norm(this%u_adj%x, this%v_adj%x, this%w_adj%x, &
1491 this%c_Xh%B, this%c_Xh%volume, n)
1492 end if
1493 norm_l2 = sqrt(this%norm_scaling) * norm_l2
1494 scaling_factor = 1.0_rp
1495
1496 ! Rescale the flow if necessary
1497 if (norm_l2 .gt. this%norm_l2_upper &
1498 .or. norm_l2 .lt. this%norm_l2_lower) then
1499 scaling_factor = this%norm_target / norm_l2
1500 call rescale_fluid(this, scaling_factor)
1501 norm_l2 = this%norm_target
1502
1503 if (tstep .eq. 1) then
1504 scaling_factor = 1.0_rp
1505 end if
1506 end if
1507
1508 ! Log the results
1509 !call neko_log%section('Power Iterations', lvl = NEKO_LOG_DEBUG)
1510 call neko_log%section('Power Iterations')
1511
1512 write (log_message, '(A7,E20.14)') 'Norm: ', norm_l2
1513 call neko_log%message(log_message, lvl = neko_log_debug)
1514 write (log_message, '(A7,E20.14)') 'Scaling: ', scaling_factor
1515 call neko_log%message(log_message, lvl = neko_log_debug)
1516
1517 ! Save to file
1518 call data_line%init(2)
1519 data_line%x = [norm_l2, scaling_factor]
1520 call this%file_output%write(data_line, t)
1521
1522 !call neko_log%end_section('Power Iterations', lvl = NEKO_LOG_DEBUG)
1523 call neko_log%end_section('Power Iterations')
1524 end subroutine power_iterations_compute
1525
1526
1527end module adjoint_fluid_pnpn
adjoint_fluid_pnpn
Adjoint Pn/Pn formulation.
Definition adjoint_fluid_pnpn.f90:34
adjoint_fluid_pnpn::adjoint_fluid_pnpn_free
subroutine adjoint_fluid_pnpn_free(this)
Definition adjoint_fluid_pnpn.f90:636
adjoint_fluid_pnpn::adjoint_fluid_pnpn_init
subroutine adjoint_fluid_pnpn_init(this, msh, lx, params, user, time_scheme)
Boundary condition factory for pressure.
Definition adjoint_fluid_pnpn.f90:284
adjoint_fluid_pnpn::adjoint_fluid_pnpn_write_boundary_conditions
subroutine adjoint_fluid_pnpn_write_boundary_conditions(this)
Write a field with boundary condition specifications.
Definition adjoint_fluid_pnpn.f90:1182
adjoint_fluid_pnpn::power_iterations_compute
subroutine power_iterations_compute(this, t, tstep)
Compute the power_iterations field.
Definition adjoint_fluid_pnpn.f90:1454
adjoint_fluid_pnpn::adjoint_fluid_pnpn_step
subroutine adjoint_fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
Advance fluid simulation in time.
Definition adjoint_fluid_pnpn.f90:720
adjoint_fluid_pnpn::adjoint_fluid_pnpn_setup_bcs
subroutine adjoint_fluid_pnpn_setup_bcs(this, user, params)
Sets up the boundary condition for the scheme.
Definition adjoint_fluid_pnpn.f90:967
adjoint_fluid_pnpn::adjoint_fluid_pnpn_restart
subroutine adjoint_fluid_pnpn_restart(this, dtlag, tlag)
Definition adjoint_fluid_pnpn.f90:483
adjoint_fluid_scheme
Fluid formulations.
Definition adjoint_fluid_scheme.f90:34
advection_adjoint_fctry
Contains the factory routine for advection_t children.
Definition advection_adjoint_fctry.f90:34
advection_adjoint_fctry::advection_adjoint_factory
subroutine, public advection_adjoint_factory(object, json, coef)
A factory for advection_t decendants.
Definition advection_adjoint_fctry.f90:58
advection_adjoint
Subroutines to add advection terms to the RHS of a transport equation.
Definition advection_adjoint.f90:34
json_utils_ext
Definition json_utils_ext.f90:1
json_utils_ext::json_key_fallback
character(len=:) function, allocatable, public json_key_fallback(json, lookup, fallback)
Create a json_string based on fallback logic. If the lookup key is present in the json object,...
Definition json_utils_ext.f90:20
user
User defined user region.
Definition user.f90:2
adjoint_fluid_pnpn::adjoint_fluid_pnpn_t
Definition adjoint_fluid_pnpn.f90:96
adjoint_fluid_scheme::adjoint_fluid_scheme_t
Base type of all fluid formulations.
Definition adjoint_fluid_scheme.f90:90
advection_adjoint::advection_adjoint_t
Base abstract type for computing the advection operator.
Definition advection_adjoint.f90:43