adjoint_fluid_pnpn.f90

! Copyright (c) 2024, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module adjoint_fluid_pnpn
  use, intrinsic :: iso_fortran_env, only: error_unit
  use coefs, only: coef_t
  use symmetry, only: symmetry_t
  use field_registry, only: neko_field_registry
  use logger, only: neko_log, log_size, neko_log_debug
  use num_types, only: rp
  use krylov, only: ksp_monitor_t
  use pnpn_residual, only: pnpn_prs_res_t, pnpn_vel_res_t, &
       pnpn_prs_res_factory, pnpn_vel_res_factory, &
       pnpn_prs_res_stress_factory, pnpn_vel_res_stress_factory
  use rhs_maker, only: rhs_maker_sumab_t, rhs_maker_bdf_t, rhs_maker_ext_t, &
       rhs_maker_oifs_t, rhs_maker_sumab_fctry, rhs_maker_bdf_fctry, &
       rhs_maker_ext_fctry, rhs_maker_oifs_fctry
  ! use adjoint_volflow, only: adjoint_volflow_t
  use adjoint_fluid_scheme, only: adjoint_fluid_scheme_t
  use device_mathops, only: device_opcolv, device_opadd2cm
  use fluid_aux, only: fluid_step_info
  use time_scheme_controller, only: time_scheme_controller_t
  use projection, only: projection_t
  use device, only : device_memcpy, host_to_device, device_event_sync,&
       device_event_create, device_event_destroy, device_stream_wait_event, &
       glb_cmd_queue
  use advection_adjoint, only: advection_adjoint_t
  use advection_adjoint_fctry, only: advection_adjoint_factory
  use profiler, only: profiler_start_region, profiler_end_region
  use json_module, only: json_file, json_core, json_value
  use json_utils, only: json_get, json_get_or_default, json_extract_item, &
       json_extract_object
  use json_module, only: json_file
  use ax_product, only: ax_t, ax_helm_factory
  use field, only: field_t
  use dirichlet, only: dirichlet_t
  use shear_stress, only: shear_stress_t
  use wall_model_bc, only: wall_model_bc_t
  use facet_normal, only: facet_normal_t
  use non_normal, only: non_normal_t
  use mesh, only: mesh_t
  use user_intf, only: user_t
  use time_step_controller, only: time_step_controller_t
  use gs_ops, only: gs_op_add
  use neko_config, only: neko_bcknd_device
  use mathops, only: opadd2cm, opcolv
  use bc_list, only: bc_list_t
  use zero_dirichlet, only: zero_dirichlet_t
  use utils, only: neko_error, neko_type_error
  use field_math, only: field_add2, field_copy
  use bc, only: bc_t
  use file, only: file_t
  use operators, only: ortho
  use vector, only: vector_t
  use device_math, only: device_vlsc3, device_cmult
  use math, only: vlsc3, cmult
  use json_utils_ext, only: json_key_fallback
  use, intrinsic :: iso_c_binding, only: c_ptr, c_null_ptr, c_associated
  use comm, only: neko_comm, mpi_real_precision
  use mpi_f08, only: mpi_sum, mpi_max, mpi_allreduce, mpi_comm_world, &
       mpi_integer, mpi_logical, mpi_lor
 
  implicit none
  private
 
 
  type, public, extends(adjoint_fluid_scheme_t) :: adjoint_fluid_pnpn_t
 
     type(field_t) :: p_res, u_res, v_res, w_res
 
     type(field_t) :: dp, du, dv, dw
 
     type(field_t), pointer :: u_b => null() 
     type(field_t), pointer :: v_b => null() 
     type(field_t), pointer :: w_b => null() 
     type(field_t), pointer :: p_b => null() 
 
     !
     ! Implicit operators, i.e. the left-hand-side of the Helmholz problem.
     !
 
     ! Coupled Helmholz operator for velocity
     class(ax_t), allocatable :: ax_vel
     ! Helmholz operator for pressure
     class(ax_t), allocatable :: ax_prs
 
     !
     ! Projections for solver speed-up
     !
 
     type(projection_t) :: proj_prs
     type(projection_t) :: proj_u
     type(projection_t) :: proj_v
     type(projection_t) :: proj_w
 
     !
     ! Special Karniadakis scheme boundary conditions in the pressure equation
     !
 
     type(facet_normal_t) :: bc_prs_surface
 
     type(facet_normal_t) :: bc_sym_surface
 
     !
     ! Boundary conditions and  lists for residuals and solution increments
     !
 
     type(zero_dirichlet_t) :: bc_vel_res
     type(zero_dirichlet_t) :: bc_du
     type(zero_dirichlet_t) :: bc_dv
     type(zero_dirichlet_t) :: bc_dw
     type(zero_dirichlet_t) :: bc_dp
 
     type(bc_list_t) :: bclst_vel_res
     type(bc_list_t) :: bclst_du
     type(bc_list_t) :: bclst_dv
     type(bc_list_t) :: bclst_dw
     type(bc_list_t) :: bclst_dp
 
 
     ! Checker for wether we have a strong pressure bc. If not, the pressure
     ! is demeaned at every time step.
     logical :: prs_dirichlet = .false.
 
 
     ! The advection operator.
     class(advection_adjoint_t), allocatable :: adv
 
     ! Time OIFS interpolation scheme for advection.
     logical :: oifs
 
     ! Time variables
     type(field_t) :: abx1, aby1, abz1
     type(field_t) :: abx2, aby2, abz2
 
     ! Advection terms for the oifs method
     type(field_t) :: advx, advy, advz
 
     class(pnpn_prs_res_t), allocatable :: prs_res
 
     class(pnpn_vel_res_t), allocatable :: vel_res
 
     class(rhs_maker_sumab_t), allocatable :: sumab
 
     class(rhs_maker_ext_t), allocatable :: makeabf
 
     class(rhs_maker_bdf_t), allocatable :: makebdf
 
     class(rhs_maker_oifs_t), allocatable :: makeoifs
 
     !  type(fluid_volflow_t) :: vol_flow
     type(c_ptr) :: event = c_null_ptr
 
     ! ======================================================================= !
     ! Addressable attributes
 
     real(kind=rp) :: norm_scaling 
     real(kind=rp) :: norm_target 
     real(kind=rp) :: norm_tolerance 
 
     ! ======================================================================= !
     ! Definition of shorthands and local variables
 
     real(kind=rp) :: norm_l2_base
 
     real(kind=rp) :: norm_l2_upper
     real(kind=rp) :: norm_l2_lower
 
     type(file_t) :: file_output
 
   contains
     procedure, pass(this) :: init => adjoint_fluid_pnpn_init
     procedure, pass(this) :: free => adjoint_fluid_pnpn_free
     procedure, pass(this) :: step => adjoint_fluid_pnpn_step
     procedure, pass(this) :: restart => adjoint_fluid_pnpn_restart
     procedure, pass(this) :: setup_bcs => adjoint_fluid_pnpn_setup_bcs
     procedure, pass(this) :: write_boundary_conditions => &
          adjoint_fluid_pnpn_write_boundary_conditions
 
     procedure, public, pass(this) :: pw_compute_ => power_iterations_compute
 
  end type adjoint_fluid_pnpn_t
 
  interface pressure_bc_factory
     module subroutine pressure_bc_factory(object, scheme, json, coef, user)
       class(bc_t), pointer, intent(inout) :: object
       type(adjoint_fluid_pnpn_t), intent(in) :: scheme
       type(json_file), intent(inout) :: json
       type(coef_t), intent(in) :: coef
       type(user_t), intent(in) :: user
     end subroutine pressure_bc_factory
  end interface pressure_bc_factory
 
  interface velocity_bc_factory
     module subroutine velocity_bc_factory(object, scheme, json, coef, user)
       class(bc_t), pointer, intent(inout) :: object
       type(adjoint_fluid_pnpn_t), intent(in) :: scheme
       type(json_file), intent(inout) :: json
       type(coef_t), intent(in) :: coef
       type(user_t), intent(in) :: user
     end subroutine velocity_bc_factory
  end interface velocity_bc_factory
 
contains
 
  subroutine adjoint_fluid_pnpn_init(this, msh, lx, params, user, time_scheme)
    class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
    type(mesh_t), target, intent(inout) :: msh
    integer, intent(in) :: lx
    type(json_file), target, intent(inout) :: params
    type(user_t), target, intent(in) :: user
    type(time_scheme_controller_t), target, intent(in) :: time_scheme
    character(len=15), parameter :: scheme = 'Adjoint (Pn/Pn)'
    real(kind=rp) :: abs_tol
    character(len=LOG_SIZE) :: log_buf
    integer :: solver_maxiter
    character(len=:), allocatable :: solver_type, precon_type
    logical :: monitor, found
    logical :: advection
    type(json_file) :: precon_params
 
    ! Temporary field pointers
    character(len=:), allocatable :: file_name
    character(len=256) :: header_line
 
    call this%free()
 
    ! Initialize base class
    call this%init_base(msh, lx, params, scheme, user, .true.)
 
    ! Add pressure field to the registry. For this scheme it is in the same
    ! Xh as the velocity
    call neko_field_registry%add_field(this%dm_Xh, 'p_adj')
    this%p_adj => neko_field_registry%get_field('p_adj')
 
    !
    ! Select governing equations via associated residual and Ax types
    !
 
    if (this%variable_material_properties .eqv. .true.) then
       ! Setup backend dependent Ax routines
       call ax_helm_factory(this%Ax_vel, full_formulation = .true.)
 
       ! Setup backend dependent prs residual routines
       call pnpn_prs_res_stress_factory(this%prs_res)
 
       ! Setup backend dependent vel residual routines
       call pnpn_vel_res_stress_factory(this%vel_res)
    else
       ! Setup backend dependent Ax routines
       call ax_helm_factory(this%Ax_vel, full_formulation = .false.)
 
       ! Setup backend dependent prs residual routines
       call pnpn_prs_res_factory(this%prs_res)
 
       ! Setup backend dependent vel residual routines
       call pnpn_vel_res_factory(this%vel_res)
    end if
 
    ! Setup Ax for the pressure
    call ax_helm_factory(this%Ax_prs, full_formulation = .false.)
 
 
    ! Setup backend dependent summation of AB/BDF
    call rhs_maker_sumab_fctry(this%sumab)
 
    ! Setup backend dependent summation of extrapolation scheme
    call rhs_maker_ext_fctry(this%makeabf)
 
    ! Setup backend depenent contributions to F from lagged BD terms
    call rhs_maker_bdf_fctry(this%makebdf)
 
    ! Setup backend dependent summations of the OIFS method
    call rhs_maker_oifs_fctry(this%makeoifs)
 
    ! Initialize variables specific to this plan
    associate(xh_lx => this%Xh%lx, xh_ly => this%Xh%ly, xh_lz => this%Xh%lz, &
         dm_xh => this%dm_Xh, nelv => this%msh%nelv)
 
      call this%p_res%init(dm_xh, "p_res")
      call this%u_res%init(dm_xh, "u_res")
      call this%v_res%init(dm_xh, "v_res")
      call this%w_res%init(dm_xh, "w_res")
      call this%abx1%init(dm_xh, "abx1")
      call this%aby1%init(dm_xh, "aby1")
      call this%abz1%init(dm_xh, "abz1")
      call this%abx2%init(dm_xh, "abx2")
      call this%aby2%init(dm_xh, "aby2")
      call this%abz2%init(dm_xh, "abz2")
      call this%advx%init(dm_xh, "advx")
      call this%advy%init(dm_xh, "advy")
      call this%advz%init(dm_xh, "advz")
    end associate
 
    call this%du%init(this%dm_Xh, 'du')
    call this%dv%init(this%dm_Xh, 'dv')
    call this%dw%init(this%dm_Xh, 'dw')
    call this%dp%init(this%dm_Xh, 'dp')
 
    ! Set up boundary conditions
    call this%setup_bcs(user, params)
 
    ! Check if we need to output boundaries
    call json_get_or_default(params, 'case.output_boundary', found, .false.)
    if (found) call this%write_boundary_conditions()
 
    ! Intialize projection space
    if (this%variable_material_properties .and. &
         this%vel_projection_dim .gt. 0) then
       call neko_error("Velocity projection not available for full stress &
       &formulation")
    end if
 
 
    call this%proj_prs%init(this%dm_Xh%size(), this%pr_projection_dim, &
         this%pr_projection_activ_step)
 
    call this%proj_u%init(this%dm_Xh%size(), this%vel_projection_dim, &
         this%vel_projection_activ_step)
    call this%proj_v%init(this%dm_Xh%size(), this%vel_projection_dim, &
         this%vel_projection_activ_step)
    call this%proj_w%init(this%dm_Xh%size(), this%vel_projection_dim, &
         this%vel_projection_activ_step)
 
 
    ! Add lagged term to checkpoint
    call this%chkp%add_lag(this%ulag, this%vlag, this%wlag)
 
    ! Determine the time-interpolation scheme
    call json_get_or_default(params, 'case.numerics.oifs', this%oifs, .false.)
 
    ! Initialize the advection factory
    call json_get_or_default(params, 'case.fluid.advection', advection, .true.)
 
    ! Todo: make sure this actually follow the forward advection factory
    call advection_adjoint_factory(this%adv, params, this%c_Xh)
 
    if (params%valid_path('case.fluid.flow_rate_force')) then
       call neko_error("Flow rate forcing not available for adjoint_fluid_pnpn")
       !  call this%vol_flow%init(this%dm_Xh, params)
    end if
 
    ! Setup pressure solver
    call neko_log%section("Pressure solver")
 
    call json_get_or_default(params, &
         'case.fluid.pressure_solver.max_iterations', &
         solver_maxiter, 800)
    call json_get(params, 'case.fluid.pressure_solver.type', solver_type)
    call json_get(params, 'case.fluid.pressure_solver.preconditioner.type', &
         precon_type)
    call json_extract_object(params, &
         'case.fluid.pressure_solver.preconditioner', precon_params)
    call json_get(params, 'case.fluid.pressure_solver.absolute_tolerance', &
         abs_tol)
    call json_get_or_default(params, 'case.fluid.velocity_solver.monitor', &
         monitor, .false.)
    call neko_log%message('Type       : ('// trim(solver_type) // &
         ', ' // trim(precon_type) // ')')
    write(log_buf, '(A,ES13.6)') 'Abs tol    :', abs_tol
    call neko_log%message(log_buf)
 
    call this%solver_factory(this%ksp_prs, this%dm_Xh%size(), &
         solver_type, solver_maxiter, abs_tol, monitor)
    call this%precon_factory_(this%pc_prs, this%ksp_prs, &
         this%c_Xh, this%dm_Xh, this%gs_Xh, this%bcs_prs, &
         precon_type, precon_params)
 
    call neko_log%end_section()
 
    ! Setup gather-scatter event (only relevant on devices)
    if (neko_bcknd_device .eq. 1) then
       call device_event_create(this%event, 2)
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Handling the rescaling and baseflow
 
    ! Read the norm scaling from the json file
    call json_get_or_default(params, 'norm_scaling', &
         this%norm_scaling, 0.5_rp)
 
    ! The baseflow is the solution to the forward.
    ! Userdefined baseflows can be invoked via setting initial conditions
    ! call neko_field_registry%add_field(this%dm_Xh, 'u')
    ! call neko_field_registry%add_field(this%dm_Xh, 'v')
    ! call neko_field_registry%add_field(this%dm_Xh, 'w')
    ! call neko_field_registry%add_field(this%dm_Xh, 'p')
    this%u_b => neko_field_registry%get_field('u')
    this%v_b => neko_field_registry%get_field('v')
    this%w_b => neko_field_registry%get_field('w')
    this%p_b => neko_field_registry%get_field('p')
 
    ! Read the json file
    call json_get_or_default(params, 'norm_target', &
         this%norm_target, -1.0_rp)
    call json_get_or_default(params, 'norm_tolerance', &
         this%norm_tolerance, 10.0_rp)
 
    ! Build the header
    call json_get_or_default(params, 'output_file', &
         file_name, 'power_iterations.csv')
    this%file_output = file_t(trim(file_name))
    write(header_line, '(A)') 'Time, Norm, Scaling'
    call this%file_output%set_header(header_line)
 
  end subroutine adjoint_fluid_pnpn_init
 
  subroutine adjoint_fluid_pnpn_restart(this, dtlag, tlag)
    class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
    real(kind=rp) :: dtlag(10), tlag(10)
    integer :: i, j, n
 
    n = this%u_adj%dof%size()
    if (allocated(this%chkp%previous_mesh%elements) .or. &
         this%chkp%previous_Xh%lx .ne. this%Xh%lx) then
       associate(u => this%u_adj, v => this%v_adj, w => this%w_adj, &
            p => this%p_adj, c_xh => this%c_Xh, ulag => this%ulag, &
            vlag => this%vlag, wlag => this%wlag)
         do concurrent(j = 1:n)
            u%x(j,1,1,1) = u%x(j,1,1,1) * c_xh%mult(j,1,1,1)
            v%x(j,1,1,1) = v%x(j,1,1,1) * c_xh%mult(j,1,1,1)
            w%x(j,1,1,1) = w%x(j,1,1,1) * c_xh%mult(j,1,1,1)
            p%x(j,1,1,1) = p%x(j,1,1,1) * c_xh%mult(j,1,1,1)
         end do
         do i = 1, this%ulag%size()
            do concurrent(j = 1:n)
               ulag%lf(i)%x(j,1,1,1) = ulag%lf(i)%x(j,1,1,1) &
                    * c_xh%mult(j,1,1,1)
               vlag%lf(i)%x(j,1,1,1) = vlag%lf(i)%x(j,1,1,1) &
                    * c_xh%mult(j,1,1,1)
               wlag%lf(i)%x(j,1,1,1) = wlag%lf(i)%x(j,1,1,1) &
                    * c_xh%mult(j,1,1,1)
            end do
         end do
       end associate
    end if
 
    if (neko_bcknd_device .eq. 1) then
       associate(u => this%u_adj, v => this%v_adj, w => this%w_adj, &
            ulag => this%ulag, vlag => this%vlag, wlag => this%wlag,&
            p => this%p_adj)
         call device_memcpy(u%x, u%x_d, u%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(v%x, v%x_d, v%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(w%x, w%x_d, w%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(p%x, p%x_d, p%dof%size(), &
              host_to_device, sync = .false.)
         call device_memcpy(ulag%lf(1)%x, ulag%lf(1)%x_d, &
              u%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(ulag%lf(2)%x, ulag%lf(2)%x_d, &
              u%dof%size(), host_to_device, sync = .false.)
 
         call device_memcpy(vlag%lf(1)%x, vlag%lf(1)%x_d, &
              v%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(vlag%lf(2)%x, vlag%lf(2)%x_d, &
              v%dof%size(), host_to_device, sync = .false.)
 
         call device_memcpy(wlag%lf(1)%x, wlag%lf(1)%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(wlag%lf(2)%x, wlag%lf(2)%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%abx1%x, this%abx1%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%abx2%x, this%abx2%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%aby1%x, this%aby1%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%aby2%x, this%aby2%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%abz1%x, this%abz1%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%abz2%x, this%abz2%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%advx%x, this%advx%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%advy%x, this%advy%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
         call device_memcpy(this%advz%x, this%advz%x_d, &
              w%dof%size(), host_to_device, sync = .false.)
       end associate
    end if
 
    ! Make sure that continuity is maintained (important for interpolation)
    ! Do not do this for lagged rhs
    ! (derivatives are not necessairly coninous across elements)
 
    if (allocated(this%chkp%previous_mesh%elements) &
         .or. this%chkp%previous_Xh%lx .ne. this%Xh%lx) then
       call this%gs_Xh%op(this%u_adj, gs_op_add)
       call this%gs_Xh%op(this%v_adj, gs_op_add)
       call this%gs_Xh%op(this%w_adj, gs_op_add)
       call this%gs_Xh%op(this%p_adj, gs_op_add)
 
       do i = 1, this%ulag%size()
          call this%gs_Xh%op(this%ulag%lf(i), gs_op_add)
          call this%gs_Xh%op(this%vlag%lf(i), gs_op_add)
          call this%gs_Xh%op(this%wlag%lf(i), gs_op_add)
       end do
    end if
 
    !! If we would decide to only restart from lagged fields instead of asving
    !! abx1, aby1 etc.
    !! Observe that one also needs to recompute the focing at the old time steps
    !u_temp = this%ulag%lf(2)
    !v_temp = this%vlag%lf(2)
    !w_temp = this%wlag%lf(2)
    !! Compute the source terms
    !call this%source_term%compute(tlag(2), -1)
    !
    !! Pre-multiply the source terms with the mass matrix.
    !if (NEKO_BCKND_DEVICE .eq. 1) then
    !   call device_opcolv(this%f_adj_x%x_d, this%f_adj_y%x_d, this%f_adj_z%x_d,
    ! this%c_Xh%B_d, this%msh%gdim, n)
    !else
    !   call opcolv(this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, this%c_Xh%B,
    ! this%msh%gdim, n)
    !end if
 
    !! Add the advection operators to the right-hand-side.
    !call this%adv%compute(u_temp, v_temp, w_temp, &
    !                      this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, &
    !                      this%Xh, this%c_Xh, this%dm_Xh%size())
    !this%abx2 = this%f_x
    !this%aby2 = this%f_y
    !this%abz2 = this%f_z
    !
    !u_temp = this%ulag%lf(1)
    !v_temp = this%vlag%lf(1)
    !w_temp = this%wlag%lf(1)
    !call this%source_term%compute(tlag(1), 0)
 
    !! Pre-multiply the source terms with the mass matrix.
    !if (NEKO_BCKND_DEVICE .eq. 1) then
    !   call device_opcolv(this%f_adj_x%x_d, this%f_adj_y%x_d, this%f_adj_z%x_d,
    ! this%c_Xh%B_d, this%msh%gdim, n)
    !else
    !   call opcolv(this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, this%c_Xh%B,
    ! this%msh%gdim, n)
    !end if
 
    !! Pre-multiply the source terms with the mass matrix.
    !if (NEKO_BCKND_DEVICE .eq. 1) then
    !   call device_opcolv(this%f_adj_x%x_d, this%f_adj_y%x_d, this%f_adj_z%x_d,
    ! this%c_Xh%B_d, this%msh%gdim, n)
    !else
    !   call opcolv(this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, this%c_Xh%B,
    ! this%msh%gdim, n)
    !end if
 
    !call this%adv%compute(u_temp, v_temp, w_temp, &
    !                      this%f_adj_x%x, this%f_adj_y%x, this%f_adj_z%x, &
    !                      this%Xh, this%c_Xh, this%dm_Xh%size())
    !this%abx1 = this%f_x
    !this%aby1 = this%f_y
    !this%abz1 = this%f_z
 
  end subroutine adjoint_fluid_pnpn_restart
 
  subroutine adjoint_fluid_pnpn_free(this)
    class(adjoint_fluid_pnpn_t), intent(inout) :: this
 
    !Deallocate velocity and pressure fields
    call this%scheme_free()
 
    call this%bc_prs_surface%free()
    call this%bc_sym_surface%free()
    call this%bclst_vel_res%free()
    call this%bclst_dp%free()
    call this%proj_prs%free()
    call this%proj_u%free()
    call this%proj_v%free()
    call this%proj_w%free()
 
    call this%p_res%free()
    call this%u_res%free()
    call this%v_res%free()
    call this%w_res%free()
 
    call this%du%free()
    call this%dv%free()
    call this%dw%free()
    call this%dp%free()
 
    call this%abx1%free()
    call this%aby1%free()
    call this%abz1%free()
 
    call this%abx2%free()
    call this%aby2%free()
    call this%abz2%free()
 
    call this%advx%free()
    call this%advy%free()
    call this%advz%free()
 
    if (allocated(this%Ax_vel)) then
       deallocate(this%Ax_vel)
    end if
 
    if (allocated(this%Ax_prs)) then
       deallocate(this%Ax_prs)
    end if
 
    if (allocated(this%prs_res)) then
       deallocate(this%prs_res)
    end if
 
    if (allocated(this%vel_res)) then
       deallocate(this%vel_res)
    end if
 
    if (allocated(this%sumab)) then
       deallocate(this%sumab)
    end if
 
    if (allocated(this%makeabf)) then
       deallocate(this%makeabf)
    end if
 
    if (allocated(this%makebdf)) then
       deallocate(this%makebdf)
    end if
 
    if (allocated(this%makeoifs)) then
       deallocate(this%makeoifs)
    end if
 
    ! call this%vol_flow%free()
    if (neko_bcknd_device .eq. 1) then
       if (c_associated(this%event)) then
          call device_event_destroy(this%event)
       end if
    end if
 
  end subroutine adjoint_fluid_pnpn_free
 
  subroutine adjoint_fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
    class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
    real(kind=rp), intent(in) :: t
    integer, intent(in) :: tstep
    real(kind=rp), intent(in) :: dt
    type(time_scheme_controller_t), intent(in) :: ext_bdf
    type(time_step_controller_t), intent(in) :: dt_controller
    ! number of degrees of freedom
    integer :: n
    ! Solver results monitors (pressure + 3 velocity)
    type(ksp_monitor_t) :: ksp_results(4)
    ! Extrapolated velocity for the pressure residual
    type(field_t), pointer :: u_e, v_e, w_e
    ! Indices for tracking temporary fields
    integer :: temp_indices(3)
 
    if (this%freeze) return
 
    n = this%dm_Xh%size()
 
    call profiler_start_region('Adjoint', 1)
    associate(u => this%u_adj, v => this%v_adj, w => this%w_adj, &
         p => this%p_adj, &
         du => this%du, dv => this%dv, dw => this%dw, dp => this%dp, &
         u_b => this%u_b, v_b => this%v_b, w_b => this%w_b, &
         u_res => this%u_res, v_res => this%v_res, w_res => this%w_res, &
         p_res => this%p_res, ax_vel => this%Ax_vel, ax_prs => this%Ax_prs, &
         xh => this%Xh, &
         c_xh => this%c_Xh, dm_xh => this%dm_Xh, gs_xh => this%gs_Xh, &
         ulag => this%ulag, vlag => this%vlag, wlag => this%wlag, &
         msh => this%msh, prs_res => this%prs_res, &
         source_term => this%source_term, &
         vel_res => this%vel_res, sumab => this%sumab, &
         makeabf => this%makeabf, makebdf => this%makebdf, &
         vel_projection_dim => this%vel_projection_dim, &
         pr_projection_dim => this%pr_projection_dim, &
         rho => this%rho, mu => this%mu, oifs => this%oifs, &
         rho_field => this%rho, mu_field => this%mu, &
         f_x => this%f_adj_x, f_y => this%f_adj_y, f_z => this%f_adj_z, &
         if_variable_dt => dt_controller%if_variable_dt, &
         dt_last_change => dt_controller%dt_last_change, &
         event => this%event)
 
      ! Get temporary arrays
      call this%scratch%request_field(u_e, temp_indices(1))
      call this%scratch%request_field(v_e, temp_indices(2))
      call this%scratch%request_field(w_e, temp_indices(3))
      call sumab%compute_fluid(u_e, v_e, w_e, u, v, w, &
           ulag, vlag, wlag, ext_bdf%advection_coeffs, ext_bdf%nadv)
 
      ! Compute the source terms
      call this%source_term%compute(t, tstep)
 
      ! Add Neumann bc contributions to the RHS
      call this%bcs_vel%apply_vector(f_x%x, f_y%x, f_z%x, &
           this%dm_Xh%size(), t, tstep, strong = .false.)
 
      if (oifs) then
         call neko_error("OIFS not implemented for adjoint")
         !  ! Add the advection operators to the right-hand-side.
         !  call this%adv%compute_adjoint(u, v, w, &
         !       this%advx, this%advy, this%advz, &
         !       Xh, this%c_Xh, dm_Xh%size(), dt)
 
         ! At this point the RHS contains the sum of the advection operator and
         ! additional source terms, evaluated using the velocity field from the
         ! previous time-step. Now, this value is used in the explicit time
         ! scheme to advance both terms in time.
         !  call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
         !       this%abx2, this%aby2, this%abz2, &
         !       f_x%x, f_y%x, f_z%x, &
         !       rho, ext_bdf%advection_coeffs, n)
 
         ! Now, the source terms from the previous time step are added to the
         ! RHS.
         !  call makeoifs%compute_fluid(this%advx%x, this%advy%x, this%advz%x, &
         !       f_x%x, f_y%x, f_z%x, &
         !       rho, dt, n)
      else
         call this%adv%compute_adjoint(u, v, w, u_b, v_b, w_b, &
              f_x, f_y, f_z, &
              xh, this%c_Xh, dm_xh%size())
 
         ! At this point the RHS contains the sum of the advection operator and
         ! additional source terms, evaluated using the velocity field from the
         ! previous time-step. Now, this value is used in the explicit time
         ! scheme to advance both terms in time.
         call makeabf%compute_fluid(this%abx1, this%aby1, this%abz1,&
              this%abx2, this%aby2, this%abz2, &
              f_x%x, f_y%x, f_z%x, &
              rho%x(1,1,1,1), ext_bdf%advection_coeffs, n)
 
         ! Add the RHS contributions coming from the BDF scheme.
         call makebdf%compute_fluid(ulag, vlag, wlag, f_x%x, f_y%x, f_z%x, &
              u, v, w, c_xh%B, rho%x(1,1,1,1), dt, &
              ext_bdf%diffusion_coeffs, ext_bdf%ndiff, n)
      end if
 
      call ulag%update()
      call vlag%update()
      call wlag%update()
 
      call this%bc_apply_vel(t, tstep, strong = .true.)
      call this%bc_apply_prs(t, tstep)
 
      ! Update material properties if necessary
      call this%update_material_properties(t, tstep)
 
      ! Compute pressure residual.
      call profiler_start_region('Pressure_residual', 18)
 
      call prs_res%compute(p, p_res, &
           u, v, w, &
           u_e, v_e, w_e, &
           f_x, f_y, f_z, &
           c_xh, gs_xh, &
           this%bc_prs_surface, this%bc_sym_surface, &
           ax_prs, ext_bdf%diffusion_coeffs(1), dt, &
           mu_field, rho_field, event)
 
      ! De-mean the pressure residual when no strong pressure boundaries present
      if (.not. this%prs_dirichlet) call ortho(p_res%x, this%glb_n_points, n)
 
      call gs_xh%op(p_res, gs_op_add, event)
      if (neko_bcknd_device .eq. 1) then
         call device_stream_wait_event(glb_cmd_queue, event, 0)
      end if
 
      ! Set the residual to zero at strong pressure boundaries.
      call this%bclst_dp%apply_scalar(p_res%x, p%dof%size(), t, tstep)
 
 
      call profiler_end_region('Pressure_residual', 18)
 
 
      call this%proj_prs%pre_solving(p_res%x, tstep, c_xh, n, dt_controller, &
           'Pressure')
 
      call this%pc_prs%update()
 
      call profiler_start_region('Pressure_solve', 3)
 
      ! Solve for the pressure increment.
      ksp_results(1) = this%ksp_prs%solve(ax_prs, dp, p_res%x, n, c_xh, &
           this%bclst_dp, gs_xh)
 
 
      call profiler_end_region('Pressure_solve', 3)
 
      call this%proj_prs%post_solving(dp%x, ax_prs, c_xh, &
           this%bclst_dp, gs_xh, n, tstep, dt_controller)
 
      ! Update the pressure with the increment. Demean if necessary.
      call field_add2(p, dp, n)
      if (.not. this%prs_dirichlet) call ortho(p%x, this%glb_n_points, n)
 
      ! Compute velocity residual.
      call profiler_start_region('Velocity_residual', 19)
      call vel_res%compute(ax_vel, u, v, w, &
           u_res, v_res, w_res, &
           p, &
           f_x, f_y, f_z, &
           c_xh, msh, xh, &
           mu_field, rho_field, ext_bdf%diffusion_coeffs(1), &
           dt, dm_xh%size())
 
      call gs_xh%op(u_res, gs_op_add, event)
      call gs_xh%op(v_res, gs_op_add, event)
      call gs_xh%op(w_res, gs_op_add, event)
      if (neko_bcknd_device .eq. 1) then
         call device_stream_wait_event(glb_cmd_queue, event, 0)
      end if
 
      ! Set residual to zero at strong velocity boundaries.
      call this%bclst_vel_res%apply(u_res, v_res, w_res, t, tstep)
 
 
      call profiler_end_region('Velocity_residual', 19)
 
      call this%proj_u%pre_solving(u_res%x, tstep, c_xh, n, dt_controller)
      call this%proj_v%pre_solving(v_res%x, tstep, c_xh, n, dt_controller)
      call this%proj_w%pre_solving(w_res%x, tstep, c_xh, n, dt_controller)
 
      call this%pc_vel%update()
 
      call profiler_start_region("Velocity_solve", 4)
      ksp_results(2:4) = this%ksp_vel%solve_coupled(ax_vel, du, dv, dw, &
           u_res%x, v_res%x, w_res%x, n, c_xh, &
           this%bclst_du, this%bclst_dv, this%bclst_dw, gs_xh, &
           this%ksp_vel%max_iter)
      call profiler_end_region("Velocity_solve", 4)
 
      call this%proj_u%post_solving(du%x, ax_vel, c_xh, &
           this%bclst_du, gs_xh, n, tstep, dt_controller)
      call this%proj_v%post_solving(dv%x, ax_vel, c_xh, &
           this%bclst_dv, gs_xh, n, tstep, dt_controller)
      call this%proj_w%post_solving(dw%x, ax_vel, c_xh, &
           this%bclst_dw, gs_xh, n, tstep, dt_controller)
 
      if (neko_bcknd_device .eq. 1) then
         call device_opadd2cm(u%x_d, v%x_d, w%x_d, &
              du%x_d, dv%x_d, dw%x_d, 1.0_rp, n, msh%gdim)
      else
         call opadd2cm(u%x, v%x, w%x, du%x, dv%x, dw%x, 1.0_rp, n, msh%gdim)
      end if
 
      if (this%forced_flow_rate) then
         call neko_error('Forced flow rate is not implemented for the adjoint')
         !call this%vol_flow%adjust( u, v, w, p, u_res, v_res, w_res, p_res, &
         !     c_Xh, gs_Xh, ext_bdf, rho, mu, dt, &
         !     this%bclst_dp, this%bclst_du, this%bclst_dv, &
         !     this%bclst_dw, this%bclst_vel_res, Ax_vel, Ax_prs, this%ksp_prs,&
         !     this%ksp_vel, this%pc_prs, this%pc_vel, this%ksp_prs%max_iter, &
         !     this%ksp_vel%max_iter)
      end if
 
      call fluid_step_info(tstep, t, dt, ksp_results, this%strict_convergence)
 
      call this%scratch%relinquish_field(temp_indices)
 
    end associate
    call profiler_end_region('Adjoint', 1)
 
    ! Compute the norm of the field and determine if we should do a rescale.
    ! TODO: HARRY
    ! we need to discuss this rescaling in the context of topology optimization
    ! vs stability analysis,
    ! for now I'm commenting this out
    !call this%PW_compute_(t, tstep)
 
  end subroutine adjoint_fluid_pnpn_step
 
  ! ========================================================================== !
 
  subroutine adjoint_fluid_pnpn_setup_bcs(this, user, params)
    use mpi_f08, only: mpi_in_place
 
    class(adjoint_fluid_pnpn_t), intent(inout) :: this
    type(user_t), target, intent(in) :: user
    type(json_file), intent(inout) :: params
    integer :: i, n_bcs, j, zone_size, global_zone_size, ierr
    class(bc_t), pointer :: bc_i
    type(json_core) :: core
    type(json_value), pointer :: bc_object
    type(json_file) :: bc_subdict
    logical :: found
    ! Monitor which boundary zones have been marked
    logical, allocatable :: marked_zones(:)
    integer, allocatable :: zone_indices(:)
    character(len=:), allocatable :: json_key
 
    ! Lists for the residuals and solution increments
    call this%bclst_vel_res%init()
    call this%bclst_du%init()
    call this%bclst_dv%init()
    call this%bclst_dw%init()
    call this%bclst_dp%init()
 
    call this%bc_vel_res%init_from_components(this%c_Xh)
    call this%bc_du%init_from_components(this%c_Xh)
    call this%bc_dv%init_from_components(this%c_Xh)
    call this%bc_dw%init_from_components(this%c_Xh)
    call this%bc_dp%init_from_components(this%c_Xh)
 
    ! Special PnPn boundary conditions for pressure
    call this%bc_prs_surface%init_from_components(this%c_Xh)
    call this%bc_sym_surface%init_from_components(this%c_Xh)
 
    json_key = 'case.adjoint_fluid.boundary_conditions'
 
    ! Populate bcs_vel and bcs_prs based on the case file
    if (params%valid_path(json_key)) then
       call params%info(json_key, n_children = n_bcs)
       call params%get_core(core)
       call params%get(json_key, bc_object, found)
 
       !
       ! Velocity bcs
       !
       call this%bcs_vel%init(n_bcs)
 
       allocate(marked_zones(size(this%msh%labeled_zones)))
       marked_zones = .false.
 
       do i = 1, n_bcs
          ! Create a new json containing just the subdict for this bc
          call json_extract_item(core, bc_object, i, bc_subdict)
 
          call json_get(bc_subdict, "zone_indices", zone_indices)
 
          ! Check that we are not trying to assing a bc to zone, for which one
          ! has already been assigned and that the zone has more than 0 size
          ! in the mesh.
          do j = 1, size(zone_indices)
             zone_size = this%msh%labeled_zones(zone_indices(j))%size
             call mpi_allreduce(zone_size, global_zone_size, 1, &
                  mpi_integer, mpi_max, neko_comm, ierr)
 
             if (global_zone_size .eq. 0) then
                write(error_unit, '(A, A, I0, A, A, I0, A)') "*** ERROR ***: ",&
                     "Zone index ", zone_indices(j), &
                     " is invalid as this zone has 0 size, meaning it does ", &
                     "not in the mesh. Check adjoint boundary condition ", &
                     i, "."
                error stop
             end if
 
             if (marked_zones(zone_indices(j)) .eqv. .true.) then
                write(error_unit, '(A, A, I0, A, A, A, A)') "*** ERROR ***: ", &
                     "Zone with index ", zone_indices(j), &
                     " has already been assigned a boundary condition. ", &
                     "Please check your boundary_conditions entry for the ", &
                     "adjoint and make sure that each zone index appears ", &
                     "only in a single boundary condition."
                error stop
             else
                marked_zones(zone_indices(j)) = .true.
             end if
          end do
 
          bc_i => null()
          call velocity_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
 
          ! Not all bcs require an allocation for velocity in particular,
          ! so we check.
          if (associated(bc_i)) then
 
             ! We need to treat mixed bcs separately because they are by
             ! convention marked weak and currently contain nested
             ! bcs, some of which are strong.
             select type (bc_i)
             type is (symmetry_t)
                ! Symmetry has 3 internal bcs, but only one actually contains
                ! markings.
                ! Symmetry's apply_scalar doesn't do anything, so we need to
                ! mark individual nested bcs to the du,dv,dw, whereas the
                ! vel_res can just get symmetry as a whole, because on this
                ! list we call apply_vector.
                ! Additionally we have to mark the special surface bc for p.
                call this%bclst_vel_res%append(bc_i)
                call this%bc_du%mark_facets(bc_i%bc_x%marked_facet)
                call this%bc_dv%mark_facets(bc_i%bc_y%marked_facet)
                call this%bc_dw%mark_facets(bc_i%bc_z%marked_facet)
 
                call this%bcs_vel%append(bc_i)
 
                call this%bc_sym_surface%mark_facets(bc_i%marked_facet)
             type is (non_normal_t)
                ! This is a bc for the residuals and increments, not the
                ! velocity itself. So, don't append to bcs_vel
                call this%bclst_vel_res%append(bc_i)
                call this%bc_du%mark_facets(bc_i%bc_x%marked_facet)
                call this%bc_dv%mark_facets(bc_i%bc_y%marked_facet)
                call this%bc_dw%mark_facets(bc_i%bc_z%marked_facet)
             type is (shear_stress_t)
                ! Same as symmetry
                call this%bclst_vel_res%append(bc_i%symmetry)
                call this%bclst_du%append(bc_i%symmetry%bc_x)
                call this%bclst_dv%append(bc_i%symmetry%bc_y)
                call this%bclst_dw%append(bc_i%symmetry%bc_z)
 
                call this%bcs_vel%append(bc_i)
             type is (wall_model_bc_t)
                ! Same as symmetry
                call this%bclst_vel_res%append(bc_i%symmetry)
                call this%bclst_du%append(bc_i%symmetry%bc_x)
                call this%bclst_dv%append(bc_i%symmetry%bc_y)
                call this%bclst_dw%append(bc_i%symmetry%bc_z)
 
                call this%bcs_vel%append(bc_i)
             class default
 
                ! For the default case we use our dummy zero_dirichlet bcs to
                ! mark the same faces as in ordinary velocity dirichlet
                ! conditions.
                ! Additionally we mark the special PnPn pressure  bc.
                if (bc_i%strong .eqv. .true.) then
                   call this%bc_vel_res%mark_facets(bc_i%marked_facet)
                   call this%bc_du%mark_facets(bc_i%marked_facet)
                   call this%bc_dv%mark_facets(bc_i%marked_facet)
                   call this%bc_dw%mark_facets(bc_i%marked_facet)
 
                   call this%bc_prs_surface%mark_facets(bc_i%marked_facet)
                end if
 
                call this%bcs_vel%append(bc_i)
             end select
          end if
       end do
 
       ! Make sure all labeled zones with non-zero size have been marked
       do i = 1, size(this%msh%labeled_zones)
          if ((this%msh%labeled_zones(i)%size .gt. 0) .and. &
               (marked_zones(i) .eqv. .false.)) then
             write(error_unit, '(A, A, I0)') "*** ERROR ***: ", &
                  "No adjoint boundary condition assigned to zone ", i
             error stop
          end if
       end do
 
       !
       ! Pressure bcs
       !
       call this%bcs_prs%init(n_bcs)
 
       do i = 1, n_bcs
          ! Create a new json containing just the subdict for this bc
          call json_extract_item(core, bc_object, i, bc_subdict)
          bc_i => null()
          call pressure_bc_factory(bc_i, this, bc_subdict, this%c_Xh, user)
 
          ! Not all bcs require an allocation for pressure in particular,
          ! so we check.
          if (associated(bc_i)) then
             call this%bcs_prs%append(bc_i)
 
             ! Mark strong bcs in the dummy dp bc to force zero change.
             if (bc_i%strong .eqv. .true.) then
                call this%bc_dp%mark_facets(bc_i%marked_facet)
             end if
 
          end if
 
       end do
    end if
 
    call this%bc_prs_surface%finalize()
    call this%bc_sym_surface%finalize()
 
    call this%bc_vel_res%finalize()
    call this%bc_du%finalize()
    call this%bc_dv%finalize()
    call this%bc_dw%finalize()
    call this%bc_dp%finalize()
 
    call this%bclst_vel_res%append(this%bc_vel_res)
    call this%bclst_du%append(this%bc_du)
    call this%bclst_dv%append(this%bc_dv)
    call this%bclst_dw%append(this%bc_dw)
    call this%bclst_dp%append(this%bc_dp)
 
    ! If we have no strong pressure bcs, we will demean the pressure
    this%prs_dirichlet = .not. this%bclst_dp%is_empty()
    call mpi_allreduce(mpi_in_place, this%prs_dirichlet, 1, &
         mpi_logical, mpi_lor, neko_comm)
 
  end subroutine adjoint_fluid_pnpn_setup_bcs
 
  subroutine adjoint_fluid_pnpn_write_boundary_conditions(this)
    use inflow, only: inflow_t
    use field_dirichlet, only: field_dirichlet_t
    use blasius, only: blasius_t
    use field_dirichlet_vector, only: field_dirichlet_vector_t
    use usr_inflow, only: usr_inflow_t
    use dong_outflow, only: dong_outflow_t
    class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
    type(dirichlet_t) :: bdry_mask
    type(field_t), pointer :: bdry_field
    type(file_t) :: bdry_file
    integer :: temp_index, i
    class(bc_t), pointer :: bci
    character(len=LOG_SIZE) :: log_buf
 
    call neko_log%section("Adjoint boundary conditions")
    write(log_buf, '(A)') &
         'Marking using integer keys in boundary_adjoint0.f00000'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') 'Condition-value pairs: '
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  no_slip                         = 1'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  velocity_value                  = 2'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  outflow, normal_outflow (+dong) = 3'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  symmetry                        = 4'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  user_velocity_pointwise         = 5'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  periodic                        = 6'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  user_velocity                   = 7'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  user_pressure                   = 8'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  shear_stress                    = 9'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  wall_modelling                  = 10'
    call neko_log%message(log_buf)
    write(log_buf, '(A)') '  blasius_profile                 = 11'
    call neko_log%message(log_buf)
    call neko_log%end_section()
 
    call this%scratch%request_field(bdry_field, temp_index)
    bdry_field = 0.0_rp
 
 
 
    call bdry_mask%init_from_components(this%c_Xh, 6.0_rp)
    call bdry_mask%mark_zone(this%msh%periodic)
    call bdry_mask%finalize()
    call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
    call bdry_mask%free()
 
    do i = 1, this%bcs_prs%size()
       bci => this%bcs_prs%get(i)
       select type (bc => bci)
       type is (zero_dirichlet_t)
          call bdry_mask%init_from_components(this%c_Xh, 3.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (dong_outflow_t)
          call bdry_mask%init_from_components(this%c_Xh, 3.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (field_dirichlet_t)
          call bdry_mask%init_from_components(this%c_Xh, 8.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       end select
    end do
 
    do i = 1, this%bcs_vel%size()
       bci => this%bcs_vel%get(i)
       select type (bc => bci)
       type is (zero_dirichlet_t)
          call bdry_mask%init_from_components(this%c_Xh, 1.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (inflow_t)
          call bdry_mask%init_from_components(this%c_Xh, 2.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (symmetry_t)
          call bdry_mask%init_from_components(this%c_Xh, 4.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (usr_inflow_t)
          call bdry_mask%init_from_components(this%c_Xh, 5.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (field_dirichlet_vector_t)
          call bdry_mask%init_from_components(this%c_Xh, 7.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (shear_stress_t)
          call bdry_mask%init_from_components(this%c_Xh, 9.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (wall_model_bc_t)
          call bdry_mask%init_from_components(this%c_Xh, 10.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       type is (blasius_t)
          call bdry_mask%init_from_components(this%c_Xh, 11.0_rp)
          call bdry_mask%mark_facets(bci%marked_facet)
          call bdry_mask%finalize()
          call bdry_mask%apply_scalar(bdry_field%x, this%dm_Xh%size())
          call bdry_mask%free()
       end select
    end do
 
 
    bdry_file = file_t('boundary_adjoint.fld')
    call bdry_file%write(bdry_field)
 
    call this%scratch%relinquish_field(temp_index)
  end subroutine adjoint_fluid_pnpn_write_boundary_conditions
 
  ! End of section to verify
  ! ========================================================================== !
 
  subroutine rescale_fluid(fluid_data, scale)
 
    class(adjoint_fluid_pnpn_t), intent(inout) :: fluid_data
    real(kind=rp), intent(in) :: scale
 
    ! Local variables
    integer :: i
 
    ! Scale the velocity fields
    if (neko_bcknd_device .eq. 1) then
       call device_cmult(fluid_data%u_adj%x_d, scale, fluid_data%u_adj%size())
       call device_cmult(fluid_data%v_adj%x_d, scale, fluid_data%v_adj%size())
       call device_cmult(fluid_data%w_adj%x_d, scale, fluid_data%w_adj%size())
    else
       call cmult(fluid_data%u_adj%x, scale, fluid_data%u_adj%size())
       call cmult(fluid_data%v_adj%x, scale, fluid_data%v_adj%size())
       call cmult(fluid_data%w_adj%x, scale, fluid_data%w_adj%size())
    end if
 
    ! Scale the right hand sides
    if (neko_bcknd_device .eq. 1) then
       call device_cmult(fluid_data%f_adj_x%x_d, scale, &
            fluid_data%f_adj_x%size())
       call device_cmult(fluid_data%f_adj_y%x_d, scale, &
            fluid_data%f_adj_y%size())
       call device_cmult(fluid_data%f_adj_z%x_d, scale, &
            fluid_data%f_adj_z%size())
       ! HARRY
       ! maybe the abx's too
       call device_cmult(fluid_data%abx1%x_d, scale, fluid_data%abx1%size())
       call device_cmult(fluid_data%aby1%x_d, scale, fluid_data%aby1%size())
       call device_cmult(fluid_data%abz1%x_d, scale, fluid_data%abz1%size())
       call device_cmult(fluid_data%abx2%x_d, scale, fluid_data%abx2%size())
       call device_cmult(fluid_data%aby2%x_d, scale, fluid_data%aby2%size())
       call device_cmult(fluid_data%abz2%x_d, scale, fluid_data%abz2%size())
 
    else
       call cmult(fluid_data%f_adj_x%x, scale, fluid_data%f_adj_x%size())
       call cmult(fluid_data%f_adj_y%x, scale, fluid_data%f_adj_y%size())
       call cmult(fluid_data%f_adj_z%x, scale, fluid_data%f_adj_z%size())
 
       call cmult(fluid_data%abx1%x, scale, fluid_data%abx1%size())
       call cmult(fluid_data%aby1%x, scale, fluid_data%aby1%size())
       call cmult(fluid_data%abz1%x, scale, fluid_data%abz1%size())
 
       call cmult(fluid_data%abx2%x, scale, fluid_data%abx2%size())
       call cmult(fluid_data%aby2%x, scale, fluid_data%aby2%size())
       call cmult(fluid_data%abz2%x, scale, fluid_data%abz2%size())
    end if
 
    ! Scale the lag terms
    if (neko_bcknd_device .eq. 1) then
       do i = 1, fluid_data%ulag%size()
          call device_cmult(fluid_data%ulag%lf(i)%x_d, &
               scale, fluid_data%ulag%lf(i)%size())
       end do
 
       do i = 1, fluid_data%vlag%size()
          call device_cmult(fluid_data%vlag%lf(i)%x_d, &
               scale, fluid_data%vlag%lf(i)%size())
       end do
 
       do i = 1, fluid_data%wlag%size()
          call device_cmult(fluid_data%wlag%lf(i)%x_d, &
               scale, fluid_data%wlag%lf(i)%size())
       end do
    else
       do i = 1, fluid_data%ulag%size()
          call cmult(fluid_data%ulag%lf(i)%x, &
               scale, fluid_data%ulag%lf(i)%size())
       end do
 
       do i = 1, fluid_data%vlag%size()
          call cmult(fluid_data%vlag%lf(i)%x, &
               scale, fluid_data%vlag%lf(i)%size())
       end do
 
       do i = 1, fluid_data%wlag%size()
          call cmult(fluid_data%wlag%lf(i)%x, &
               scale, fluid_data%wlag%lf(i)%size())
       end do
    end if
 
  end subroutine rescale_fluid
 
  function norm(x, y, z, B, volume, n)
    use mpi_f08, only: mpi_in_place
 
    integer, intent(in) :: n
    real(kind=rp), dimension(n), intent(in) :: x, y, z
    real(kind=rp), dimension(n), intent(in) :: b
    real(kind=rp), intent(in) :: volume
 
    real(kind=rp) :: norm
 
    norm = vlsc3(x, x, b, n) + vlsc3(y, y, b, n) + vlsc3(z, z, b, n)
 
    call mpi_allreduce(mpi_in_place, norm, 1, &
         mpi_real_precision, mpi_sum, mpi_comm_world)
 
    norm = sqrt(norm / volume)
  end function norm
 
  function device_norm(x_d, y_d, z_d, B_d, volume, n)
    use mpi_f08, only: mpi_in_place
 
    type(c_ptr), intent(in) :: x_d, y_d, z_d
    type(c_ptr), intent(in) :: B_d
    real(kind=rp), intent(in) :: volume
    integer, intent(in) :: n
 
    real(kind=rp) :: device_norm
 
    device_norm = device_vlsc3(x_d, x_d, b_d, n) + &
         device_vlsc3(y_d, y_d, b_d, n) + &
         device_vlsc3(z_d, z_d, b_d, n)
 
    call mpi_allreduce(mpi_in_place, device_norm, 1, &
         mpi_real_precision, mpi_sum, mpi_comm_world)
 
    device_norm = sqrt(device_norm / volume)
 
  end function device_norm
 
  subroutine power_iterations_compute(this, t, tstep)
    class(adjoint_fluid_pnpn_t), target, intent(inout) :: this
    real(kind=rp), intent(in) :: t
    integer, intent(in) :: tstep
 
    ! Local variables
    real(kind=rp) :: scaling_factor
    real(kind=rp) :: norm_l2, norm_l2_base
    character(len=256) :: log_message
    type(vector_t) :: data_line
    integer :: n
 
    n = this%c_Xh%dof%size()
    if (tstep .eq. 1) then
       if (neko_bcknd_device .eq. 1) then
          norm_l2_base = device_norm(this%u_adj%x_d, this%v_adj%x_d, &
               this%w_adj%x_d, &
               this%c_Xh%B_d, this%c_Xh%volume, n)
       else
          norm_l2_base = this%norm_scaling * norm(this%u_adj%x, this%v_adj%x, &
               this%w_adj%x, &
               this%c_Xh%B, this%c_Xh%volume, n)
       end if
       if (this%norm_target .lt. 0.0_rp) then
          this%norm_target = norm_l2_base
       end if
 
       this%norm_l2_upper = this%norm_tolerance * this%norm_target
       this%norm_l2_lower = this%norm_target / this%norm_tolerance
 
    end if
 
    ! Compute the norm of the velocity field and eigenvalue estimate
    if (neko_bcknd_device .eq. 1) then
       norm_l2 = device_norm(this%u_adj%x_d, this%v_adj%x_d, this%w_adj%x_d, &
            this%c_Xh%B_d, this%c_Xh%volume, n)
    else
       norm_l2 = norm(this%u_adj%x, this%v_adj%x, this%w_adj%x, &
            this%c_Xh%B, this%c_Xh%volume, n)
    end if
    norm_l2 = sqrt(this%norm_scaling) * norm_l2
    scaling_factor = 1.0_rp
 
    ! Rescale the flow if necessary
    if (norm_l2 .gt. this%norm_l2_upper &
         .or. norm_l2 .lt. this%norm_l2_lower) then
       scaling_factor = this%norm_target / norm_l2
       call rescale_fluid(this, scaling_factor)
       norm_l2 = this%norm_target
 
       if (tstep .eq. 1) then
          scaling_factor = 1.0_rp
       end if
    end if
 
    ! Log the results
    !call neko_log%section('Power Iterations', lvl = NEKO_LOG_DEBUG)
    call neko_log%section('Power Iterations')
 
    write (log_message, '(A7,E20.14)') 'Norm: ', norm_l2
    call neko_log%message(log_message, lvl = neko_log_debug)
    write (log_message, '(A7,E20.14)') 'Scaling: ', scaling_factor
    call neko_log%message(log_message, lvl = neko_log_debug)
 
    ! Save to file
    call data_line%init(2)
    data_line%x = [norm_l2, scaling_factor]
    call this%file_output%write(data_line, t)
 
    !call neko_log%end_section('Power Iterations', lvl = NEKO_LOG_DEBUG)
    call neko_log%end_section('Power Iterations')
  end subroutine power_iterations_compute
 
 
end module adjoint_fluid_pnpn