adjoint_scalar_scheme.f90

! Copyright (c) 2022-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_scalar_scheme
  use gather_scatter, only : gs_t
  use checkpoint, only : chkp_t
  use num_types, only: rp
  use field, only : field_t
  use field_list, only: field_list_t
  use space, only : space_t
  use dofmap, only : dofmap_t
  use krylov, only : ksp_t, krylov_solver_factory, ksp_max_iter
  use coefs, only : coef_t
  use dirichlet, only : dirichlet_t
  use neumann, only : neumann_t
  use jacobi, only : jacobi_t
  use device_jacobi, only : device_jacobi_t
  use sx_jacobi, only : sx_jacobi_t
  use hsmg, only : hsmg_t
  use bc, only : bc_t
  use bc_list, only : bc_list_t
  use precon, only : pc_t, precon_factory, precon_destroy
  use field_dirichlet, only: field_dirichlet_t, field_dirichlet_update
  use mesh, only : mesh_t, neko_msh_max_zlbls, neko_msh_max_zlbl_len
  use facet_zone, only : facet_zone_t
  use time_scheme_controller, only : time_scheme_controller_t
  use logger, only : neko_log, log_size, neko_log_verbose
  use field_registry, only : neko_field_registry
  use usr_scalar, only : usr_scalar_t, usr_scalar_bc_eval
  use json_utils, only : json_get, json_get_or_default, json_extract_item, &
       json_extract_object
  use json_module, only : json_file
  use user_intf, only : user_t, dummy_user_material_properties, &
       user_material_properties
  use utils, only : neko_error, neko_warning
  use comm, only: neko_comm, mpi_integer, mpi_sum
  use scalar_source_term, only : scalar_source_term_t
  use field_series, only : field_series_t
  use math, only : cfill, add2s2
  use field_math, only : field_add2s2
  use device_math, only : device_cfill, device_add2s2
  use neko_config, only : neko_bcknd_device
  use field_series, only : field_series_t
  use time_step_controller, only : time_step_controller_t
  use gradient_jump_penalty, only : gradient_jump_penalty_t
  use json_utils_ext, only: json_key_fallback
  use scalar_scheme, only: scalar_scheme_precon_factory, &
       scalar_scheme_solver_factory
  use scratch_registry, only : neko_scratch_registry
  use time_state, only : time_state_t
  use device, only : device_memcpy, device_to_host
  use field_math, only : field_col3, field_cmult2, field_add2, field_cfill
  implicit none
 
  type, abstract :: adjoint_scalar_scheme_t
     character(len=:), allocatable :: name
     type(field_t), pointer :: u
     type(field_t), pointer :: v
     type(field_t), pointer :: w
     type(field_t), pointer :: s
     type(field_t), pointer :: s_adj
     type(field_series_t) :: s_adj_lag
     type(space_t), pointer :: xh
     type(dofmap_t), pointer :: dm_xh
     type(gs_t), pointer :: gs_xh
     type(coef_t), pointer :: c_xh
     type(field_t), pointer :: f_xh => null()
     type(scalar_source_term_t) :: source_term
     class(ksp_t), allocatable :: ksp
     integer :: ksp_maxiter
     integer :: projection_dim
     
     integer :: projection_activ_step
     class(pc_t), allocatable :: pc
     type(bc_list_t) :: bcs
     type(json_file), pointer :: params
     type(mesh_t), pointer :: msh => null()
     type(chkp_t), pointer :: chkp => null()
     character(len=:), allocatable :: nut_field_name
     type(field_t), pointer :: rho => null()
     type(field_t) :: lambda
     type(field_t) :: cp
     real(kind=rp) :: pr_turb
     type(field_list_t) :: material_properties
     logical :: variable_material_properties = .false.
     procedure(user_material_properties), nopass, pointer :: &
          user_material_properties => null()
     logical :: if_gradient_jump_penalty
     type(gradient_jump_penalty_t) :: gradient_jump_penalty
   contains
     procedure, pass(this) :: scheme_init => adjoint_scalar_scheme_init
     procedure, pass(this) :: scheme_free => adjoint_scalar_scheme_free
     procedure, pass(this) :: validate => adjoint_scalar_scheme_validate
     procedure, pass(this) :: set_material_properties => &
          adjoint_scalar_scheme_set_material_properties
     procedure, pass(this) :: update_material_properties => &
          adjoint_scalar_scheme_update_material_properties
     procedure(adjoint_scalar_scheme_init_intrf), pass(this), deferred :: init
     procedure(adjoint_scalar_scheme_free_intrf), pass(this), deferred :: free
     procedure(adjoint_scalar_scheme_step_intrf), pass(this), deferred :: step
     procedure(adjoint_scalar_scheme_restart_intrf), pass(this), deferred :: &
          restart
  end type adjoint_scalar_scheme_t
 
  abstract interface
     subroutine adjoint_scalar_scheme_init_intrf(this, msh, coef, gs, params, &
          user, ulag, vlag, wlag, time_scheme, rho)
       import adjoint_scalar_scheme_t
       import json_file
       import coef_t
       import gs_t
       import mesh_t
       import user_t
       import field_series_t, field_t
       import time_scheme_controller_t
       import rp
       class(adjoint_scalar_scheme_t), target, intent(inout) :: this
       type(mesh_t), target, intent(in) :: msh
       type(coef_t), target, intent(in) :: coef
       type(gs_t), target, intent(inout) :: gs
       type(json_file), target, intent(inout) :: params
       type(user_t), target, intent(in) :: user
       type(field_series_t), target, intent(in) :: ulag, vlag, wlag
       type(time_scheme_controller_t), target, intent(in) :: time_scheme
       type(field_t), target, intent(in) :: rho
     end subroutine adjoint_scalar_scheme_init_intrf
  end interface
 
  abstract interface
     subroutine adjoint_scalar_scheme_restart_intrf(this, dtlag, tlag)
       import adjoint_scalar_scheme_t
       import chkp_t
       import rp
       class(adjoint_scalar_scheme_t), target, intent(inout) :: this
       real(kind=rp) :: dtlag(10), tlag(10)
     end subroutine adjoint_scalar_scheme_restart_intrf
  end interface
 
  abstract interface
     subroutine adjoint_scalar_scheme_free_intrf(this)
       import adjoint_scalar_scheme_t
       class(adjoint_scalar_scheme_t), intent(inout) :: this
     end subroutine adjoint_scalar_scheme_free_intrf
  end interface
 
  abstract interface
     subroutine adjoint_scalar_scheme_step_intrf(this, t, tstep, dt, ext_bdf, &
          dt_controller)
       import adjoint_scalar_scheme_t
       import time_scheme_controller_t
       import time_step_controller_t
       import rp
       class(adjoint_scalar_scheme_t), 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
     end subroutine adjoint_scalar_scheme_step_intrf
  end interface
 
contains
 
  subroutine adjoint_scalar_scheme_init(this, msh, c_Xh, gs_Xh, params, &
       scheme, user, rho)
    class(adjoint_scalar_scheme_t), target, intent(inout) :: this
    type(mesh_t), target, intent(in) :: msh
    type(coef_t), target, intent(in) :: c_Xh
    type(gs_t), target, intent(inout) :: gs_Xh
    type(json_file), target, intent(inout) :: params
    character(len=*), intent(in) :: scheme
    type(user_t), target, intent(in) :: user
    type(field_t), target, intent(in) :: rho
    ! IO buffer for log output
    character(len=LOG_SIZE) :: log_buf
    ! Variables for retrieving json parameters
    logical :: logical_val
    real(kind=rp) :: solver_abstol
    integer :: integer_val
    character(len=:), allocatable :: solver_type, solver_precon
    ! real(kind=rp) :: GJP_param_a, GJP_param_b
    character(len=:), allocatable :: json_key
    type(json_file) :: precon_params
 
    this%u => neko_field_registry%get_field('u')
    this%v => neko_field_registry%get_field('v')
    this%w => neko_field_registry%get_field('w')
    this%s => neko_field_registry%get_field('s')
    this%rho => rho
 
    call neko_log%section('Adjoint scalar')
    ! In principle this could be different,
    ! but I think we should use the same as the forward
    call json_get(params, 'case.scalar.solver.type', solver_type)
    json_key = json_key_fallback(params, &
         'case.adjoint_scalar.solver.preconditioner.type', &
         'case.scalar.solver.preconditioner.type')
    call json_get(params, json_key, solver_precon)
 
    json_key = json_key_fallback(params, &
         'case.adjoint_scalar.solver.preconditioner', &
         'case.scalar.solver.preconditioner')
    call json_extract_object(params, json_key, precon_params)
 
    json_key = json_key_fallback(params, &
         'case.adjoint_scalar.solver.absolute_tolerance', &
         'case.scalar.solver.absolute_tolerance')
    call json_get(params, json_key, solver_abstol)
 
    json_key = json_key_fallback(params, &
         'case.adjoint_scalar.solver.projection_space_size', &
         'case.scalar.solver.projection_space_size')
    call json_get_or_default(params, json_key, this%projection_dim, 0)
 
    json_key = json_key_fallback(params, &
         'case.adjoint_scalar.solver.projection_hold_steps', &
         'case.scalar.solver.projection_hold_steps')
    call json_get_or_default(params, json_key, this%projection_activ_step, 5)
 
 
    write(log_buf, '(A, A)') 'Type       : ', trim(scheme)
    call neko_log%message(log_buf)
    call neko_log%message('Ksp adjoint scalar : ('// trim(solver_type) // &
         ', ' // trim(solver_precon) // ')')
    write(log_buf, '(A,ES13.6)') ' `-abs tol :', solver_abstol
    call neko_log%message(log_buf)
 
    this%Xh => this%u%Xh
    this%dm_Xh => this%u%dof
    this%params => params
    this%msh => msh
    if (.not. neko_field_registry%field_exists('s_adj')) then
       call neko_field_registry%add_field(this%dm_Xh, 's_adj')
    end if
    this%s_adj => neko_field_registry%get_field('s_adj')
 
    call this%s_adj_lag%init(this%s_adj, 2)
 
    this%gs_Xh => gs_xh
    this%c_Xh => c_xh
 
    !
    ! Material properties
    !
    call this%set_material_properties(params, user)
    !
    ! Turbulence modelling and variable material properties
    !
    if (params%valid_path('case.scalar.variable_material_properties')) then
       call json_get(params, 'variable_material_properties', &
            this%variable_material_properties)
 
       ! Warn, no variable properties, but nut_field
       if ((params%valid_path('case.scalar.nut_field')) .and. &
            (this%variable_material_properties .eqv. .false.)) then
          call neko_warning("You set variable_material properties to " // &
               "false, the nut_field setting will have no effect.")
       end if
 
       ! Warn, no variable properties, but user routine associated
       if ((.not. associated(user%material_properties, &
            dummy_user_material_properties)) .and. &
            (this%variable_material_properties .eqv. .false.)) then
          call neko_warning("You set variable_material properties to " // &
               "false, you can only vary rho and mu in time in the user file.")
       end if
    else if (params%valid_path('case.scalar.nut_field')) then
       call json_get(params, 'case.scalar.Pr_t', this%pr_turb)
       call json_get(params, 'case.scalar.nut_field', this%nut_field_name)
       this%variable_material_properties = .true.
    else if (.not. associated(user%material_properties, &
         dummy_user_material_properties)) then
       this%nut_field_name = ""
       this%variable_material_properties = .true.
    end if
 
    write(log_buf, '(A,L1)') 'LES        : ', this%variable_material_properties
    call neko_log%message(log_buf)
 
 
    !
    ! Setup right-hand side field.
    !
    allocate(this%f_Xh)
    call this%f_Xh%init(this%dm_Xh, fld_name = "adjoint_scalar_rhs")
 
    ! Initialize the source term
    call this%source_term%init(this%f_Xh, this%c_Xh, user)
    ! We should ONLY read the adjoint
    call this%source_term%add(params, 'case.adjoint_scalar.source_terms')
 
    ! todo parameter file ksp tol should be added
    json_key = json_key_fallback(params, &
         'case.adjoint_fluid.velocity_solver.max_iterations', &
         'case.fluid.velocity_solver.max_iterations')
    call json_get_or_default(params, json_key, integer_val, ksp_max_iter)
    json_key = json_key_fallback(params, &
         'case.adjoint_fluid.velocity_solver.monitor', &
         'case.fluid.velocity_solver.monitor')
    call json_get_or_default(params, json_key, logical_val, .false.)
    call scalar_scheme_solver_factory(this%ksp, this%dm_Xh%size(), &
         solver_type, integer_val, solver_abstol, logical_val)
    call scalar_scheme_precon_factory(this%pc, this%ksp, &
         this%c_Xh, this%dm_Xh, this%gs_Xh, this%bcs, &
         solver_precon, precon_params)
 
    ! Initiate gradient jump penalty
    json_key = json_key_fallback(params, &
         'case.adjoint_scalar.gradient_jump_penalty.enabled', &
         'case.scalar.gradient_jump_penalty.enabled')
    call json_get_or_default(params, json_key, &
         this%if_gradient_jump_penalty, .false.)
 
    if (this%if_gradient_jump_penalty .eqv. .true.) then
       call neko_error("gradient jump not implemented for adjoint scalar")
       ! if ((this%dm_Xh%xh%lx - 1) .eq. 1) then
       !    call json_get_or_default(params, &
       !         'case.scalar.gradient_jump_penalty.tau',&
       !         GJP_param_a, 0.02_rp)
       !    GJP_param_b = 0.0_rp
       ! else
       !    call json_get_or_default(params, &
       !         'case.scalar.gradient_jump_penalty.scaling_factor',&
       !         GJP_param_a, 0.8_rp)
       !    call json_get_or_default(params, &
       !         'case.scalar.gradient_jump_penalty.scaling_exponent',&
       !         GJP_param_b, 4.0_rp)
       ! end if
       ! call this%gradient_jump_penalty%init(params, this%dm_Xh, this%c_Xh, &
       !      GJP_param_a, GJP_param_b)
    end if
 
    call neko_log%end_section()
 
  end subroutine adjoint_scalar_scheme_init
 
 
  subroutine adjoint_scalar_scheme_free(this)
    class(adjoint_scalar_scheme_t), intent(inout) :: this
 
    nullify(this%Xh)
    nullify(this%dm_Xh)
    nullify(this%gs_Xh)
    nullify(this%c_Xh)
    nullify(this%params)
 
    if (allocated(this%ksp)) then
       call this%ksp%free()
       deallocate(this%ksp)
    end if
 
    if (allocated(this%pc)) then
       call precon_destroy(this%pc)
       deallocate(this%pc)
    end if
 
    call this%source_term%free()
 
    call this%bcs%free()
 
    call this%cp%free()
    call this%lambda%free()
 
    ! Free gradient jump penalty
    if (this%if_gradient_jump_penalty .eqv. .true.) then
       call this%gradient_jump_penalty%free()
    end if
 
  end subroutine adjoint_scalar_scheme_free
 
  subroutine adjoint_scalar_scheme_validate(this)
    class(adjoint_scalar_scheme_t), target, intent(inout) :: this
 
    if ( (.not. allocated(this%u%x)) .or. &
         (.not. allocated(this%v%x)) .or. &
         (.not. allocated(this%w%x)) .or. &
         (.not. allocated(this%s%x)) .or. &
         (.not. allocated(this%s_adj%x))) then
       call neko_error('Fields are not allocated')
    end if
 
    if (.not. allocated(this%ksp)) then
       call neko_error('No Krylov solver for velocity defined')
    end if
 
    if (.not. associated(this%Xh)) then
       call neko_error('No function space defined')
    end if
 
    if (.not. associated(this%dm_Xh)) then
       call neko_error('No dofmap defined')
    end if
 
    if (.not. associated(this%c_Xh)) then
       call neko_error('No coefficients defined')
    end if
 
    if (.not. associated(this%f_Xh)) then
       call neko_error('No rhs allocated')
    end if
 
    if (.not. associated(this%params)) then
       call neko_error('No parameters defined')
    end if
 
    !
    ! Setup checkpoint structure (if everything is fine)
    !
!    @todo no io for now
!    call this%chkp%init(this%u, this%v, this%w, this%p)
 
  end subroutine adjoint_scalar_scheme_validate
 
  subroutine adjoint_scalar_scheme_update_material_properties(t, tstep, this)
    class(adjoint_scalar_scheme_t), intent(inout) :: this
    real(kind=rp),intent(in) :: t
    integer, intent(in) :: tstep
    type(field_t), pointer :: nut
    integer :: index
    ! Factor to transform nu_t to lambda_t
    type(field_t), pointer :: lambda_factor
 
    call this%user_material_properties(t, tstep, this%name, &
         this%material_properties)
 
    ! factor = rho * cp / pr_turb
    if (this%variable_material_properties .and. &
         len(trim(this%nut_field_name)) > 0) then
       nut => neko_field_registry%get_field(this%nut_field_name)
 
       ! lambda = lambda + rho * cp * nut / pr_turb
       call neko_scratch_registry%request_field(lambda_factor, index)
 
       call field_col3(lambda_factor, this%cp, this%rho)
       call field_cmult2(lambda_factor, nut, 1.0_rp / this%pr_turb)
       call field_add2(this%lambda, lambda_factor)
       call neko_scratch_registry%relinquish_field(index)
    end if
 
    ! Since cp is a fields and we use the %x(1,1,1,1) of the
    ! host array data to pass constant  material properties
    ! to some routines, we need to make sure that the host
    ! values are also filled
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%cp%x, this%cp%x_d, this%cp%size(), &
            device_to_host, sync=.false.)
    end if
 
  end subroutine adjoint_scalar_scheme_update_material_properties
 
  subroutine adjoint_scalar_scheme_set_material_properties(this, params, user)
    class(adjoint_scalar_scheme_t), intent(inout) :: this
    type(json_file), intent(inout) :: params
    type(user_t), target, intent(in) :: user
    character(len=LOG_SIZE) :: log_buf
    ! A local pointer that is needed to make Intel happy
    procedure(user_material_properties), pointer :: dummy_mp_ptr
    real(kind=rp) :: const_cp, const_lambda
 
    ! TODO
    ! I think everything should be read from .scalar. and we shouldnt
    ! change anything, but this should be come back to.
    dummy_mp_ptr => dummy_user_material_properties
 
    ! Fill lambda field with the physical value
    call this%lambda%init(this%dm_Xh, "lambda")
    call this%cp%init(this%dm_Xh, "cp")
 
    call this%material_properties%init(2)
    call this%material_properties%assign_to_field(1, this%cp)
    call this%material_properties%assign_to_field(2, this%lambda)
 
    if (.not. associated(user%material_properties, dummy_mp_ptr)) then
 
       write(log_buf, '(A)') "Material properties must be set in the user " // &
            "file!"
       call neko_log%message(log_buf)
       this%user_material_properties => user%material_properties
       call user%material_properties(0.0_rp, 0, this%name, &
            this%material_properties)
    else
       this%user_material_properties => dummy_user_material_properties
       if (params%valid_path('case.scalar.Pe') .and. &
            (params%valid_path('case.scalar.lambda') .or. &
            params%valid_path('case.scalar.cp'))) then
          call neko_error("To set the material properties for the scalar, " // &
               "either provide Pe OR lambda and cp in the case file.")
          ! Non-dimensional case
       else if (params%valid_path('case.scalar.Pe')) then
          write(log_buf, '(A)') 'Non-dimensional scalar material properties' //&
               ' input.'
          call neko_log%message(log_buf, lvl = neko_log_verbose)
          write(log_buf, '(A)') 'Specific heat capacity will be set to 1,'
          call neko_log%message(log_buf, lvl = neko_log_verbose)
          write(log_buf, '(A)') 'conductivity to 1/Pe. Assumes density is 1.'
          call neko_log%message(log_buf, lvl = neko_log_verbose)
 
          ! Read Pe into lambda for further manipulation.
          call json_get(params, 'case.scalar.Pe', const_lambda)
          write(log_buf, '(A,ES13.6)') 'Pe         :', const_lambda
          call neko_log%message(log_buf)
 
          ! Set cp and rho to 1 since the setup is non-dimensional.
          const_cp = 1.0_rp
          ! Invert the Pe to get conductivity
          const_lambda = 1.0_rp/const_lambda
          ! Dimensional case
       else
          call json_get(params, 'case.scalar.lambda', const_lambda)
          call json_get(params, 'case.scalar.cp', const_cp)
       end if
    end if
    ! We need to fill the fields based on the parsed const values
    ! if the user routine is not used.
    if (associated(user%material_properties, dummy_mp_ptr)) then
       ! Fill mu and rho field with the physical value
       call field_cfill(this%lambda, const_lambda)
       call field_cfill(this%cp, const_cp)
 
       write(log_buf, '(A,ES13.6)') 'lambda     :', const_lambda
       call neko_log%message(log_buf)
       write(log_buf, '(A,ES13.6)') 'cp         :', const_cp
       call neko_log%message(log_buf)
    end if
 
    ! Since cp is a field and we use the %x(1,1,1,1) of the
    ! host array data to pass constant material properties
    ! to some routines, we need to make sure that the host
    ! values are also filled
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%cp%x, this%cp%x_d, this%cp%size(), &
            device_to_host, sync=.false.)
    end if
  end subroutine adjoint_scalar_scheme_set_material_properties
 
end module adjoint_scalar_scheme