Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add: include new +turbulence subpackage #1004

Merged
merged 7 commits into from
Nov 22, 2024
Merged

Add: include new +turbulence subpackage #1004

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ec71dc2
Solve: typo
AlbertoCuadra Nov 11, 2024
28fbd7a
Add: include VelocityField object
AlbertoCuadra Nov 22, 2024
fb2a469
Add: include method to convert a struct or a 4D matrix into a Velocit…
AlbertoCuadra Nov 22, 2024
d4caa06
Add: new routines into the +math subpackge needed for the new +turbul…
AlbertoCuadra Nov 22, 2024
30987d2
Add: include new HelmholtzSolver class to perform Helmholtz decomposi…
AlbertoCuadra Nov 22, 2024
8fced25
Add: include new TurbulenceSpectra class to calculate the energy spec…
AlbertoCuadra Nov 22, 2024
57cac31
Add: include example demonstrating usage of `HelmholtzSolver` and `Tu…
AlbertoCuadra Nov 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 37 additions & 1 deletion +combustiontoolbox/+common/@Units/Units.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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
341 changes: 341 additions & 0 deletions +combustiontoolbox/+turbulence/@HelmholtzSolver/HelmholtzSolver.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Loading
Loading