brinkman_source_term.f90

Go to the documentation of this file.

! Copyright (c) 2023-2026, 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 brinkman_source_term
  use aabb, only : aabb_t, get_aabb
  use coefs, only : coef_t
  use device, only : device_memcpy, host_to_device, device_to_host
  use field, only : field_t
  use field_list, only : field_list_t
  use math, only : cfill_mask, copy
  use device_math, only : device_cfill_mask
  use field_math, only : field_pwmax2, field_subcol3, field_copy
  use registry, only : neko_registry
  use scratch_registry, only : neko_scratch_registry
  use mappings, only : smooth_step_field, step_function_field, &
       permeability_field
  use file, only : file_t
  use json_module, only : json_file, json_core, json_value
  use json_utils, only : json_get, json_get_or_default, json_extract_item, &
       json_get_or_lookup, json_get_or_lookup_or_default
  use logger, only : neko_log, log_size, neko_log_debug
  use tri_mesh, only : tri_mesh_t
  use neko_config, only : neko_bcknd_device
  use num_types, only : rp, dp, sp
  use point_zone, only : point_zone_t
  use point_zone_registry, only : neko_point_zone_registry
  use profiler, only : profiler_start_region, profiler_end_region
  use signed_distance, only : signed_distance_field
  use source_term, only : source_term_t
  use utils, only : neko_error, filename_suffix
  use filter, only : filter_t
  use pde_filter, only : pde_filter_t
  use field_output, only : field_output_t
  use fld_file_output, only : fld_file_output_t
  use fld_file_data, only : fld_file_data_t
  use time_state, only : time_state_t
 
  use global_interpolation, only: global_interpolation_t
  use interpolation, only: interpolator_t
  use space, only: space_t, gll
  implicit none
  private
 
  type, public, extends(source_term_t) :: brinkman_source_term_t
     private
 
     type(field_t), pointer :: unfiltered => null()
     type(field_t), pointer :: indicator => null()
     type(field_t), pointer :: brinkman => null()
     class(filter_t), allocatable :: filter
 
   contains
     procedure, public, pass(this) :: init => &
          brinkman_source_term_init_from_json
     procedure, public, pass(this) :: free => brinkman_source_term_free
     procedure, public, pass(this) :: compute_ => brinkman_source_term_compute
 
     ! ----------------------------------------------------------------------- !
     ! Internal subroutines for adding a given object type to the Brinkman
     ! source term
 
     procedure, pass(this) :: add_boundary_mesh
     procedure, pass(this) :: add_point_zone
     procedure, pass(this) :: add_file
     procedure, pass(this) :: add_field
 
  end type brinkman_source_term_t
 
contains
 
  ! ========================================================================== !
  ! Public methods
 
  subroutine brinkman_source_term_init_from_json(this, json, fields, coef, &
       variable_name)
    class(brinkman_source_term_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    type(field_list_t), intent(in), target :: fields
    type(coef_t), intent(in), target :: coef
    character(len=*), intent(in) :: variable_name
    real(kind=rp) :: start_time, end_time
    character(len=LOG_SIZE) :: error_msg
 
    character(len=:), allocatable :: filter_type
    real(kind=rp), dimension(:), allocatable :: brinkman_limits
    real(kind=rp) :: brinkman_penalty
 
    character(len=:), allocatable :: object_type
    type(json_file) :: object_settings, filter_settings
    integer :: n_objects, i
 
    type(file_t) :: output
    type(field_list_t) :: output_fields
    logical :: output_enable
    character(len=:), allocatable :: output_path, output_format, &
         output_precision, fname, suffix
    integer :: precision
 
    ! ------------------------------------------------------------------------ !
    ! Read the general options for the Brinkman source term
 
    call neko_log%section('Brinkman source term')
 
    ! Read the options for the permeability field
    call json_get_or_lookup(json, 'limits', brinkman_limits)
    call json_get_or_lookup(json, 'penalty', brinkman_penalty)
 
    ! Mandatory fields for the general source term
    call json_get_or_lookup_or_default(json, "start_time", start_time, 0.0_rp)
    call json_get_or_lookup_or_default(json, "end_time", end_time, huge(0.0_rp))
 
    ! Output settings
    call json_get_or_default(json, 'output.enable', output_enable, .false.)
    call json_get_or_default(json, 'output.precision', output_precision, 'sp')
    call json_get_or_default(json, 'output.path', output_path, './')
    call json_get_or_default(json, 'output.format', output_format, 'fld')
 
    if (size(brinkman_limits) .ne. 2) then
       call neko_error('brinkman_limits must be a 2 element array of reals')
    end if
 
    if (output_path(len_trim(output_path):len_trim(output_path)) .ne. '/') then
       output_path = trim(output_path) // '/'
    end if
 
    call this%free()
    call this%init_base(fields, coef, start_time, end_time)
 
    ! ------------------------------------------------------------------------ !
    ! Allocate the permeability and indicator field
 
    call neko_registry%add_field(coef%dof, 'brinkman_indicator', .true.)
    call neko_registry%add_field(coef%dof, 'brinkman_permeability', .true.)
 
    this%indicator => neko_registry%get_field('brinkman_indicator')
    this%brinkman => neko_registry%get_field('brinkman_permeability')
 
    ! ------------------------------------------------------------------------ !
    ! Select which constructor should be called
 
    call json%info('objects', n_children = n_objects)
 
    do i = 1, n_objects
       call json_extract_item(json, "objects", i, object_settings)
       call json_get(object_settings, 'type', object_type)
 
       select case (object_type)
       case ('boundary_mesh')
          call this%add_boundary_mesh(object_settings)
       case ('point_zone')
          call this%add_point_zone(object_settings)
       case ('field')
          call this%add_field(object_settings)
       case ('file')
          call this%add_file(object_settings)
 
       case default
          write(error_msg, '(A,I0,A,A)') 'Brinkman object: ', i, &
               ' unknown type: ', trim(object_type)
          call neko_error(trim(error_msg))
       end select
 
    end do
 
    ! ------------------------------------------------------------------------ !
    ! Filter the indicator field
 
    call json_get_or_default(json, 'filter.type', filter_type, 'none')
    select case (filter_type)
    case ('PDE')
 
       ! Copy the current indicator to unfiltered (essentially a rename)
       call neko_registry%add_field(coef%dof, 'brinkman_unfiltered', .true.)
       this%unfiltered => neko_registry%get_field('brinkman_unfiltered')
       call field_copy(this%unfiltered, this%indicator)
 
       ! Allocate a PDE filter
       allocate(pde_filter_t::this%filter)
       call json_get(json, 'filter', filter_settings)
       call this%filter%init(filter_settings, coef)
       call this%filter%apply(this%indicator, this%unfiltered)
 
    case ('none')
       ! Do nothing
    case default
       call neko_error('Brinkman source term unknown filter type')
    end select
 
    ! ------------------------------------------------------------------------ !
    ! Compute the permeability field
 
    this%brinkman = this%indicator
    call permeability_field(this%brinkman, &
         brinkman_limits(1), brinkman_limits(2), brinkman_penalty)
 
    ! ------------------------------------------------------------------------ !
    ! Set up output the brinkman fields if enabled
 
    if (output_enable) then
       fname = trim(output_path) // 'brinkman'
       select case (trim(output_format))
       case ('nek5000')
          suffix = '.fld'
       case ('adios2')
          suffix = '.bp'
       case ('hdf5')
          suffix = '.h5'
       case default
          suffix = '.' // trim(output_format)
       end select
       fname = trim(fname) // trim(suffix)
 
       select case (trim(output_precision))
       case ('sp', 'single')
          precision = sp
       case ('dp', 'double')
          precision = dp
       case default
          call neko_error('Unknown output precision')
       end select
 
       call output%init(fname, precision = precision)
 
       call output_fields%init(2)
       call output_fields%assign_to_field(1, this%indicator)
       call output_fields%assign_to_field(2, this%brinkman)
 
       call neko_log%message('Brinkman output')
       call neko_log%message('  1: Indicator')
       call neko_log%message('  2: Permeability')
       if (associated(this%unfiltered)) then
          call output_fields%append(this%unfiltered)
          call neko_log%message('  3: Unfiltered Indicator')
       end if
 
       call output%write(output_fields)
 
       call output%free()
       call output_fields%free()
    end if
 
    call neko_log%end_section()
  end subroutine brinkman_source_term_init_from_json
 
  subroutine brinkman_source_term_free(this)
    class(brinkman_source_term_t), intent(inout) :: this
 
    if (associated(this%indicator)) nullify(this%indicator)
    if (associated(this%unfiltered)) nullify(this%unfiltered)
    if (associated(this%brinkman)) nullify(this%brinkman)
 
    if (allocated(this%filter)) then
       call this%filter%free()
       deallocate(this%filter)
    end if
    call this%free_base()
 
  end subroutine brinkman_source_term_free
 
  subroutine brinkman_source_term_compute(this, time)
    class(brinkman_source_term_t), intent(inout) :: this
    type(time_state_t), intent(in) :: time
    type(field_t), pointer :: u, v, w, fu, fv, fw
    integer :: n
 
    n = this%fields%item_size(1)
 
    u => neko_registry%get_field('u')
    v => neko_registry%get_field('v')
    w => neko_registry%get_field('w')
 
    fu => this%fields%get(1)
    fv => this%fields%get(2)
    fw => this%fields%get(3)
 
    call field_subcol3(fu, u, this%brinkman, n)
    call field_subcol3(fv, v, this%brinkman, n)
    call field_subcol3(fw, w, this%brinkman, n)
 
    nullify(u, v, w, fu, fv, fw)
 
  end subroutine brinkman_source_term_compute
 
  ! ========================================================================== !
  ! Private methods
 
  subroutine add_boundary_mesh(this, json)
    class(brinkman_source_term_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
 
    ! Options
    character(len=:), allocatable :: mesh_file_name
    character(len=:), allocatable :: distance_transform
    character(len=:), allocatable :: filter_type
    character(len=:), allocatable :: mesh_transform
 
    ! Read the options for the boundary mesh
    type(file_t) :: mesh_file
    type(tri_mesh_t) :: boundary_mesh
    real(kind=rp) :: scalar_r
    real(kind=dp) :: scalar_d
    logical :: cache, cache_exist
    character(len=:), allocatable :: cache_filename
    type(file_t) :: cache_file
    type(fld_file_output_t) :: cache_output
    type(fld_file_data_t) :: cache_data
    type(global_interpolation_t) :: global_interp
    type(space_t) :: prev_Xh
    type(interpolator_t) :: space_interp
 
    ! Mesh transform options variables
    real(kind=dp), dimension(:), allocatable :: box_min, box_max
    logical :: keep_aspect_ratio
    real(kind=dp), dimension(3) :: scaling
    real(kind=dp), dimension(3) :: translation
    type(field_t), pointer :: temp_field
    type(aabb_t) :: mesh_box, target_box
    integer :: idx_p, temp_idx
    character(len=LOG_SIZE) :: log_msg
 
    ! ------------------------------------------------------------------------ !
    ! Read the options for the boundary mesh
 
    call json_get(json, 'name', mesh_file_name)
    call json_get_or_default(json, 'cache', cache, .false.)
 
    ! Settings on how to filter the design field
    call json_get(json, 'distance_transform.type', distance_transform)
 
    ! ------------------------------------------------------------------------ !
    ! Check if we can load from cache
    if (cache) then
       call json_get(json, 'cache_file', cache_filename)
 
       inquire(file=trim(cache_filename) // "0.nek5000", exist=cache_exist)
       write(log_msg, '(A)') "Checking for Brinkman source term cache."
       call neko_log%message(log_msg, neko_log_debug)
 
       if (cache_exist) then
          write(log_msg, '(A)') "Loading Brinkman source term from cache."
          call neko_log%message(log_msg, neko_log_debug)
 
          call neko_scratch_registry%request_field(temp_field, temp_idx, .true.)
 
          call cache_data%init()
          call cache_file%init(cache_filename // "0.fld")
          call cache_file%set_counter(0)
          call cache_file%read(cache_data)
          call cache_data%import_fields(p = temp_field)
 
          ! Update the global indicator field by max operator
          call field_pwmax2(this%indicator, temp_field)
 
          ! Clean up
          call neko_scratch_registry%relinquish(temp_idx)
          call cache_data%free()
          call cache_file%free()
          return
       end if
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Load the immersed boundary mesh
 
    call mesh_file%init(mesh_file_name)
    call mesh_file%read(boundary_mesh)
 
    if (boundary_mesh%nelv .eq. 0) then
       call neko_error('No elements in the boundary mesh')
    end if
 
    ! ------------------------------------------------------------------------ !
    ! Transform the mesh if specified.
 
    call json_get_or_default(json, 'mesh_transform.type', &
         mesh_transform, 'none')
 
    select case (mesh_transform)
    case ('none')
       ! Do nothing
    case ('bounding_box')
       call json_get_or_lookup(json, 'mesh_transform.box_min', box_min)
       call json_get_or_lookup(json, 'mesh_transform.box_max', box_max)
       call json_get_or_default(json, 'mesh_transform.keep_aspect_ratio', &
            keep_aspect_ratio, .true.)
 
       if (size(box_min) .ne. 3 .or. size(box_max) .ne. 3) then
          call neko_error('Case file: mesh_transform. &
          &box_min and box_max must be 3 element arrays of reals')
       end if
 
       call target_box%init(box_min, box_max)
 
       mesh_box = get_aabb(boundary_mesh)
 
       scaling = target_box%get_diagonal() / mesh_box%get_diagonal()
       if (keep_aspect_ratio) then
          scaling = minval(scaling)
       end if
 
       translation = - scaling * mesh_box%get_min() + target_box%get_min()
 
       do idx_p = 1, boundary_mesh%mpts
          boundary_mesh%points(idx_p)%x = &
               scaling * boundary_mesh%points(idx_p)%x + translation
       end do
 
       ! Report the transformation applied
       write(log_msg, '(A)') "The following transformation was applied:"
       call neko_log%message(log_msg)
       write(log_msg, '(A, 3F12.6)') "Scaling: ", scaling
       call neko_log%message(log_msg)
       write(log_msg, '(A, 3F12.6)') "Translation: ", translation
       call neko_log%message(log_msg)
 
    case default
       call neko_error('Unknown mesh transform')
    end select
 
    ! ------------------------------------------------------------------------ !
    ! Compute the permeability field
 
    ! Assign the signed distance field to all GLL points in the permeability
    ! field. Initially we just run a brute force loop over all GLL points and
    ! compute the signed distance function. This should be replaced with a
    ! more efficient method, such as a tree search.
 
    call neko_scratch_registry%request_field(temp_field, temp_idx, .true.)
 
    ! Select how to transform the distance field to a design field
    select case (distance_transform)
    case ('smooth_step')
       call json_get_or_lookup(json, 'distance_transform.value', scalar_d)
       scalar_r = real(scalar_d, kind=rp)
 
       call signed_distance_field(temp_field, boundary_mesh, scalar_d)
       call smooth_step_field(temp_field, scalar_r, 0.0_rp)
 
    case ('step')
 
       call json_get_or_lookup(json, 'distance_transform.value', scalar_d)
       scalar_r = real(scalar_d, kind=rp)
 
       call signed_distance_field(temp_field, boundary_mesh, scalar_d)
       call step_function_field(temp_field, scalar_r, 1.0_rp, 0.0_rp)
 
    case default
       call neko_error('Unknown distance transform')
    end select
 
    ! Write the field to cache
    if (cache) then
       write(log_msg, '(A)') "Writing Brinkman source term to cache."
       call neko_log%message(log_msg, neko_log_debug)
       call cache_output%init(dp, cache_filename, 1)
       call cache_output%fields%assign_to_field(1, temp_field)
       call cache_output%sample(0.0_rp)
    end if
 
    ! Update the global indicator field by max operator
    call field_pwmax2(this%indicator, temp_field)
 
    call neko_scratch_registry%relinquish(temp_idx)
 
  end subroutine add_boundary_mesh
 
  subroutine add_point_zone(this, json)
    class(brinkman_source_term_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
 
    ! Options
    character(len=:), allocatable :: zone_name
 
    type(field_t), pointer :: temp_field
    class(point_zone_t), pointer :: zone
    integer :: i, temp_idx
 
    ! ------------------------------------------------------------------------ !
    ! Read the options for the point zone
 
    call json_get(json, 'name', zone_name)
 
    ! Compute the indicator field
    call neko_scratch_registry%request_field(temp_field, temp_idx, .true.)
 
    zone => neko_point_zone_registry%get_point_zone(zone_name)
 
    if (neko_bcknd_device .eq. 1) then
       call device_cfill_mask(temp_field%x_d, 1.0_rp, temp_field%size(), &
            zone%mask%get_d(), zone%size)
    else
       call cfill_mask(temp_field%x, 1.0_rp, temp_field%size(), &
            zone%mask%get(), zone%size)
    end if
 
    ! Update the global indicator field by max operator
    call field_pwmax2(this%indicator, temp_field)
 
    call neko_scratch_registry%relinquish(temp_idx)
  end subroutine add_point_zone
 
  subroutine add_file(this, json)
    class(brinkman_source_term_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    type(file_t) :: file
    type(field_t), pointer :: temp_field
    character(len=:), allocatable :: file_name, field_name, tmp_str
    character(len=80) :: suffix
    integer :: file_idx, temp_idx
 
    call json_get(json, 'file_name', file_name)
    call json_get_or_default(json, 'field_name', field_name, &
         'brinkman_indicator')
    call json_get_or_default(json, 'file_index', file_idx, 0)
 
    call neko_scratch_registry%request_field(temp_field, temp_idx, .true.)
 
    call file%init(file_name)
    call file%set_counter(file_idx)
 
    call filename_suffix(file_name, suffix)
    select case (trim(suffix))
    case ('fld')
       block
         type(fld_file_data_t) :: fld_data
         type(field_list_t) :: fld_fields
         integer :: idx(1)
 
         call fld_data%init()
         call file%read(fld_data)
         select case (field_name(1:1))
         case ('p')
            call fld_data%import_fields(p = temp_field)
         case ('u')
            call fld_data%import_fields(u = temp_field)
         case ('v')
            call fld_data%import_fields(v = temp_field)
         case ('w')
            call fld_data%import_fields(w = temp_field)
         case ('t')
            call fld_data%import_fields(t = temp_field)
         case ('s')
 
            if (len_trim(field_name) .eq. 3) then
               read(field_name(2:3), '(I2)') idx(1)
            else if (len_trim(field_name) .eq. 2) then
               read(field_name(2:2), '(I1)') idx(1)
            else
               call neko_error('For fields with prefix s, the field name ' // &
                    'must be in the format sXX, where XX is the index of ' // &
                    'the field in the fld file')
            end if
 
            call fld_fields%init(1)
            call fld_fields%assign(1, temp_field)
 
            call fld_data%import_fields(s_target_list = fld_fields, &
                 s_index_list = idx)
 
            call fld_fields%free()
         case default
            call neko_error('Unknown field prefix in field name: ' // &
                 trim(field_name))
         end select
 
         call fld_data%free()
       end block
 
    case ('vtkhdf')
 
       ! VTKHDF will read the name of the field object.
       tmp_str = trim(temp_field%name)
       temp_field%name = trim(field_name)
 
       call file%read(temp_field)
 
       temp_field%name = trim(tmp_str)
 
    case default
       call neko_error("Brinkman cannot read file: " // trim(file_name))
    end select
 
    ! Update the global indicator field by max operator
    call field_pwmax2(this%indicator, temp_field)
 
    call neko_scratch_registry%relinquish(temp_idx)
    call file%free()
 
  end subroutine add_file
 
  subroutine add_field(this, json)
    class(brinkman_source_term_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    character(len=:), allocatable :: field_name
    type(field_t), pointer :: temp_field
 
    ! Read the field name
    call json_get(json, 'name', field_name)
    temp_field => neko_registry%get_field(field_name)
 
    ! Update the global indicator field by max operator
    call field_pwmax2(this%indicator, temp_field)
 
    if (associated(temp_field)) nullify(temp_field)
 
  end subroutine add_field
 
end module brinkman_source_term