diff --git a/+combustiontoolbox/+common/@Units/Units.m b/+combustiontoolbox/+common/@Units/Units.m index 5be54caf3..446c084ef 100644 --- a/+combustiontoolbox/+common/@Units/Units.m +++ b/+combustiontoolbox/+common/@Units/Units.m @@ -31,7 +31,7 @@ import combustiontoolbox.common.Units % Get the conversion factor property name - conversion_factor_name = [unit_in,'2',unit_out]; + conversion_factor_name = [unit_in, '2', unit_out]; % Check conversion factor exist assert(isprop(Units, conversion_factor_name), 'Conversion from %s to %s is not defined.', unit_in, unit_out); @@ -69,6 +69,42 @@ moles = weightPercentage ./ W; end + function velocity = convertData2VelocityField(data) + % Convert the input to a VelocityField object + % + % Args: + % data: Either a VelocityField object, a struct, or a 4D matrix + % + % Returns: + % velocity (VelocityField): VelocityField object + + % Import packages + import combustiontoolbox.turbulence.VelocityField + + if isa(data, 'VelocityField') + % Input is already a VelocityField object + velocity = data; + return + end + + if isstruct(data) && all(isfield(data, {'u', 'v', 'w'})) + % Input is a struct; convert to VelocityField + velocity = VelocityField(data.u, data.v, data.w); + return; + end + + if ndims(data) == 4 && size(data, 4) == 3 + % Input is a 4D matrix; convert to VelocityField + velocity = VelocityField(data(:, :, :, 1), ... + data(:, :, :, 2), ... + data(:, :, :, 3)); + return; + end + + % Unsupported format + error('Unsupported velocity input format. Expected a VelocityField object, a struct, or a 4D matrix.'); + end + end end \ No newline at end of file diff --git a/+combustiontoolbox/+turbulence/@HelmholtzSolver/HelmholtzSolver.m b/+combustiontoolbox/+turbulence/@HelmholtzSolver/HelmholtzSolver.m new file mode 100644 index 000000000..201314f34 --- /dev/null +++ b/+combustiontoolbox/+turbulence/@HelmholtzSolver/HelmholtzSolver.m @@ -0,0 +1,341 @@ +classdef HelmholtzSolver < handle + % The :mat:func:`HelmholtzSolver` class computes the Helmholtz-Hodge + % decomposition of a 3D velocity field into its solenoidal and + % dilatational parts using fast Fourier transform [1]. + % + % The decomposition is performed with the spectral method, which is + % only suitable for relatively smooth fields, i.e., with little power + % on small scales. The code assumes that the grid is uniform with + % dx = dy = dz. + % + % This code is based on Ref. [2] and has been rewritten in MATLAB + % with some modifications. + % + % Notes: + % For even NX, NY, and NZ, decomposed fields can be complex, + % with the imaginary part coming from the real part of the kmode + % at Nyquist frequency. In principle, the Nyquist frequency + % kmode should be dropped when doing the first derivatives to + % maintain symmetry. See footnote on page 4 of [2]. However, + % when the field is smooth enough, the imaginary part caused by + % the Nyquist frequency kmode should be negligible. + % + % Args: + % file_location (char): Path to the data .hdf file + % file_location_nodes (char): Path to the grid .hdf file + % T_ref (float): Temperature of reference [K] + % mu_ref (float): Dynamic viscosity of reference [kg/(m-s)] or [Pa-s] + % + % Example: + % solver = HelmholtzSolver(); + % + % References: + % [1] Johnson, S. G. (2011). Notes on FFT-based differentiation. + % MIT Applied Mathematics, Tech. Rep. + % Available: http://math.mit.edu/~stevenj/fft-deriv.pdf + % [2] Xun Shi, Helmholtz-Hodge decomposition using fft (Python), + % Available: https://github.com/shixun22/helmholtz + % + % + % @author: Alberto Cuadra Lara + % Postdoctoral researcher - Group Fluid Mechanics + % Universidad Carlos III de Madrid + + properties + tol0 = 1e-3 % Tolerance for checks + FLAG_CHECKS = true % Flag to perform checks + FLAG_TIME = true % Flag to print elapsed time + FLAG_REPORT = false % Flag to print predefined plots + time % Elapsed time + plotConfig % PlotConfig object + end + + methods + + function obj = HelmholtzSolver(varargin) + % Constructor + defaultPlotConfig = combustiontoolbox.utils.display.PlotConfig(); + + % Parse input arguments + p = inputParser; + addParameter(p, 'tol0', obj.tol0, @(x) isnumeric(x)); + addParameter(p, 'FLAG_CHECKS', obj.FLAG_CHECKS, @(x) islogical(x)); + addParameter(p, 'FLAG_TIME', obj.FLAG_TIME, @(x) islogical(x)); + addParameter(p, 'FLAG_REPORT', obj.FLAG_REPORT, @(x) islogical(x)); + addParameter(p, 'plotConfig', defaultPlotConfig, @(x) isa(x, 'combustiontoolbox.utils.display.PlotConfig')); + parse(p, varargin{:}); + + % Set properties + obj.tol0 = p.Results.tol0; + obj.FLAG_CHECKS = p.Results.FLAG_CHECKS; + obj.FLAG_TIME = p.Results.FLAG_TIME; + obj.FLAG_REPORT = p.Results.FLAG_REPORT; + obj.plotConfig = p.Results.plotConfig; + end + + function obj = set(obj, property, value, varargin) + % Set properties of the HelmholtzSolver object + % + % Args: + % obj (HelmholtzSolver): HelmholtzSolver object + % property (char): Property name + % value (float): Property value + % + % Optional Args: + % * property (char): Property name + % * value (float): Property value + % + % Returns: + % obj (HelmholtzSolver): HelmholtzSolver object with updated properties + % + % Examples: + % * set(HelmholtzSolver(), 'tol0', 1e-10) + % * set(HelmholtzSolver(), 'tol0', 1e-10, 'FLAG_CHECKS', false) + + varargin = [{property, value}, varargin{:}]; + + for i = 1:2:length(varargin) + % Assert that the property exists + assert(isprop(obj, varargin{i}), 'Property not found'); + + % Set property + obj.(varargin{i}) = varargin{i + 1}; + end + + end + + function [solenoidal, dilatational, STOP] = solve(obj, velocity, varargin) + % Compute Helmholtz decomposition of the velocity field + % + % Args: + % obj (HelmholtzSolver): HelmholtzSolver object + % velocity (VelocityField): Velocity field as a VelocityField object, struct, or 4D matrix + % + % Optional Args: + % * rho (float): Density field + % + % Returns: + % solenoidal (VelocityField): Struct with fields (u, v, w) containing the solenoidal velocity components + % dilatational (VelocityField): Struct with fields (u, v, w) containing the dilatational velocity components + % STOP (float): Relative error doing the decomposition + % + % Examples: + % * [solenoidal, dilatational, STOP] = solve(obj, velocity) + % * [solenoidal, dilatational, STOP] = solve(obj, velocity, 'rho', rho) + + % Import packages + import combustiontoolbox.common.Units.convertData2VelocityField + + % Timer + obj.time = tic; + + % Parse input arguments + p = inputParser; + addOptional(p, 'rho', [], @(x) isnumeric(x)); + parse(p, varargin{:}); + + % Set properties + rho = p.Results.rho; + + % Reshape velocity input and compute fluctuations + velocity = convertData2VelocityField(velocity); + velocity = obj.computeFluctuations(velocity, rho); + + % Solve the Helmholtz equation + [solenoidal, dilatational] = obj.decomposition(velocity); + + % Check results + STOP = obj.check(velocity, solenoidal, dilatational); + + % Time elapsed + obj.time = toc(obj.time); + + % Print elapsed time + printTime(obj); + end + + function [solenoidal, dilatational] = decomposition(obj, velocity) + % Compute Helmholtz decomposition of the velocity field + % + % Args: + % velocity (VelocityField): VelocityField instance with fields (u, v, w) containing the velocity components + % + % Returns: + % solenoidal (VelocityField): VelocityField instance with fields (u, v, w) containing the solenoidal velocity components + % dilatational (VelocityField): VelocityField instance with fields (u, v, w) containing the dilatational velocity components + + % Import packages + import combustiontoolbox.turbulence.VelocityField + + % Get N-D fast Fourier transform (fft) + U = fftn(velocity.u); + V = fftn(velocity.v); + W = fftn(velocity.w); + + % Compute wave numbers + [KX, KY, KZ] = obj.computeWaveNumbers(size(velocity.u)); + + % Compute k^2, avoiding division by zero + K2 = KX.^2 + KY.^2 + KZ.^2; + K2(K2 == 0) = 1; + + % Compute velocity divergence + div = (U .* KX + V .* KY + W .* KZ); + clear U V W % Free memory + + % Compute the Helmholtz decomposition (dilatational) + H = div ./ K2; + clear div K2 % Free memory + + % Compute dilatational components + dilatational = VelocityField(real(ifftn(H .* KX)), ... + real(ifftn(H .* KY)), ... + real(ifftn(H .* KZ))); + clear H KX KY KZ; % Free memory + + % Compute solenoidal components + solenoidal = VelocityField(velocity.u - dilatational.u, ... + velocity.v - dilatational.v, ... + velocity.w - dilatational.w); + end + + function [STOP, status] = check(obj, velocity, solenoidal, dilatational) + % Check the results of the Helmholtz decomposition + % + % Args: + % obj (HelmholtzSolver): HelmholtzSolver object + % velocity (VelocityField): VelocityField instance with fields (u, v, w) containing the velocity components + % solenoidal (VelocityField): VelocityField instance with fields (u, v, w) containing the solenoidal velocity components + % dilatational (VelocityField): VelocityField instance with fields (u, v, w) containing the dilatational velocity components + % + % Returns: + % STOP (float): Relative error doing the decomposition + % status (bool): Flag indicating if the checks passed + + if ~obj.FLAG_CHECKS + status = []; + STOP = []; + return + end + + % Validate Helmholtz decomposition results + fprintf('Performing checks... '); + + % Compute wave numbers + [KX, KY, KZ] = obj.computeWaveNumbers(size(velocity.u)); + + % Compute magnitude of the original velocity field for relative error + velocityMagnitude = globalMagnitude(velocity); + + % Check if the solenoidal part is divergence-free + divSolenoidal = ifftn((fftn(solenoidal.u) .* KX + ... + fftn(solenoidal.v) .* KY + ... + fftn(solenoidal.w) .* KZ)); + + divError = max(abs(divSolenoidal(:))) / velocityMagnitude; + + if divError > obj.tol0 + fprintf('Error!\nSolenoidal part divergence too high: %.2e\n', divError); + STOP = divError; + return; + end + + % Check if the dilatational part is curl-free + curlDilatational = ifftn((fftn(dilatational.w) .* KY - fftn(dilatational.v) .* KZ + ... + fftn(dilatational.u) .* KZ - fftn(dilatational.w) .* KX + ... + fftn(dilatational.v) .* KX - fftn(dilatational.u) .* KY) * 1i * 2 * pi); + + curlError = max(abs(curlDilatational(:))) / velocityMagnitude; + + if curlError > obj.tol0 + fprintf('Error!\nDilatational part curl too high: %.2e\n', curlError); + STOP = curlError; + return; + end + + % Check if the solenoidal and dilatational parts sum up to the original field + uRecon = solenoidal.u + dilatational.u; + vRecon = solenoidal.v + dilatational.v; + wRecon = solenoidal.w + dilatational.w; + + diffError = sqrt(sum((uRecon(:) - velocity.u(:)).^2 + ... + (vRecon(:) - velocity.v(:)).^2 + ... + (wRecon(:) - velocity.w(:)).^2)) / velocityMagnitude; + + if diffError > obj.tol0 + fprintf('Error!\nReconstructed field does not match original field: %.2e\n', diffMax); + STOP = diffMax; + return; + end + + % Set flag to pass and stop + STOP = max([divError, curlError, diffError]); + status = true; + fprintf('OK!\n'); + end + + function printTime(obj) + % Print execution time + % + % Args: + % obj (EquilibriumSolver): Object of the class EquilibriumSolver + + if ~obj.FLAG_TIME + return + end + + % Definitions + operationName = 'Helmholtz decomposition'; + + % Print elapsed time + fprintf('\nElapsed time for %s: %.5f seconds\n', operationName, obj.time); + end + + end + + methods (Static) + + function velocity = computeFluctuations(velocity, rho) + % Compute fluctuating velocity components + % + % Args: + % velocity (struct): Struct with fields (u, v, w) + % rho (float): Density field + % + % Returns: + % velocity (struct): Struct with fields (u, v, w) containing the fluctuating velocity components + + % For compressible flows + if ~isempty(rho) + rhoMean = mean(rho, 'all'); + rhou = mean(rho .* velocity.u, 'all') / rhoMean; + rhov = mean(rho .* velocity.v, 'all') / rhoMean; + rhow = mean(rho .* velocity.w, 'all') / rhoMean; + + velocity.u = sqrt(rho) .* (velocity.u - rhou); + velocity.v = sqrt(rho) .* (velocity.v - rhov); + velocity.w = sqrt(rho) .* (velocity.w - rhow); + return + end + + % For incompressible flows + velocity.u = velocity.u - mean(velocity.u, 'all'); + velocity.v = velocity.v - mean(velocity.v, 'all'); + velocity.w = velocity.w - mean(velocity.w, 'all'); + end + + function [KX, KY, KZ] = computeWaveNumbers(sz) + % Compute wave number grids for FFT + + % Import packages + import combustiontoolbox.utils.math.fftfreq + + kx = fftfreq(sz(1)); + ky = fftfreq(sz(2)); + kz = fftfreq(sz(3)); + [KX, KY, KZ] = ndgrid(kx, ky, kz); + end + + end + +end \ No newline at end of file diff --git a/+combustiontoolbox/+turbulence/@TurbulenceSpectra/TurbulenceSpectra.m b/+combustiontoolbox/+turbulence/@TurbulenceSpectra/TurbulenceSpectra.m new file mode 100644 index 000000000..da5b7b0d6 --- /dev/null +++ b/+combustiontoolbox/+turbulence/@TurbulenceSpectra/TurbulenceSpectra.m @@ -0,0 +1,310 @@ +classdef TurbulenceSpectra < handle + % The :mat:func:`TurbulenceSpectra` class provides methods for computing turbulence spectra. + % Supports spherical and cross-plane averaging. + + properties + averaging = 'spherical' % Type of averaging ('spherical', 'crossplane') + axis = 'x' % Axis for cross-plane averaging ('x', 'y', or 'z') + time % Elapsed time + plotConfig % PlotConfig object + FLAG_TIME = true % Flag to print elapsed time + FLAG_REPORT = false % Flag to print predefined plots + end + + methods + + function obj = TurbulenceSpectra(varargin) + % Constructor for TurbulenceSpectra + % + % Optional Args:d + % averaging (char): Type of averaging ('spherical' or 'crossplane') + % axis (char): Axis for cross-plane averaging ('x', 'y', 'z') + % + % Example: + % analyzer = TurbulenceSpectra('averaging', 'crossplane', 'axis', 'z'); + + % Default + defaultPlotConfig = combustiontoolbox.utils.display.PlotConfig(); + defaultPlotConfig.xscale = 'log'; defaultPlotConfig.yscale = 'log'; + defaultPlotConfig.innerposition = [0.15 0.15 0.35 0.5]; + defaultPlotConfig.outerposition = [0.15 0.15 0.35 0.5]; + + % Parse input arguments + p = inputParser; + addParameter(p, 'averaging', obj.averaging, @(x) ischar(x) && ismember(lower(x), {'spherical', 'crossplane'})); + addParameter(p, 'axis', obj.axis, @(x) ischar(x) && ismember(lower(x), {'x', 'y', 'z'})); + addParameter(p, 'plotConfig', defaultPlotConfig, @(x) isa(x, 'combustiontoolbox.utils.display.PlotConfig')); + addParameter(p, 'FLAG_TIME', obj.FLAG_TIME, @islogical); + addParameter(p, 'FLAG_REPORT', obj.FLAG_REPORT, @islogical); + parse(p, varargin{:}); + + obj.averaging = lower(p.Results.averaging); + obj.axis = lower(p.Results.axis); + obj.plotConfig = p.Results.plotConfig; + obj.FLAG_TIME = p.Results.FLAG_TIME; + obj.FLAG_REPORT = p.Results.FLAG_REPORT; + end + + function obj = set(obj, property, value, varargin) + % Set properties of the TurbulenceSpectra object + % + % Args: + % obj (TurbulenceSpectra): TurbulenceSpectra object + % property (char): Property name + % value (float): Property value + % + % Optional Args: + % * property (char): Property name + % * value (float): Property value + % + % Returns: + % obj (TurbulenceSpectra): TurbulenceSpectra object with updated properties + % + % Examples: + % * set(TurbulenceSpectra(), 'averaging', 'crossplane'); + % * set(TurbulenceSpectra(), 'averaging', 'crossplane', 'axis', 'z'); + % * set(TurbulenceSpectra(), 'averaging', 'crossplane', 'axis', 'z', 'FLAG_TIME', false); + + varargin = [{property, value}, varargin{:}]; + + for i = 1:2:length(varargin) + % Assert that the property exists + assert(isprop(obj, varargin{i}), 'Property not found'); + + % Set property + obj.(varargin{i}) = varargin{i + 1}; + end + + end + + function [EK_avg, k] = getEnergySpectra(obj, velocity) + % Compute the energy spectra of a 3D velocity field. + % + % Args: + % velocity (VelocityField): Velocity field as a VelocityField object, struct, or 4D matrix + % + % Returns: + % EK_avg (float): Averaged energy spectra + % k (float): Wavenumber vector + + % Import packages + import combustiontoolbox.common.Units.convertData2VelocityField + + % Start timer + obj.time = tic; + + % Convert input to VelocityField object + velocity = convertData2VelocityField(velocity); + + % Set the appropriate averaging function + switch obj.averaging + case 'spherical' + % Use 3D FFT for spherical averaging + U = fftn(velocity.u) / numel(velocity.u); + V = fftn(velocity.v) / numel(velocity.v); + W = fftn(velocity.w) / numel(velocity.w); + EK = 0.5 * fftshift(abs(U).^2 + abs(V).^2 + abs(W).^2); + % Set the function for spherical averaging + averagingFunction = @obj.getSphericallyAveragedSpectra; + case 'crossplane' + % Use 2D FFT for cross-plane averaging + EK = obj.getCrossplaneFFT(velocity, obj.axis); + % Set the function for cross-plane averaging + averagingFunction = @obj.getCrossplaneAveragedSpectra; + end + + % Compute the energy spectra based on the averaging type + [EK_avg, k] = averagingFunction(EK); + + % Time elapsed + obj.time = toc(obj.time); + + % Print elapsed time + printTime(obj); + end + + function ax = plot(obj, k, EK, varargin) + % Plot results + % + % Args: + % obj (TurbluenceSpectra): TurbulenceSpectra object + % k (float): Wavenumber vector + % EK (float): Energy spectra + % + % Optional Args: + % * EK_i (float): Array of energy spectra + % + % Returns: + % ax (Axes): Axes object + % + % Examples: + % * ax = plot(TurbulenceSpectra(), k, EK); + % * ax = plot(TurbulenceSpectra(), k, EK1, EK2); + + % Import packages + import combustiontoolbox.utils.display.* + + % Definitions + numCases = length(varargin) + 1; + + % Plot spectra + ax = setFigure(obj.plotConfig); + plotFigure('k', k, 'E(k)', EK(1:length(k)), 'color', 'auto', 'ax', ax); + + for i = 1:numCases-1 + EK = varargin{i}; + plotFigure('k', k, 'E(k)', EK(1:length(k)), 'color', 'auto', 'ax', ax); + end + + % legend(ax, {'solenoidal', 'dilatational', 'total'}, 'interpreter', 'latex'); + % plotFigure('k', k, 'E(k)', 0.1 * k.^(-5/3), 'linestyle', '--', 'color', [0 0 0], 'ax', ax); + end + + function printTime(obj) + % Print execution time + % + % Args: + % obj (EquilibriumSolver): Object of the class EquilibriumSolver + + if ~obj.FLAG_TIME + return + end + + % Definitions + operationName = 'energy spectra'; + + % Print elapsed time + fprintf('\nElapsed time for %s: %.5f seconds\n', operationName, obj.time); + end + + end + + methods (Static, Access = private) + + function EK = getCrossplaneFFT(velocity, axisType) + % Compute 2D FFT on homogeneous planes based on axisType + % + % Args: + % u (float): 3D velocity field in the x-direction + % v (float): 3D velocity field in the y-direction + % w (float): 3D velocity field in the z-direction + % axisType (char): Axis for cross-plane averaging ('x', 'y', or 'z') + % + % Returns: + % EK (float): 3D array representing the 2D energy spectra computed on the + % homogeneous planes, sliced along the inhomogeneous axis. + % The first two dimensions correspond to the homogeneous plane, + % while the third dimension represents the slices along the + % inhomogeneous axis. + % + % Examples: + % * EK = getCrossplaneFFT(velocity, 'x'); + % * EK = getCrossplaneFFT(velocity, 'y'); + % * EK = getCrossplaneFFT(velocity, 'z'); + + % Map axis to slicing dimension + sliceDim = find('xyz' == axisType); + dims = 1:3; + homogeneousDims = setdiff(dims, sliceDim); + + % Permute dimensions to bring slice dimension last + permutedU = permute(velocity.u, [homogeneousDims, sliceDim]); + permutedV = permute(velocity.v, [homogeneousDims, sliceDim]); + permutedW = permute(velocity.w, [homogeneousDims, sliceDim]); + + % Extract slice sizes and initialize EK + [size1, size2, numSlices] = size(permutedU); + EK = zeros(size1, size2, numSlices); + + % Compute 2D FFT slice by slice + for i = 1:numSlices + % Extract slices + U = fft2(squeeze(permutedU(:, :, i))) / numel(permutedU(:, :, i)); + V = fft2(squeeze(permutedV(:, :, i))) / numel(permutedV(:, :, i)); + W = fft2(squeeze(permutedW(:, :, i))) / numel(permutedW(:, :, i)); + + % Compute energy spectra for the slice + EK(:, :, i) = 0.5 * fftshift(abs(U).^2 + abs(V).^2 + abs(W).^2); + end + + end + + function [EK_avg, k] = getCrossplaneAveragedSpectra(EK) + % Compute the cross-plane averaged energy spectra of a 3D field. + % The homogeneous plane are the first two dimensions of the 3D field. + % + % Args: + % EK (float): 3D energy spectra (combined from all velocity components) + % + % Returns: + % EK_avg (float): Radially averaged energy spectra for the slices + % k (float): Wavenumber vector along the homogeneous direction + % + % Example: + % [EK_avg, k] = getCrossplaneAveragedSpectra(EK); + + % Definitions + sliceDim = 3; % Axis for slicing + numSlices = size(EK, sliceDim); + + % Get the center of the homogeneous plane + [N1, N2, ~] = size(EK); + centerX1 = floor(N1 / 2); + centerX2 = floor(N2 / 2); + [X1, X2] = ndgrid(1:N1, 1:N2); + + % Compute radii for radial averaging + radii = round(sqrt((X1 - centerX1).^2 + (X2 - centerX2).^2)) + 1; + maxRadius = max(radii(:)); + + % Radially average slices + EK_avg = zeros(numSlices, maxRadius); + for i = 1:numSlices + % Extract slice + sliceData = squeeze(EK(:, :, i)); + + % Radially average + EK_avg(i, :) = accumarray(radii(:), sliceData(:), [maxRadius, 1], @mean, 1e-50); + end + + % Compute wavenumber vector + fftResult = fft(EK(:, 1, 1)); + numFreqs = ceil(length(fftResult) / 2); + k = 0:numFreqs - 1; + end + + function [EK_avg, k] = getSphericallyAveragedSpectra(EK) + % Compute the spherically averaged energy spectra of a 3D field + % + % Args: + % EK (float): 3D energy spectra + % + % Returns: + % EK_avg (float): Spherically averaged energy spectra + % + % Example: + % EK_avg = getSphericallyAveragedSpectra(EK); + + % Get the center of the 3D field + [NX, NY, NZ] = size(EK); + centerX = floor(NX / 2); + centerY = floor(NY / 2); + centerZ = floor(NZ / 2); + + % Compute radii for radial averaging + [X, Y, Z] = ndgrid(1:NX, 1:NY, 1:NZ); + radii = round(sqrt((X - centerX).^2 + (Y - centerY).^2 + (Z - centerZ).^2)) + 1; + maxRadius = max(radii(:)); + + % Radially average the energy spectra + EK_avg = accumarray(radii(:), EK(:), [maxRadius, 1], @sum, 1e-50); + + % Compute wavenumber vector + fftResult = fft(EK(:, 1, 1)); + numFreqs = ceil(length(fftResult) / 2); + k = 0:numFreqs - 1; + end + + end + +end diff --git a/+combustiontoolbox/+turbulence/@VelocityField/VelocityField.m b/+combustiontoolbox/+turbulence/@VelocityField/VelocityField.m new file mode 100644 index 000000000..fbe53deca --- /dev/null +++ b/+combustiontoolbox/+turbulence/@VelocityField/VelocityField.m @@ -0,0 +1,94 @@ +classdef VelocityField < handle + % This :mat:func:`VelocityField` class provides methods for computing properties of a velocity field. + + properties + u % x-component of the velocity field + v % y-component of the velocity field + w % z-component of the velocity field + end + + methods + + function obj = VelocityField(u, v, w) + % Constructor for VelocityField + % + % Args: + % u (float): x-component of velocity + % v (float): y-component of velocity + % w (float): z-component of velocity + % + % Returns: + % obj (VelocityField): VelocityField object + + % Parse inputs + p = inputParser; + addRequired(p, 'u', @isnumeric); + addRequired(p, 'v', @isnumeric); + addRequired(p, 'w', @isnumeric); + parse(p, u, v, w); + + % Assign properties + obj.u = p.Results.u; + obj.v = p.Results.v; + obj.w = p.Results.w; + end + + function magnitudeField = pointwiseMagnitude(obj) + % Compute the pointwise magnitude of the velocity field + % + % Args: + % obj (VelocityField): VelocityField object + % + % Returns: + % magnitudeField (float): 3D matrix of velocity magnitudes + % + % Example: + % magnitudeField = pointwiseMagnitude(obj); + + magnitudeField = sqrt(obj.u.^2 + obj.v.^2 + obj.w.^2); + end + + function globalMagnitude = globalMagnitude(obj) + % Compute the global magnitude (Frobenius norm) of the velocity field + % + % Args: + % obj (VelocityField): VelocityField object + % + % Returns: + % globalMagnitude (float): Scalar magnitude of the velocity field + % + % Example: + % globalMagnitude = globalMagnitude(obj); + + globalMagnitude = sqrt(sum(obj.u(:).^2 + obj.v(:).^2 + obj.w(:).^2)); + end + + function obj = normalize(obj) + % Normalize the velocity field by its magnitude + % + % Returns: + % obj: Normalized velocity field + + magnitude = obj.computeMagnitude(); + obj.u = obj.u ./ magnitude; + obj.v = obj.v ./ magnitude; + obj.w = obj.w ./ magnitude; + end + + function obj = scale(obj, factor) + % Scale the velocity field by a factor + % + % Args: + % factor: Scaling factor + % + % Returns: + % obj: Scaled velocity field + + obj.u = obj.u * factor; + obj.v = obj.v * factor; + obj.w = obj.w * factor; + end + + end + +end \ No newline at end of file diff --git a/+combustiontoolbox/+utils/+math/fftfreq.m b/+combustiontoolbox/+utils/+math/fftfreq.m new file mode 100644 index 000000000..7c77c7974 --- /dev/null +++ b/+combustiontoolbox/+utils/+math/fftfreq.m @@ -0,0 +1,40 @@ +function f = fftfreq(n, varargin) + % Return the Discrete Fourier Transform sample frequencies + % + % The returned float array `f` contains the frequency bin centers in cycles + % per unit of the sample spacing (with zero at the start). For instance, if + % the sample spacing is in seconds, then the frequency unit is cycles/second. + % + % Given a window length `n` and a sample spacing `d`: + % + % f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even + % f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd + % + % Args: + % n (float): Window length + % + % Optional Args: + % d (float): Sample spacing (inverse of the sampling rate) [default: 1] + % + % Returns: + % value (float): Array of length `n//2 + 1` containing the sample frequencies + + % Import packages + import combustiontoolbox.utils.math.floorDiv + + % Definitions + d = 1; + n = round(n); + + % Unpack additional inputs + if (nargin > 1) + d = varargin{1}; + end + + % Calculations + val = 1.0 / (n * d); + N = floorDiv(n - 1, 2) + 1; + p1 = 0:N - 1; + p2 = -floorDiv(n, 2):-1; + f = [p1, p2] * val; +end \ No newline at end of file diff --git a/+combustiontoolbox/+utils/+math/floorDiv.m b/+combustiontoolbox/+utils/+math/floorDiv.m new file mode 100644 index 000000000..b2e884184 --- /dev/null +++ b/+combustiontoolbox/+utils/+math/floorDiv.m @@ -0,0 +1,15 @@ +function x = floorDiv(a, b) + % Round the result of division toward negative infinity + % + % Args: + % a (float): Value a + % b (float): Value b + % + % Returns: + % x (float): Floor division + % + % Example: + % x = floorDiv(5, 2); + + x = floor(a ./ b); +end \ No newline at end of file diff --git a/+combustiontoolbox/+utils/+math/gradientPeriodic.m b/+combustiontoolbox/+utils/+math/gradientPeriodic.m new file mode 100644 index 000000000..1b059e059 --- /dev/null +++ b/+combustiontoolbox/+utils/+math/gradientPeriodic.m @@ -0,0 +1,29 @@ +function [grad_x, grad_y, grad_z] = gradientPeriodic(phi, x, y, z) + % Computes the gradient of a three-dimensional scalar field using + % second-order central finite differences considering a periodic box + % + % Args: + % phi (float): Scalar field to compute the gradient of + % x (float): 3D array with the x-coordinate of the grid + % y (float): 3D array with the y-coordinate of the grid + % z (float): 3D array with the z-coordinate of the grid + % + % Returns: + % grad_x (float): 3D array with the gradient of phi over x + % grad_y (float): 3D array with the gradient of phi over y + % grad_z (float): 3D array with the gradient of phi over z + % + % Example: + % [grad_x, grad_y, grad_z] = gradientPeriodic(phi, x, y, z); + + % Get grid spacing in the three directions + hx = x(2) - x(1); + hy = y(2) - y(1); + hz = z(2) - z(1); + + % Compute the gradient in each direction using central finite differences + % considering periodic boundary conditions + grad_x = (circshift(phi, [-1, 0, 0]) - circshift(phi, [1, 0, 0])) ./ (2 * hx); + grad_y = (circshift(phi, [ 0, -1, 0]) - circshift(phi, [0, 1, 0])) ./ (2 * hy); + grad_z = (circshift(phi, [ 0, 0, -1]) - circshift(phi, [0, 0, 1])) ./ (2 * hz); +end \ No newline at end of file diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index f931b475a..0691ecae6 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -11,7 +11,7 @@ To contribute to this repository, please follow these steps: 2. Clone it to your local machine. ```console -git clone https://github.com/your-github-username/combustion_toolbox.git +https://github.com/CombustionToolbox/combustion_toolbox.git ``` 3. Create a new branch for your changes: @@ -42,4 +42,4 @@ Please ensure that any MATLAB code you contribute follows the same code style as If you encounter a bug in this repository, please submit a detailed bug report to the Issues section of this repository. Please include steps to reproduce the bug, expected vs. actual behavior, and any relevant code snippets or screenshots. ## Contact Us -If you have any questions or concerns about contributing to this repository, please feel free to contact us at [acuadra@ing.uc3m.es](mailto:acuadra@ing.uc3m.es). We would be happy to help! \ No newline at end of file +If you have any questions or concerns about contributing to this repository, please feel free to contact us at [acuadra@ing.uc3m.es](mailto:acuadra@ing.uc3m.es). We would be happy to help! diff --git a/examples/Example_TURBULENCE_HELMHOLTZ_DECOMPOSTION.m b/examples/Example_TURBULENCE_HELMHOLTZ_DECOMPOSTION.m new file mode 100644 index 000000000..b744b199c --- /dev/null +++ b/examples/Example_TURBULENCE_HELMHOLTZ_DECOMPOSTION.m @@ -0,0 +1,55 @@ +% ------------------------------------------------------------------------- +% EXAMPLE: Helmholtz decomposition +% +% Perform Helmholtz decomposition of a velocity field and analyze the +% turbulence spectra of the original, solenoidal, and dilatational fields. +% +% Note: The data used is a downsampled velocity field from a Direct Numerical +% Simulation (DNS) of a Homogeneous Isotropic Turbulence (HIT) case carried +% out with the HTR solver [1]. Take the data as an example to illustrate +% the use of the HelmholtzSolver and TurbulenceSpectra classes. +% +% References: +% [1] Di Renzo, M., Fu, L., & Urzay, J. (2020). HTR solver: An +% open-source exascale-oriented task-based multi-GPU high-order +% code for hypersonic aerothermodynamics. Computer Physics +% Communications, 255, 107262. +% +% +% @author: Alberto Cuadra Lara +% Postdoctoral researcher - Group Fluid Mechanics +% Universidad Carlos III de Madrid +% +% Last update Nov 22 2024 +% ------------------------------------------------------------------------- + +% Import required packages +import combustiontoolbox.turbulence.* + +% Load field variables +filename = 'velocityField.h5'; +filePath = which(filename); +rho = double(h5read(filePath, ['/', 'density'])); +u = double(h5read(filePath, ['/', 'velocity/u'])); +v = double(h5read(filePath, ['/', 'velocity/v'])); +w = double(h5read(filePath, ['/', 'velocity/w'])); + +% Convert velocity to VelocityField object +velocityField = VelocityField(u, v, w); + +% Initialize HelmholtzSolver +solver = HelmholtzSolver(); + +% Perform Helmholtz decomposition +[solenoidal, dilatational, STOP] = solver.solve(velocityField, 'rho', rho); + +% Analyze turbulence spectra using TurbulenceSpectra +analyzer = TurbulenceSpectra(); + +% Compute energy spectra for the original, solenoidal, and dilatational fields +[EK1, k1] = analyzer.getEnergySpectra(velocityField); % Original field +[EK2, k2] = analyzer.getEnergySpectra(solenoidal); % Solenoidal component +[EK3, k3] = analyzer.getEnergySpectra(dilatational); % Dilatational component + +% Plot results +analyzer.plot(k1, EK1, EK2, EK3); \ No newline at end of file diff --git a/validations/htr/data/hit/velocityField.h5 b/validations/htr/data/hit/velocityField.h5 new file mode 100644 index 000000000..e649297b7 Binary files /dev/null and b/validations/htr/data/hit/velocityField.h5 differ