case.f90

Go to the documentation of this file.

! Copyright (c) 2020-2023, 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 case
  use num_types, only : rp, sp, dp
  use fluid_pnpn, only : fluid_pnpn_t
  use fluid_scheme, only : fluid_scheme_t, fluid_scheme_factory
  use fluid_output, only : fluid_output_t
  use chkp_output, only : chkp_output_t
  use mean_sqr_flow_output, only : mean_sqr_flow_output_t
  use mean_flow_output, only : mean_flow_output_t
  use fluid_stats_output, only : fluid_stats_output_t
  use mesh_field, only : mesh_fld_t, mesh_field_init, mesh_field_free
  use parmetis, only : parmetis_partmeshkway
  use redist, only : redist_mesh
  use sampler, only : sampler_t
  use flow_ic, only : set_flow_ic
  use scalar_ic, only : set_scalar_ic
  use stats, only : stats_t
  use file, only : file_t
  use utils, only : neko_error
  use mesh, only : mesh_t
  use comm
  use time_scheme_controller, only : time_scheme_controller_t
  use logger, only : neko_log, neko_log_quiet, log_size
  use jobctrl, only : jobctrl_set_time_limit
  use user_intf, only : user_t
  use scalar_pnpn, only : scalar_pnpn_t
  use json_module, only : json_file, json_core, json_value
  use json_utils, only : json_get, json_get_or_default
  use scratch_registry, only : scratch_registry_t, neko_scratch_registry
  use point_zone_registry, only: neko_point_zone_registry
  use material_properties, only : material_properties_t
  implicit none
  private
 
  type, public :: case_t
     type(mesh_t) :: msh
     type(json_file) :: params
     type(time_scheme_controller_t) :: ext_bdf
     real(kind=rp), dimension(10) :: tlag
     real(kind=rp), dimension(10) :: dtlag
     real(kind=rp) :: dt
     real(kind=rp) :: end_time
     type(sampler_t) :: s
     type(fluid_output_t) :: f_out
     type(fluid_stats_output_t) :: f_stats_output
     type(chkp_output_t) :: f_chkp
     type(mean_flow_output_t) :: f_mf
     type(mean_sqr_flow_output_t) :: f_msqrf
     type(stats_t) :: q
     type(user_t) :: usr
     class(fluid_scheme_t), allocatable :: fluid
     type(scalar_pnpn_t), allocatable :: scalar
     type(material_properties_t) :: material_properties
  end type case_t
 
  interface case_init
     module procedure case_init_from_file, case_init_from_json
  end interface case_init
 
  public :: case_init, case_free
 
contains
 
  subroutine case_init_from_file(C, case_file)
    type(case_t), target, intent(inout) :: C
    character(len=*), intent(in) :: case_file
    integer :: ierr, integer_val
    character(len=:), allocatable :: json_buffer
 
    call neko_log%section('Case')
    call neko_log%message('Reading case file ' // trim(case_file), &
                          neko_log_quiet)
 
    if (pe_rank .eq. 0) then
       call c%params%load_file(filename = trim(case_file))
       call c%params%print_to_string(json_buffer)
       integer_val = len(json_buffer)
    end if
 
    call mpi_bcast(integer_val, 1, mpi_integer, 0, neko_comm, ierr)
    if (pe_rank .ne. 0) allocate(character(len=integer_val)::json_buffer)
    call mpi_bcast(json_buffer, integer_val, mpi_character, 0, neko_comm, ierr)
    call c%params%load_from_string(json_buffer)
 
    deallocate(json_buffer)
 
    call case_init_common(c)
 
  end subroutine case_init_from_file
 
  subroutine case_init_from_json(C, case_json)
    type(case_t), target, intent(inout) :: C
    type(json_file), intent(in) :: case_json
 
    call neko_log%section('Case')
    call neko_log%message('Creating case from JSON object', neko_log_quiet)
 
    c%params = case_json
 
    call case_init_common(c)
 
  end subroutine case_init_from_json
 
  subroutine case_init_common(C)
    type(case_t), target, intent(inout) :: C
    character(len=:), allocatable :: output_directory
    integer :: lx = 0
    logical :: scalar = .false.
    type(file_t) :: msh_file, bdry_file, part_file
    type(mesh_fld_t) :: msh_part, parts
    logical :: found, logical_val
    integer :: integer_val
    real(kind=rp) :: real_val
    character(len=:), allocatable :: string_val
    real(kind=rp) :: stats_start_time, stats_output_val
    integer :: stats_sampling_interval
    integer :: output_dir_len
    integer :: precision
 
    !
    ! Load mesh
    !
    call json_get_or_default(c%params, 'case.mesh_file', string_val, 'no mesh')
    if (trim(string_val) .eq. 'no mesh') then
       call neko_error('The mesh_file keyword could not be found in the .' // &
                       'case file. Often caused by incorrectly formatted json.')
    end if
    msh_file = file_t(string_val)
 
    call msh_file%read(c%msh)
 
    !
    ! Load Balancing
    !
    call json_get_or_default(c%params, 'case.load_balancing', logical_val,&
                             .false.)
 
    if (pe_size .gt. 1 .and. logical_val) then
       call neko_log%section('Load Balancing')
       call parmetis_partmeshkway(c%msh, parts)
       call redist_mesh(c%msh, parts)
       call neko_log%end_section()
    end if
 
    !
    ! Time step
    !
    call c%params%get('case.variable_timestep', logical_val, found)
    if (.not. logical_val) then
       call json_get(c%params, 'case.timestep', c%dt)
    else
       ! randomly set an initial dt to get cfl when dt is variable
       c%dt = 1.0_rp
    end if
 
    !
    ! End time
    !
    call json_get(c%params, 'case.end_time', c%end_time)
 
    !
    ! Initialize point_zones registry
    !
    call neko_point_zone_registry%init(c%params, c%msh)
 
    !
    ! Setup user defined functions
    !
    call c%usr%init()
    call c%usr%user_mesh_setup(c%msh)
 
    !
    ! Set order of timestepper
    !
    call json_get(c%params, 'case.numerics.time_order', integer_val)
    call c%ext_bdf%init(integer_val)
 
 
    !
    ! Material properties
    !
    call c%material_properties%init(c%params, c%usr)
 
    !
    ! Setup fluid scheme
    !
    call json_get(c%params, 'case.fluid.scheme', string_val)
    call fluid_scheme_factory(c%fluid, trim(string_val))
 
    call json_get(c%params, 'case.numerics.polynomial_order', lx)
    lx = lx + 1 ! add 1 to get number of gll points
    c%fluid%chkp%tlag => c%tlag
    c%fluid%chkp%dtlag => c%dtlag
    call c%fluid%init(c%msh, lx, c%params, c%usr, c%material_properties, &
                      c%ext_bdf)
    select type (f => c%fluid)
    type is (fluid_pnpn_t)
       f%chkp%abx1 => f%abx1
       f%chkp%abx2 => f%abx2
       f%chkp%aby1 => f%aby1
       f%chkp%aby2 => f%aby2
       f%chkp%abz1 => f%abz1
       f%chkp%abz2 => f%abz2
    end select
 
 
    !
    ! Setup scratch registry
    !
    neko_scratch_registry = scratch_registry_t(c%fluid%dm_Xh, 10, 10)
 
    !
    ! Setup scalar scheme
    !
    ! @todo no scalar factory for now, probably not needed
    if (c%params%valid_path('case.scalar')) then
       call json_get_or_default(c%params, 'case.scalar.enabled', scalar,&
                                .true.)
    end if
 
    if (scalar) then
       allocate(c%scalar)
       c%scalar%chkp%tlag => c%tlag
       c%scalar%chkp%dtlag => c%dtlag
       call c%scalar%init(c%msh, c%fluid%c_Xh, c%fluid%gs_Xh, c%params, c%usr,&
                          c%material_properties, c%fluid%ulag, c%fluid%vlag, &
                          c%fluid%wlag, c%ext_bdf)
       call c%fluid%chkp%add_scalar(c%scalar%s)
       c%fluid%chkp%abs1 => c%scalar%abx1
       c%fluid%chkp%abs2 => c%scalar%abx2
       c%fluid%chkp%slag => c%scalar%slag
    end if
 
    !
    ! Setup user defined conditions
    !
    if (c%params%valid_path('case.fluid.inflow_condition')) then
       call json_get(c%params, 'case.fluid.inflow_condition.type',&
                     string_val)
       if (trim(string_val) .eq. 'user') then
          call c%fluid%set_usr_inflow(c%usr%fluid_user_if)
       end if
    end if
 
    ! Setup user boundary conditions for the scalar.
    if (scalar) then
       call c%scalar%set_user_bc(c%usr%scalar_user_bc)
    end if
 
    !
    ! Setup initial conditions
    !
    call json_get(c%params, 'case.fluid.initial_condition.type',&
                  string_val)
    if (trim(string_val) .ne. 'user') then
       call set_flow_ic(c%fluid%u, c%fluid%v, c%fluid%w, c%fluid%p, &
            c%fluid%c_Xh, c%fluid%gs_Xh, string_val, c%params)
    else
       call set_flow_ic(c%fluid%u, c%fluid%v, c%fluid%w, c%fluid%p, &
            c%fluid%c_Xh, c%fluid%gs_Xh, c%usr%fluid_user_ic, c%params)
    end if
 
    if (scalar) then
       call json_get(c%params, 'case.scalar.initial_condition.type', string_val)
       if (trim(string_val) .ne. 'user') then
          call set_scalar_ic(c%scalar%s, &
            c%scalar%c_Xh, c%scalar%gs_Xh, string_val, c%params)
       else
          call set_scalar_ic(c%scalar%s, &
            c%scalar%c_Xh, c%scalar%gs_Xh, c%usr%scalar_user_ic, c%params)
       end if
    end if
 
    ! Add initial conditions to BDF scheme (if present)
    select type (f => c%fluid)
    type is (fluid_pnpn_t)
       call f%ulag%set(f%u)
       call f%vlag%set(f%v)
       call f%wlag%set(f%w)
    end select
 
    !
    ! Validate that the case is properly setup for time-stepping
    !
    call c%fluid%validate
 
    if (scalar) then
       call c%scalar%slag%set(c%scalar%s)
       call c%scalar%validate
    end if
 
    !
    ! Get and process output directory
    !
    call json_get_or_default(c%params, 'case.output_directory',&
                             output_directory, '')
 
    output_dir_len = len(trim(output_directory))
    if (output_dir_len .gt. 0) then
       if (output_directory(output_dir_len:output_dir_len) .ne. "/") then
          output_directory = trim(output_directory)//"/"
          if (pe_rank .eq. 0) then
             call execute_command_line('mkdir -p '//output_directory)
          end if
       end if
    end if
 
    !
    ! Save boundary markings for fluid (if requested)
    !
    call json_get_or_default(c%params, 'case.output_boundary',&
                             logical_val, .false.)
    if (logical_val) then
       bdry_file = file_t(trim(output_directory)//'bdry.fld')
       call bdry_file%write(c%fluid%bdry)
    end if
 
    !
    ! Save mesh partitions (if requested)
    !
    call json_get_or_default(c%params, 'case.output_partitions',&
                             logical_val, .false.)
    if (logical_val) then
       call mesh_field_init(msh_part, c%msh, 'MPI_Rank')
       msh_part%data = pe_rank
       part_file = file_t(trim(output_directory)//'partitions.vtk')
       call part_file%write(msh_part)
       call mesh_field_free(msh_part)
    end if
 
    !
    ! Setup output precision of the field files
    !
    call json_get_or_default(c%params, 'case.output_precision', string_val,&
         'single')
 
    if (trim(string_val) .eq. 'double') then
       precision = dp
    else
       precision = sp
    end if
 
    !
    ! Setup sampler
    !
    call c%s%init(c%end_time)
    if (scalar) then
       c%f_out = fluid_output_t(precision, c%fluid, c%scalar, &
            path = trim(output_directory))
    else
       c%f_out = fluid_output_t(precision, c%fluid, &
            path = trim(output_directory))
    end if
 
    call json_get_or_default(c%params, 'case.fluid.output_control',&
                             string_val, 'org')
 
    if (trim(string_val) .eq. 'org') then
       ! yes, it should be real_val below for type compatibility
       call json_get(c%params, 'case.nsamples', real_val)
       call c%s%add(c%f_out, real_val, 'nsamples')
    else if (trim(string_val) .eq. 'never') then
       ! Fix a dummy 0.0 output_value
       call json_get_or_default(c%params, 'case.fluid.output_value', real_val, &
                                0.0_rp)
       call c%s%add(c%f_out, 0.0_rp, string_val)
    else
       call json_get(c%params, 'case.fluid.output_value', real_val)
       call c%s%add(c%f_out, real_val, string_val)
    end if
 
    !
    ! Save checkpoints (if nothing specified, default to saving at end of sim)
    !
    call json_get_or_default(c%params, 'case.output_checkpoints',&
                             logical_val, .true.)
    if (logical_val) then
       call json_get_or_default(c%params, 'case.checkpoint_format', &
            string_val, "chkp")
       c%f_chkp = chkp_output_t(c%fluid%chkp, path = output_directory, &
            fmt = trim(string_val))
       call json_get_or_default(c%params, 'case.checkpoint_control', &
            string_val, "simulationtime")
       call json_get_or_default(c%params, 'case.checkpoint_value', real_val,&
            1e10_rp)
       call c%s%add(c%f_chkp, real_val, string_val)
    end if
 
    !
    ! Setup statistics
    !
 
    ! Always init, so that we can call eval in simulation.f90 with no if.
    ! Note, don't use json_get_or_default here, because that will break the
    ! valid_path if statement below (the path will become valid always).
    call c%params%get('case.statistics.start_time', stats_start_time,&
                           found)
    if (.not. found) stats_start_time = 0.0_rp
 
    call c%params%get('case.statistics.sampling_interval', &
                           stats_sampling_interval, found)
    if (.not. found) stats_sampling_interval = 10
 
    call c%q%init(stats_start_time, stats_sampling_interval)
 
    found = c%params%valid_path('case.statistics')
    if (found) then
       call json_get_or_default(c%params, 'case.statistics.enabled',&
                                logical_val, .true.)
       if (logical_val) then
          call c%q%add(c%fluid%mean%u)
          call c%q%add(c%fluid%mean%v)
          call c%q%add(c%fluid%mean%w)
          call c%q%add(c%fluid%mean%p)
 
          c%f_mf = mean_flow_output_t(c%fluid%mean, stats_start_time, &
                                      path = output_directory)
 
          call json_get(c%params, 'case.statistics.output_control', &
                        string_val)
          call json_get(c%params, 'case.statistics.output_value', &
                        stats_output_val)
 
          call c%s%add(c%f_mf, stats_output_val, string_val)
          call c%q%add(c%fluid%stats)
 
          c%f_stats_output = fluid_stats_output_t(c%fluid%stats, &
            stats_start_time, path = output_directory)
          call c%s%add(c%f_stats_output, stats_output_val, string_val)
       end if
    end if
 
    !
    ! Setup joblimit
    !
    if (c%params%valid_path('case.job_timelimit')) then
       call json_get(c%params, 'case.job_timelimit', string_val)
       call jobctrl_set_time_limit(string_val)
    end if
 
    call neko_log%end_section()
 
  end subroutine case_init_common
 
  subroutine case_free(C)
    type(case_t), intent(inout) :: c
 
    if (allocated(c%fluid)) then
       call c%fluid%free()
       deallocate(c%fluid)
    end if
 
    if (allocated(c%scalar)) then
       call c%scalar%free()
       deallocate(c%scalar)
    end if
 
    call c%msh%free()
 
    call c%s%free()
 
    call c%q%free()
 
  end subroutine case_free
 
end module case