!> @file initial_conditions.f90 !> @brief Initial-condition setup for all supported problem types. !! !! Provides a single public entry point `apply_initial_condition` that fills !! the conserved-state array `state%ub` (interior + halo) for the chosen !! problem. For boundary conditions that use early-return in apply_bcs !! (dirichlet / inflow / supersonic_inlet), this routine also pre-fills the !! halo cells with the prescribed boundary state. !! !! Supported problems: !! 'sod' — Sod shock tube !! 'shu_osher' — Shu-Osher shock/entropy-wave interaction !! 'smooth_wave' — Smooth advecting density wave (periodic) !! 'linear_advection' — Sinusoidal density on uniform background flow (periodic) !! 'woodward_colella' — Woodward-Colella interacting blast waves !! 'lax' — Lax shock tube !! 'from_file' — Read primitives (x, rho, u, p) from an external data file !! 'udf' — Call a user-defined Fortran subroutine compiled at runtime !! 'acoustic_pulse' — Gaussian acoustic pulse (quiescent background; NRBC benchmark) module initial_conditions use precision, only: wp use config, only: config_t use solver_state, only: solver_state_t, neq use boundary_conditions, only: apply_bcs use logger, only: log_info use option_registry, only: problem_sod, problem_shu_osher, problem_smooth_wave, & problem_linear_advection, problem_woodward_colella, & problem_lax, problem_from_file, problem_udf, & problem_acoustic_pulse, method_fvm use iso_fortran_env, only: iostat_end use iso_c_binding, only: c_char, c_double, c_int, c_null_char implicit none private public :: apply_initial_condition, quote_shell_arg integer(c_int), parameter :: host_platform_windows = 1_c_int integer(c_int), parameter :: host_platform_macos = 2_c_int integer(c_int), parameter :: host_platform_linux = 3_c_int interface function cfd_solver_host_platform() bind(C, name="cfd_solver_host_platform") import :: c_int integer(c_int) :: cfd_solver_host_platform end function cfd_solver_host_platform function cfd_solver_call_udf(libpath, n, x, rho, rho_u, e_arr) bind(C, name="cfd_solver_call_udf") import :: c_char, c_double, c_int integer(c_int) :: cfd_solver_call_udf character(kind=c_char), intent(in) :: libpath(*) integer(c_int), value, intent(in) :: n real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: rho(n) real(c_double), intent(out) :: rho_u(n) real(c_double), intent(out) :: e_arr(n) end function cfd_solver_call_udf function cfd_solver_get_pid() bind(C, name="cfd_solver_get_pid") import :: c_int integer(c_int) :: cfd_solver_get_pid end function cfd_solver_get_pid end interface contains !> Return the host platform detected by the C portability shim. function get_host_platform() result(platform) integer(c_int) :: platform platform = cfd_solver_host_platform() end function get_host_platform !> Return the preferred runtime-UDF output directory on Windows. function get_windows_temp_dir() result(temp_dir) character(len=512) :: temp_dir integer :: status, value_len temp_dir = '' call get_environment_variable('TEMP', temp_dir, length=value_len, status=status) if (status == 0 .and. value_len > 0) then temp_dir = temp_dir(:value_len) return end if temp_dir = '' call get_environment_variable('TMP', temp_dir, length=value_len, status=status) if (status == 0 .and. value_len > 0) then temp_dir = temp_dir(:value_len) return end if temp_dir = '.' end function get_windows_temp_dir !> Join a directory and leaf name with a platform-appropriate separator. function join_path(dirpath, leafname, separator) result(full_path) character(len=*), intent(in) :: dirpath, leafname character(len=1), intent(in) :: separator character(len=len_trim(dirpath) + len_trim(leafname) + 2) :: full_path integer :: dir_len character(len=1) :: last_char dir_len = len_trim(dirpath) if (dir_len == 0) then full_path = trim(leafname) return end if last_char = dirpath(dir_len:dir_len) if (trim(dirpath) == '.') then full_path = '.'//separator//trim(leafname) else if (last_char == '/' .or. last_char == '\') then full_path = trim(dirpath)//trim(leafname) else full_path = trim(dirpath)//separator//trim(leafname) end if end function join_path !> Wrap a filesystem path for the shell used by execute_command_line. !! !! Wraps the path in double quotes, blocking word-splitting on whitespace. !! Paths containing `"` are rejected because the double-quote cannot be !! escaped inside double-quoted shell arguments on all POSIX shells. !! Paths containing `$` or a backtick (`` ` ``) are also rejected because !! both characters expand inside double quotes (covering `${`, `$(`, and !! command-substitution forms), making the path unsafe to quote-and-invoke. function quote_shell_arg(raw_path) result(quoted_path) character(len=*), intent(in) :: raw_path character(len=len_trim(raw_path) + 2) :: quoted_path if (index(raw_path, '"') /= 0 .or. index(raw_path, '$') /= 0 .or. & index(raw_path, '`') /= 0) & error stop 'initial_conditions: path contains a shell metacharacter ("/$`) and cannot be quoted safely' quoted_path = '"'//trim(raw_path)//'"' end function quote_shell_arg !> Fill `state%ub` (interior + halo) with the initial condition selected by !! `cfg%problem_type`. Boundary halo cells are also pre-populated for !! Dirichlet-style BCs (dirichlet/inflow/supersonic_inlet) which use early !! return in apply_bcs. !! !! @param[inout] state Solver state (must have ub allocated, schemes bound) !! @param[in] cfg Runtime configuration subroutine apply_initial_condition(state, cfg) type(solver_state_t), intent(inout), target :: state ! target required for x_coord pointer type(config_t), intent(in) :: cfg real(wp), parameter :: pi = 4.0_wp * atan(1.0_wp) integer :: ipt real(wp) :: x, rho_so real(wp) :: q_l(neq), q_r(neq) logical :: is_fvm ! Coordinate pointer: points to x_cell (FVM) or x_node (FDM) so that ! every per-point loop body reads x = x_coord(ipt) without duplicating the ! select-case. smooth_wave + FVM uses a separate quadrature path below. real(wp), pointer :: x_coord(:) is_fvm = (trim(state % blocks(1) % method) == method_fvm) if (is_fvm) then x_coord => state % mesh % x_cell else x_coord => state % mesh % x_node end if select case (trim(cfg % problem_type)) case (problem_sod, problem_lax) ! Sod (1978) / Lax (1954) shock tubes — identical setup; the two ! problems differ only in their default left/right primitives, which ! are supplied via the namelist. ! References: Sod, J. Comput. Phys. 27:1-31, 1978; ! Lax, Commun. Pure Appl. Math. 7:159-193, 1954. call left_right_states_from_namelist(cfg, state % cfg % gam, q_l, q_r) do ipt = 1, state % n_pt x = x_coord(ipt) if (x < cfg % x_diaphragm) then state % ub(:, ipt) = q_l else if (x > cfg % x_diaphragm) then state % ub(:, ipt) = q_r else state % ub(:, ipt) = 0.5_wp * (q_l + q_r) end if end do call seed_dirichlet_halos(state, q_l, q_r) case (problem_shu_osher) ! Shu-Osher shock/entropy-wave interaction. ! See Shu & Osher (1989), J. Comput. Phys. 83:32-78. ! Left (x < -4): rho=3.857143, u=2.629369, p=10.333333 (post-shock) ! Right (x >= -4): rho = 1 + 0.2*sin(5*x), u = 0, p = 1 ! Left ghost: constant post-shock state (supersonic inflow) q_l(1) = 3.857143_wp q_l(2) = 3.857143_wp * 2.629369_wp q_l(3) = 10.333333_wp / (state % cfg % gam - 1.0_wp) & + 0.5_wp * 3.857143_wp * 2.629369_wp**2 ! Right ghost: undisturbed far-field (rho=1, u=0, p=1). q_r(1) = 1.0_wp q_r(2) = 0.0_wp q_r(3) = 1.0_wp / (state % cfg % gam - 1.0_wp) do ipt = 1, state % n_pt x = x_coord(ipt) if (x < -4.0_wp) then state % ub(:, ipt) = q_l else rho_so = 1.0_wp + 0.2_wp * sin(5.0_wp * x) state % ub(1, ipt) = rho_so state % ub(2, ipt) = 0.0_wp state % ub(3, ipt) = 1.0_wp / (state % cfg % gam - 1.0_wp) end if end do call seed_dirichlet_halos(state, q_l, q_r) case (problem_smooth_wave) ! Smooth advecting density wave on periodic domain [0, 1]. ! FVM path: cell-averaged IC using 4-point Gauss-Legendre quadrature over ! each cell [x_cell(i)-dx/2, x_cell(i)+dx/2]. Primitives are averaged ! (rho, u, p) then converted to conserved. For this problem u=1 and p=1 ! are constant, so only rho needs quadrature. The 4-point GL rule is exact ! for polynomials of degree ≤ 7, giving essentially machine-precision ! cell averages of the sine wave on any reasonable grid. ! FDM path: unchanged nodal sampling at x_node(ipt). if (is_fvm) then block ! 4-point Gauss-Legendre nodes and weights on [-1, 1]. real(wp), parameter :: tgl(4) = [ & -0.8611363115940526_wp, -0.3399810435848563_wp, & 0.3399810435848563_wp, 0.8611363115940526_wp] real(wp), parameter :: wgl(4) = [ & 0.3478548451374538_wp, 0.6521451548625461_wp, & 0.6521451548625461_wp, 0.3478548451374538_wp] real(wp) :: xc, dxh, rho_c, xi_k integer :: k do ipt = 1, state % n_pt xc = state % mesh % x_cell(ipt) dxh = 0.5_wp * state % mesh % dx_cell(ipt) ! rho_avg = 0.5 * sum_k(wgl(k) * rho(xc + dxh*tgl(k))) ! (factor 0.5 maps the [-1,1] quadrature weight to [a,b]/dx) rho_c = 0.0_wp do k = 1, 4 xi_k = xc + dxh * tgl(k) rho_c = rho_c + wgl(k) * (1.0_wp + 0.2_wp * sin(2.0_wp * pi * xi_k)) end do rho_c = 0.5_wp * rho_c state % ub(1, ipt) = rho_c state % ub(2, ipt) = rho_c ! rho*u (u = 1) state % ub(3, ipt) = 1.0_wp / (state % cfg % gam - 1.0_wp) & + 0.5_wp * rho_c ! E = p/(g-1) + 0.5*rho*u^2 end do end block else do ipt = 1, state % n_pt x = x_coord(ipt) state % ub(1, ipt) = 1.0_wp + 0.2_wp * sin(2.0_wp * pi * x) state % ub(2, ipt) = state % ub(1, ipt) ! rho*u (u = 1) state % ub(3, ipt) = 1.0_wp / (state % cfg % gam - 1.0_wp) & + 0.5_wp * state % ub(1, ipt) ! E = p/(g-1) + 0.5*rho*u^2 end do end if case (problem_linear_advection) ! Sinusoidal density wave on a uniform background flow, periodic domain. do ipt = 1, state % n_pt x = x_coord(ipt) state % ub(1, ipt) = 1.0_wp + 0.5_wp & * sin(2.0_wp * pi * (x - state % cfg % x_left) & / (state % cfg % x_right - state % cfg % x_left)) state % ub(2, ipt) = state % ub(1, ipt) * cfg % u_left ! rho * u_adv state % ub(3, ipt) = 1.0_wp / (state % cfg % gam - 1.0_wp) & + 0.5_wp * state % ub(1, ipt) * cfg % u_left**2 end do case (problem_woodward_colella) ! Woodward-Colella interacting blast waves. ! See Woodward & Colella (1984), J. Comput. Phys., 54:115-173. ! Three-region IC: rho=1, u=0 everywhere; pressures differ. ! BCs: reflecting walls. do ipt = 1, state % n_pt x = x_coord(ipt) state % ub(1, ipt) = 1.0_wp state % ub(2, ipt) = 0.0_wp if (x < 0.1_wp) then state % ub(3, ipt) = 1000.0_wp / (state % cfg % gam - 1.0_wp) else if (x < 0.9_wp) then state % ub(3, ipt) = 0.01_wp / (state % cfg % gam - 1.0_wp) else state % ub(3, ipt) = 100.0_wp / (state % cfg % gam - 1.0_wp) end if end do case (problem_from_file) call ic_from_file(state, cfg) case (problem_udf) call ic_udf(state, cfg) case (problem_acoustic_pulse) ! Gaussian acoustic pulse on a quiescent uniform background. ! Canonical benchmark for non-reflecting boundary conditions. block real(wp) :: x0, dp, alpha, c0sq, rho0, p0 rho0 = 1.0_wp p0 = 1.0_wp c0sq = state % cfg % gam * p0 / rho0 ! c0^2 dp = 1.0e-3_wp ! pulse amplitude [Pa] alpha = 200.0_wp ! half-width: sigma = sqrt(ln2/alpha) x0 = 0.5_wp * (state % cfg % x_left + state % cfg % x_right) ! domain centre do ipt = 1, state % n_pt x = x_coord(ipt) associate (q => state % ub(:, ipt)) q(1) = rho0 + (dp / c0sq) * exp(-alpha * (x - x0)**2) ! rho q(2) = 0.0_wp ! rho*u = 0 q(3) = (p0 + dp * exp(-alpha * (x - x0)**2)) & ! E & / (state % cfg % gam - 1.0_wp) end associate end do ! Seed halos with uniform background (NRBC re-derives the ghost each ! call, but the very first halo_exchange/apply_bcs pass on an ! uninitialised halo would feed garbage into the FVS precompute). q_l = [rho0, 0.0_wp, p0 / (state % cfg % gam - 1.0_wp)] q_r = q_l call seed_dirichlet_halos(state, q_l, q_r) end block case default error stop 'initial_conditions: unknown problem_type' end select ! Final BC pass so that the very first compute_resid sees a consistent ! halo (apply_bcs is a no-op on inner ranks at this stage because ! state%decomp is fully populated by solver_runtime). call apply_bcs(state) end subroutine apply_initial_condition ! --------------------------------------------------------------------------- ! Private helpers ! --------------------------------------------------------------------------- !> Convert namelist primitives (rho_left, u_left, p_left / rho_right, ...) to !! conserved state vectors. pure subroutine left_right_states_from_namelist(cfg, gam, q_l, q_r) type(config_t), intent(in) :: cfg real(wp), intent(in) :: gam real(wp), intent(out) :: q_l(neq), q_r(neq) q_l(1) = cfg % rho_left q_l(2) = cfg % rho_left * cfg % u_left q_l(3) = cfg % p_left / (gam - 1.0_wp) & + 0.5_wp * cfg % rho_left * cfg % u_left**2 q_r(1) = cfg % rho_right q_r(2) = cfg % rho_right * cfg % u_right q_r(3) = cfg % p_right / (gam - 1.0_wp) & + 0.5_wp * cfg % rho_right * cfg % u_right**2 end subroutine left_right_states_from_namelist !> Fill the boundary halo cells with the prescribed Dirichlet states. !! !! Only edge ranks have boundary halos that BC code-paths leave alone — those !! ranks need a pre-filled halo so apply_bcs's early-return (for dirichlet / !! inflow / supersonic_inlet) reads valid data. Inner ranks: this is a !! harmless overwrite that halo_exchange will replace on the first residual !! evaluation. subroutine seed_dirichlet_halos(state, q_l, q_r) type(solver_state_t), intent(inout) :: state real(wp), intent(in) :: q_l(neq), q_r(neq) integer :: k, h, n_local h = state % decomp % halo_width n_local = state % decomp % n_local do k = 1, h state % ub(:, 1 - k) = q_l state % ub(:, n_local + k) = q_r end do end subroutine seed_dirichlet_halos !> Load the initial condition from a data file containing columns: x rho u p. !! !! If the file grid matches the solver grid (same number of points and x-coordinates !! within a small tolerance), the data is loaded directly. Otherwise, each primitive !! is linearly interpolated onto the solver grid. Interpolation must be enabled via !! `cfg%ic_interp`; a grid mismatch with `ic_interp = .false.` causes an error stop. !! Extrapolation (solver domain wider than file domain) is always an error stop. !! !! After filling `state%ub`, halo states are seeded from the outermost cells. !! !! @param[inout] state Solver state (must have ub allocated, schemes bound) !! @param[in] cfg Runtime configuration (uses cfg%ic_file, cfg%ic_interp, cfg%gam) subroutine ic_from_file(state, cfg) type(solver_state_t), intent(inout) :: state type(config_t), intent(in) :: cfg integer :: u, info, n_file_pt, i, ipt, lo real(wp) :: x_f_dummy, rho_dummy, u_dummy, p_dummy real(wp), allocatable :: x_f(:), rho_f(:), u_f(:), p_f(:) real(wp) :: x_i, t, rho_i, u_i, p_i real(wp), parameter :: tol = 1.0e-10_wp logical :: exact_match character(len=512) :: msg ! --- 1. Open file --- open (newunit=u, file=trim(cfg % ic_file), status='old', action='read', iostat=info) if (info /= 0) & error stop 'initial_conditions: cannot open ic_file "'//trim(cfg % ic_file)//'"' ! --- 2. Count lines --- n_file_pt = 0 do read (u, *, iostat=info) x_f_dummy, rho_dummy, u_dummy, p_dummy if (info == iostat_end) exit if (info /= 0) & error stop 'initial_conditions: read error while counting lines in ic_file' n_file_pt = n_file_pt + 1 end do if (n_file_pt < 2) & error stop 'initial_conditions: ic_file must contain at least 2 data lines' ! --- 3. Read data --- allocate (x_f(n_file_pt), rho_f(n_file_pt), u_f(n_file_pt), p_f(n_file_pt), stat=info) if (info /= 0) error stop 'initial_conditions: allocate failed for ic_from_file arrays' rewind (u) do i = 1, n_file_pt read (u, *, iostat=info) x_f(i), rho_f(i), u_f(i), p_f(i) if (info /= 0) & error stop 'initial_conditions: read error in ic_file data' end do close (u, iostat=info) if (info /= 0) error stop 'initial_conditions: file close failed for ic_file' ! Verify monotonicity of file x-grid. do i = 2, n_file_pt if (x_f(i) <= x_f(i - 1)) & error stop 'initial_conditions: ic_file x-coordinates must be strictly increasing' end do ! Check whether solver domain fits within file domain (no extrapolation). if (state % cfg % x_left < x_f(1) - tol .or. & state % cfg % x_right > x_f(n_file_pt) + tol) then error stop 'initial_conditions: solver domain extends beyond ic_file domain (extrapolation not allowed)' end if ! --- 4. Detect exact match --- exact_match = (n_file_pt == state % n_pt) if (exact_match) then do i = 1, n_file_pt x_i = state % mesh % x_node(i) if (abs(x_f(i) - x_i) > tol) then exact_match = .false. exit end if end do end if if (.not. exact_match) then if (.not. cfg % ic_interp) & error stop 'initial_conditions: ic_file grid does not match solver grid and ic_interp = .false.' write (msg, '(A,I0,A,I0,A)') & 'IC file: ', n_file_pt, ' pts interpolated onto ', state % n_pt, '-pt solver grid' call log_info(trim(msg)) end if ! --- 5. Fill state%ub --- lo = 1 ! lower bracket index for sequential search do ipt = 1, state % n_pt x_i = state % mesh % x_node(ipt) if (exact_match) then rho_i = rho_f(ipt) u_i = u_f(ipt) p_i = p_f(ipt) else ! Advance lo until x_f(lo+1) >= x_i (sequential scan, O(n) total). do while (lo < n_file_pt - 1 .and. x_f(lo + 1) < x_i - tol) lo = lo + 1 end do ! Linear interpolation weight. t = (x_i - x_f(lo)) / (x_f(lo + 1) - x_f(lo)) t = max(0.0_wp, min(1.0_wp, t)) ! clamp for floating-point edge cases at boundary rho_i = rho_f(lo) + t * (rho_f(lo + 1) - rho_f(lo)) u_i = u_f(lo) + t * (u_f(lo + 1) - u_f(lo)) p_i = p_f(lo) + t * (p_f(lo + 1) - p_f(lo)) end if ! Convert primitives to conserved. state % ub(1, ipt) = rho_i state % ub(2, ipt) = rho_i * u_i state % ub(3, ipt) = p_i / (state % cfg % gam - 1.0_wp) + 0.5_wp * rho_i * u_i**2 end do deallocate (x_f, rho_f, u_f, p_f, stat=info) if (info /= 0) error stop 'initial_conditions: deallocate failed for ic_from_file arrays' ! --- 6. Seed boundary halos from the outermost loaded cells. --- call seed_dirichlet_halos(state, state % ub(:, 1), state % ub(:, state % n_pt)) end subroutine ic_from_file !> Compile a user-defined Fortran source file to a shared library at runtime !! and call its `ic_udf` subroutine to fill the initial state. !! !! The source file must contain exactly one subroutine with the interface: !! !! subroutine ic_udf(n, x, rho, rho_u, E) bind(C, name="ic_udf") !! use iso_c_binding, only: c_int, c_double !! integer(c_int), value, intent(in) :: n !! real(c_double), intent(in) :: x(n) ! grid-point coordinates [m] !! real(c_double), intent(out) :: rho(n) ! density [kg/m^3] !! real(c_double), intent(out) :: rho_u(n) ! momentum [kg/m^2/s] !! real(c_double), intent(out) :: E(n) ! total energy [J/m^3] !! end subroutine !! !! `bind(C)` requires C-interoperable argument types: use `c_int` / `c_double` !! from `iso_c_binding` regardless of the solver's `wp` kind parameter. !! See `example/udf_sod.f90` for a complete working example. !! !! @param[inout] state Solver state (must have ub allocated, schemes bound) !! @param[in] cfg Runtime configuration (uses cfg%ic_udf_src) subroutine ic_udf(state, cfg) type(solver_state_t), intent(inout) :: state type(config_t), intent(in) :: cfg integer(c_int) :: host_platform, udf_status integer :: istat, i, n character(len=1536) :: cmd character(len=512) :: lib_path, msg, temp_dir character(len=20) :: pid_str integer(c_int) :: pid real(c_double), allocatable :: x(:), rho(:), rho_u(:), e_arr(:) ! ---- 1. Compile source to a shared library ---- ! Use a per-process (PID-unique) library name on every platform so that ! concurrent UDF compiles on one machine — e.g. several self-hosted CI ! runner services sharing /tmp — never collide on the same output file. host_platform = get_host_platform() pid = cfd_solver_get_pid() write (pid_str, '(I0)') pid select case (host_platform) case (host_platform_windows) temp_dir = get_windows_temp_dir() lib_path = join_path(trim(temp_dir), 'cfd_solver_udf_'//trim(pid_str)//'.dll', '\') cmd = 'gfortran -shared -o '//quote_shell_arg(trim(lib_path))//' '//quote_shell_arg(trim(cfg % ic_udf_src)) case (host_platform_linux) lib_path = '/tmp/cfd_solver_udf_'//trim(pid_str)//'.so' cmd = 'gfortran -shared -fPIC -o '//quote_shell_arg(trim(lib_path))//' '//quote_shell_arg(trim(cfg % ic_udf_src)) case (host_platform_macos) lib_path = '/tmp/cfd_solver_udf_'//trim(pid_str)//'.dylib' cmd = 'gfortran -dynamiclib -o '//quote_shell_arg(trim(lib_path))//' '//quote_shell_arg(trim(cfg % ic_udf_src)) case default error stop 'initial_conditions: unsupported host platform for runtime UDF initial conditions' end select call execute_command_line(trim(cmd), exitstat=istat) if (istat /= 0) & error stop 'initial_conditions: UDF compilation failed — check ic_udf_src path and syntax' ! ---- 2. Build grid-point coordinate array (same convention as all preset ICs) ---- n = state % n_pt allocate (x(n), rho(n), rho_u(n), e_arr(n), stat=istat) if (istat /= 0) error stop 'initial_conditions: allocate failed in ic_udf' do i = 1, n x(i) = real(state % mesh % x_node(i), c_double) end do ! ---- 3. Load the shared library, resolve the symbol, and call the UDF ---- udf_status = cfd_solver_call_udf(trim(lib_path)//c_null_char, int(n, c_int), x, rho, rho_u, e_arr) select case (udf_status) case (0_c_int) continue case (1_c_int) error stop 'initial_conditions: runtime UDF loader failed — cannot load shared library' case (2_c_int) error stop 'initial_conditions: runtime UDF loader failed — "ic_udf" symbol not found in shared library' case (3_c_int) error stop 'initial_conditions: runtime UDF loader failed — cannot close shared library' case default error stop 'initial_conditions: unknown runtime UDF loader error' end select ! ---- 4. Pack into conserved state array ---- do i = 1, n state % ub(1, i) = rho(i) state % ub(2, i) = rho_u(i) state % ub(3, i) = e_arr(i) end do deallocate (x, rho, rho_u, e_arr, stat=istat) if (istat /= 0) error stop 'initial_conditions: deallocate failed in ic_udf' ! ---- 5. Seed boundary halos from the outermost cells. ---- call seed_dirichlet_halos(state, state % ub(:, 1), state % ub(:, state % n_pt)) write (msg, '(A,A,A)') 'UDF IC loaded from "', trim(cfg % ic_udf_src), '"' call log_info(trim(msg)) end subroutine ic_udf end module initial_conditions