!> @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` and the ghost-state vectors !! (`state%q_left`, `state%q_right`, and the boundary split fluxes for FVS !! schemes) based on the `problem_type` field in the supplied configuration. !! !! 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 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 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. !! !! Paths containing a double-quote character are rejected with error stop !! because they cannot be safely represented in the generated command string. 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) & error stop 'initial_conditions: path contains a double-quote character and cannot be quoted safely' quoted_path = '"'//trim(raw_path)//'"' end function quote_shell_arg !> Fill `state%ub` with the initial condition selected by `cfg%problem_type`. !! !! Also sets `state%q_left`, `state%q_right`, and — for FVS schemes — !! pre-computes the boundary split fluxes `fl_pos`, `fl_neg`, `fr_pos`, !! `fr_neg`. !! !! @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) :: state type(config_t), intent(in) :: cfg real(wp), parameter :: pi = 4.0_wp * atan(1.0_wp) integer :: ipt real(wp) :: x, rho_so select case (trim(cfg % problem_type)) case (problem_sod) ! Sod (1978) shock tube. Left/right primitives from namelist. ! Reference: Sod, J. Comput. Phys. 27:1-31, 1978. call set_left_right_states(state, cfg) do ipt = 1, state % n_pt x = state % cfg % x_left + state % dx * real(ipt - 1, wp) if (x < cfg % x_diaphragm) then state % ub(:, ipt) = state % q_left else if (x > cfg % x_diaphragm) then state % ub(:, ipt) = state % q_right else state % ub(:, ipt) = 0.5_wp * (state % q_left + state % q_right) end if end do 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) state % q_left(1) = 3.857143_wp state % q_left(2) = 3.857143_wp * 2.629369_wp state % q_left(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). ! At t=1.8 the rightmost wave has not reached x=+5. state % q_right(1) = 1.0_wp state % q_right(2) = 0.0_wp state % q_right(3) = 1.0_wp / (state % cfg % gam - 1.0_wp) call precompute_boundary_fluxes(state) do ipt = 1, state % n_pt x = state % cfg % x_left + state % dx * real(ipt - 1, wp) if (x < -4.0_wp) then ! Post-shock region: constant state state % ub(:, ipt) = state % q_left else ! Pre-shock region: sinusoidal density perturbation 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 case (problem_smooth_wave) ! Smooth advecting density wave on periodic domain [0, 1]. ! IC: rho = 1 + 0.2*sin(2*pi*x), u = 1, p = 1. ! Exact: rho(x,t) = 1 + 0.2*sin(2*pi*(x - t)), u = 1, p = 1. ! Ghost state vectors (fl_pos etc.) are unused for periodic BCs. do ipt = 1, state % n_pt x = state % cfg % x_left + state % dx * real(ipt - 1, wp) 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 case (problem_linear_advection) ! Sinusoidal density wave on a uniform background flow, periodic domain. ! The advection velocity u_adv is taken from cfg%u_left. ! IC: rho(x) = 1 + 0.5*sin(2*pi*(x - x_left)/L), u = u_adv, p = 1 ! Ghost state vectors (fl_pos etc.) are unused for periodic BCs. do ipt = 1, state % n_pt x = state % cfg % x_left + state % dx * real(ipt - 1, wp) 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: ! Region 1 (0 <= x < 0.1): p = 1000 ! Region 2 (0.1 <= x < 0.9): p = 0.01 ! Region 3 (0.9 <= x <= 1): p = 100 ! BCs: reflecting walls. state % q_left(1) = 1.0_wp state % q_left(2) = 0.0_wp state % q_left(3) = 1000.0_wp / (state % cfg % gam - 1.0_wp) state % q_right(1) = 1.0_wp state % q_right(2) = 0.0_wp state % q_right(3) = 100.0_wp / (state % cfg % gam - 1.0_wp) call precompute_boundary_fluxes(state) do ipt = 1, state % n_pt x = state % cfg % x_left + state % dx * real(ipt - 1, wp) 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_lax) ! Lax (1954) shock tube. ! Left: rho=0.445, u=0.698, p=3.528 ! Right: rho=0.5, u=0, p=0.571 call set_left_right_states(state, cfg) do ipt = 1, state % n_pt x = state % cfg % x_left + state % dx * real(ipt - 1, wp) if (x < cfg % x_diaphragm) then state % ub(:, ipt) = state % q_left else if (x > cfg % x_diaphragm) then state % ub(:, ipt) = state % q_right else state % ub(:, ipt) = 0.5_wp * (state % q_left + state % q_right) end if end do case (problem_from_file) ! Read primitives (x, rho, u, p) from cfg%ic_file and interpolate if needed. call ic_from_file(state, cfg) case (problem_udf) ! Compile cfg%ic_udf_src to a shared library at runtime and call it. call ic_udf(state, cfg) case (problem_acoustic_pulse) ! Gaussian acoustic pulse on a quiescent uniform background. ! Canonical benchmark for non-reflecting boundary conditions. ! Reference: Thompson (1987), J. Comput. Phys. 68, 1-24; ! Bogey & Bailly (2002), Acta Acustica 88, 463-471. ! ! Background: rho0 = 1, u0 = 0, p0 = 1 (gamma = 1.4, c0 = sqrt(gamma)) ! Perturbation: Gaussian pressure pulse of amplitude dp = 1e-3 and ! half-width parameter alpha = 200, centred at x = 0.5. ! Density perturbation is isentropic: drho = dp/c0^2. ! Velocity perturbation is zero (the pulse splits into symmetric L/R waves). 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 = state % cfg % x_left + state % dx * real(ipt - 1, wp) 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 ! Ghost states: set to uniform background (non-reflecting BC will handle them). state % q_left = [rho0, 0.0_wp, p0 / (state % cfg % gam - 1.0_wp)] state % q_right = [rho0, 0.0_wp, p0 / (state % cfg % gam - 1.0_wp)] call precompute_boundary_fluxes(state) end block case default error stop 'initial_conditions: unknown problem_type' end select ! Prime ghost states for dynamically-computed BC types so that any FVS ! split fluxes at the ghost reflect the correct physical state before the ! first call to compute_resid. ! apply_bcs dispatches each side independently (side='L'/'R'), so any ! combination of left and right BC types is handled correctly: ! 'dirichlet'/'inflow'/'supersonic_inlet' — no-op (early return) ! 'neumann'/'neumann_gradient'/'subsonic_inlet'/'subsonic_outlet' — ghost updated call apply_bcs(state) call precompute_boundary_fluxes(state) end subroutine apply_initial_condition ! --------------------------------------------------------------------------- ! Private helpers ! --------------------------------------------------------------------------- !> Convert namelist primitives (rho_left, u_left, p_left / rho_right, …) to !! conserved ghost state vectors and optionally pre-compute boundary split !! fluxes for FVS schemes. subroutine set_left_right_states(state, cfg) type(solver_state_t), intent(inout) :: state type(config_t), intent(in) :: cfg state % q_left(1) = cfg % rho_left state % q_left(2) = cfg % rho_left * cfg % u_left state % q_left(3) = cfg % p_left / (state % cfg % gam - 1.0_wp) & + 0.5_wp * cfg % rho_left * cfg % u_left**2 state % q_right(1) = cfg % rho_right state % q_right(2) = cfg % rho_right * cfg % u_right state % q_right(3) = cfg % p_right / (state % cfg % gam - 1.0_wp) & + 0.5_wp * cfg % rho_right * cfg % u_right**2 call precompute_boundary_fluxes(state) end subroutine set_left_right_states !> Pre-compute the boundary split fluxes fl_pos/fl_neg/fr_pos/fr_neg for FVS !! schemes. No-op when the FDS path is active (`state%use_fds = .true.`). subroutine precompute_boundary_fluxes(state) type(solver_state_t), intent(inout) :: state if (.not. state % use_fds) then call state % flux_split(state % q_left, state % fl_pos, 1, state % cfg % gam) call state % flux_split(state % q_left, state % fl_neg, -1, state % cfg % gam) call state % flux_split(state % q_right, state % fr_pos, 1, state % cfg % gam) call state % flux_split(state % q_right, state % fr_neg, -1, state % cfg % gam) end if end subroutine precompute_boundary_fluxes !> 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`, ghost states are set from the outermost loaded cells and !! `precompute_boundary_fluxes` is called. !! !! @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 % cfg % x_left + state % dx * real(i - 1, wp) 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 % cfg % x_left + state % dx * real(ipt - 1, wp) 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. Ghost states from outermost loaded cells --- state % q_left = state % ub(:, 1) state % q_right = state % ub(:, state % n_pt) call precompute_boundary_fluxes(state) 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. !! (The default `wp=real64` happens to match `c_double`, but this is not !! guaranteed if `wp` is changed to `real32` or `real128`; the UDF ABI is !! fixed to `c_int`/`c_double` and must not rely on `wp`.) This also !! suppresses gfortran -Wc-binding-type warnings. !! See `example/udf_sod.f90` for a complete working example. !! !! The source is compiled via `execute_command_line` and the resulting shared !! library is written to `/tmp/` on POSIX hosts or `%TEMP%` / `%TMP%` on !! Windows. Linux uses `-shared -fPIC` with a `.so` suffix; macOS uses !! `-dynamiclib` with `.dylib`; Windows uses `-shared` with `.dll`. The !! library is loaded through a small C portability shim that calls !! `dlopen`/`dlsym` on POSIX hosts and `LoadLibrary`/`GetProcAddress` on !! Windows. Requires gfortran on `PATH`. !! !! @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 ---- host_platform = get_host_platform() select case (host_platform) case (host_platform_windows) pid = cfd_solver_get_pid() write (pid_str, '(I0)') pid 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.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.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 % cfg % x_left + state % dx * real(i - 1, wp), 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. Ghost states from outermost cells ---- state % q_left = state % ub(:, 1) state % q_right = state % ub(:, state % n_pt) call precompute_boundary_fluxes(state) 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