% KIRCHHOFF_CD_SIMULATION.M Kirchhoff Ellipse Contour Dynamics Simulation Program
%
% Version: 1.0, 19 March 2008
% Authors: Travis B. Mitchell and Louis F. Rossi
% Email: travis*b*mitchell=gmail*com, rossi=math*udel*edu
% Real_email = regexprep(Email,{'=','*'},{'@','.'})
%
% This code was used in 'The Evolution of Kirchhoff Elliptic Vortices' by T.B. Mitchell and L. F. Rossi,
% which will appear in Physics of Fluids, April 2008. A preprint is included in the readme file.
%
% The routines 1) create a perturbed Kirchhoff vortex, 2) follow its evolution forward in time using a
% contour dynamics algorithm, then 3) fit the contours to modes and write the results.
%
% INITIALIZATION: Kirchhoff ellipse base state is specified (a, b, zeta), then one mode
% is applied. Depending on i_pert_control, seed mode can be specified as per equations
% (15a), (15b), or you can specify the amplitude of a purely growing or decaying eigenmode
% where the phase will be set as per equation (21).
%
% EVOLUTION: Most accurate option is to run to desired time, however if contour points get too
% close or too far, simulation will hang. So it is often useful to break the run into smaller
% steps to monitor things, and the run can later be done in one go once the largest time of
% interest is known. Options are set with i_integrator_control.
%
% FITTING: Two fitting approaches, described in section IIE, are implemented. The 'Linear'
% one is faster since it fits to fewer quantities. The 'Nonlinear' approach can follow the mode
% evolution to larger amplitude. Options are set with i_fit_control.
% SUPPLEMENTAL MATLAB ROUTINES USED
%
% ControlString.m v1.0 Generates diagnostic string
% ScanParameters.m v1.0 Allows variation of parameter for scan
% CalcPertAreaChange.m v1.0 Area diagnostic
% FitLoveEnd.m v1.0 Trims tips of large amplitude seed, see IVC of paper
% dZdt.m v1.0 Calculates velocity, implementation of Eq. 28 of paper
% PiecewiseIntegration.m v1.0 Breaks simulation run into time steps
% CalcLinEvolution.m v1.0 Calculates the predicted linear solution of a perturbed Kirchhoff vortex, see IIC of paper
% LoadExpData.m v1.0 Loads external file of experimental data
% LoveLinearFit.m v1.0 Fits Love modes, minimizes distance^2 along perpendicular direction, see IIE of paper
% LoveNonlinearFit.m v1.0 Fits Love modes, minimizes distance^2 between contour and ellipse
% FitString.m v1.0 Generates string with results of a fit
% ReadAndPlotResults.m v1.0 Sample post-simulation routine to open and plot fit results
%
% EXTERNAL MATLAB ROUTINES USED
%
% ellipsefit.m Ellipse fitter 2005-08-09 D.C. Hanselman
% www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=7012
% curvspace.m Generates evenly spaced points along curve 22 Mar 2005 Yo Fukushima
% www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=7233
% intersections.m Computes intersections where 2 curves meet v1.10 Douglas M. Schwarz
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=11837
% nearestneighbor.m Calculates nearest point 2006 Richard Brown
% www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=12574
% LoadAscii.m Loads an ascii file, external routine of unknown origin
% BASE STATE PARAMETERS
aspect_ratio = 3.5; % convention is a/b: semi-major_axis/semi-minor_axis
semi_major_axis = 1.0; % a=1 is convention used in paper, choice impacts perturbation specifications
zeta = 1.; % Note: neg. vorticity confuses fitting routines b/c CCW ellipse rotation expected
n_ellipse_points = 500; % points along circumference
% PERTURBATION PARAMETERS
m_1 = 3; % This specifies the m-number of the seeded perturbation
m_2 = 4; % Fit control. If m_2=0, will fit to only the m_1 mode
m_3 = 0; % Fit control. If m_3=0, will fit to only m_1 and m_2 modes
% Nonlinear fit only: if m_3<0, allow a and b to vary
n_vary = 1 ; % Can loop to scan over quantities: e.g. perturbation phase, aspect ratio
i_pert_control = 1; % 0 to not apply any perturbation
% 1 to specify in form of Eq. 15b (alpha_m, beta_m)
% 2 to specify in form of Eq. 16 (alpha_m, beta_m, linearized)
% 3 to specify in form of Eq. 15a (A_m, B_m)
% CASE 1 and 2 parameters
alpha_m = 2.e-3; % mode amplitude, 1E-5 value gives good agreement with linear theory
beta_m = 0.0; % specified mode angle in radians
i_beta_control = 0; % 0 to use specified beta_m
% 1 to use angle of growing unstable eigenfunction
% 2 to use angle of decaying unstable eigenfunction
% CASE 3 parameters
A_m = 5.e-3; % mode amplitude
B_m = 0.0; % mode amplitude
i_pert_end_control = 0; % 0 to do nothing
% 1 to trim end of perturbed initial condition (see section IVC of paper)
% ITEGRATOR CONTROL PARAMETERS
i_integrator_control = 1; % 0 to not run - use when integration has been done before and wish to fit
% 1 to run in one time block - most accurate, but will hang if contour too perturbed
% 2 to break run into time steps - useful if hanging occurs
% 3 to create 'fake data' based on linear theory predictions
% 4 to read in 'real data' from experiments, a and b are fitted in this case
% GENERAL PARAMETERS
sim_end_time = 10.0; % simulation run time
frame_steps = 0.5; % gives time steps at which simulation values are written
RelTolValue = 1.E-8; % Integrator tolerance
AbsTolValue = 1.E-12; % Numerical noise is not reduced below these value
% CASE 2 PARAMETERS
n_pre_evolve = 0.0; % Set to non-zero to integrate that amount of time prior to beginning segmented integration
sim_steps = 1.; % simulation time steps for segmented integration
safe_max_del_r_mult = 10.; % when inter-particle spacing increases by this, halt simulation (filamentation)
safe_max_min_r_mult = 1.20; % when min r increases by this, halt (fissioning)
safe_max_max_r = 1.40; % when max r reaches this, halt simulation (filamentation)
safe_min_r = .002; % when particle radius below this, halt (fissioning)
safe_mult_time = 10; % halt if time required for step too large
i_resample = 0; % equalizes spacing between nodes
t_resample_start = 10.0; % time to begin resampling node spacing
t_resample = 2.0; % resample every t_resample times, once begun
% CASE 4 PARAMETERS
exp_scale_factor = 455.; % When fitting experimental data, set to normalize major axis near to 1
exp_head = 'I:\Andor\2006\1024\FZ\FZ_'; % where contours.dat and times.dat files can be read
exp_T_R = 10.9; % experimental value of rotation period
% FITTING CONTROL PARAMETERS
i_fit_control = 1; % 0 to do no fitting
% 1 for linear fitting to 1, 2 or 3 modes
% 2 for nonlinear fitting to 2 modes
i_fit_range_control = 1; % 0 to fit one point at n_fit_param
% 1 to fit everything
% 2 to begin from end minus n_fit_param
% 3 start at beginning, do n_fit_param fits
n_fit_param = 1; % quantity used to specify fit range
i_fit_step = 1; % steps between fits
% OUTPUT OF RESULTS
i_write_fit_output = 1; % 1 to write mode fit output
file_head = '0319\FZ_'; % directory to write fit results to
drive_string = 'I:\Andor\2008\';
% DIAGNOSTICS
i_calc_pert_area_change = 0; % 1 to calculate and display change in area from perturbation
i_plot_final_state = 0; % 1 to see end state after integration
i_plot_perp_matrix = 0; % 1 to see perp matrix, if nonlinear fit instead serves as flag to print int_out_vec
i_plot_fit = 1; % 1 to plot fit at each time
i_plot_amp_phase = 0; % 1 to plot amp and phase versus time for entire run
i_plot_safe_calc = 1; % 1 to plot after each step if doing piecewise integration
% INITIALIZE OUTPUT MATRIX
n_time_points_sim = double(int32(sim_end_time / frame_steps + 1));
if i_integrator_control > 0 % re-initialize if doing new simulation run
if i_integrator_control == 4 % with experimental data, open to determine length
fid = fopen([exp_head,'contours.dat'], 'r');
x = fscanf(fid, '%f ', [3000 inf]);
fclose(fid);
n_time_points_sim = size(x,2) / 2 ;
end
out_matrix = zeros(14,n_vary,n_time_points_sim);
% t, t/T_R, a, b, phi, alpha_1, beta_1, alpha_2, beta_2, alpha_3, beta_3, min_r, max_r, chi out_matrix
% 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14
end
% START OF SIMULATION
warning('off','MATLAB:divideByZero')
ControlString(i_integrator_control, i_pert_control, i_beta_control, i_pert_end_control, i_fit_control, m_2, m_3) % Print run information
s = linspace(0,2*pi,n_ellipse_points+1);
s = s(1:n_ellipse_points);
for i_vary = 1:1:n_vary % FOR LOOP OVER VARIED QUANTITIES
% ******************* CALCULATE UNPERTURBED ELLIPSE PROPERTIES **********************************
dummy = ScanParameters(n_vary); % can use n_vary loop to vary a parameter for a scan, not used
if n_vary == 1
aspect_ratio = aspect_ratio;
else
aspect_ratio = aspect_ratio %+ (n_vary - 1) * 0.1; for example
end
a = semi_major_axis;
b = semi_major_axis / aspect_ratio;
c = (a^2 - b^2)^0.5;
omega = zeta * a*b/(a+b)^2; % Eq. 12, Kirchhoff expression for ellipse rotation frequency
T_R = 2 * pi / omega;
u_0 = acosh(a/c); % Use (u,v) for cylindrical elliptical coordinate, rather than (psi,eta)
h_0_squared = 1. ./ (b^2*cos(s).^2 + a^2.*sin(s).^2); % This quantity is 1 / Jacobian, see Eq. 14, also Love (1893)
gamma = aspect_ratio; % Linear theory predictions, see Love (1893) and Guo et al. (2004)
mu_m_plus = zeta / 2 * ( ((gamma - 1)/(gamma + 1))^m_1 + (2 * m_1 * gamma / (gamma + 1)^2 - 1) ); % Eq. 19a
mu_m_neg = zeta / 2 * ( ((gamma - 1)/(gamma + 1))^m_1 - (2 * m_1 * gamma / (gamma + 1)^2 - 1) ); % Eq. 19b
lambda_m = sqrt( abs(mu_m_plus * mu_m_neg)); % growth rate / oscillation frequency, real, equivalent to Eq. 17
beta_m_plus = (-1) * atan(sqrt(mu_m_neg / mu_m_plus)) / m_1; % growing eigenfunction phase, Eq. 21a
beta_m_neg = atan(sqrt(mu_m_neg / mu_m_plus)) / m_1; % decaying eigenfunction phase, Eq. 21b
% ************** CREATE PERTURBED INITIAL CONDITION, LOAD INTO Z0_START ***************************
switch i_pert_control
case 0 % no pert, but add noise to avoid fitting singularity problem
u_pert = rand(size(s)) * 1.E-12; alpha_m_0 = 0; beta_m_0 = 0;
Z0_Love = c * cosh(u_0 + u_pert) .* cos(s) + i * c * sinh(u_0 + u_pert ) .* sin(s);
case 1 % Eq. 15b form of perturbation
switch i_beta_control
case 0
u_pert = alpha_m * cos(m_1*(s-beta_m)) .* h_0_squared; beta_m_0 = beta_m;
case 1
u_pert = alpha_m * cos(m_1*(s-beta_m_plus)) .* h_0_squared; beta_m_0 = beta_m_plus;
case 2
u_pert = alpha_m * cos(m_1*(s-beta_m_neg)) .* h_0_squared; beta_m_0 = beta_m_neg;
end
alpha_m_0 = alpha_m;
Z0_Love = c * cosh(u_0 + u_pert) .* cos(s) + i * c * sinh(u_0 + u_pert ) .* sin(s);
case 2 % Eq. 16 form of perturbation (linearized)
switch i_beta_control
case 0
Z0_Love = a*cos(s) + b*i*sin(s) + alpha_m*(b*cos(s) + a*i*sin(s)).* cos(m_1*(s-beta_m)) .* h_0_squared; beta_m_0 = beta_m;
case 1
Z0_Love = a*cos(s) + b*i*sin(s) + alpha_m*(b*cos(s) + a*i*sin(s)).* cos(m_1*(s-beta_m_plus)) .* h_0_squared; beta_m_0 = beta_m_plus;
case 2
Z0_Love = a*cos(s) + b*i*sin(s) + alpha_m*(b*cos(s) + a*i*sin(s)).* cos(m_1*(s-beta_m_neg)) .* h_0_squared; beta_m_0 = beta_m_neg;
end
alpha_m_0 = alpha_m;
case 3 % Eq. 15a form of perturbation
u_pert = (A_m * cos(m_1*(s)) + B_m * sin(m_1*(s))) .* h_0_squared;
alpha_m_0 = sqrt(A_m^2 + B_m^2); beta_m_0 = atan(B_m / A_m) / m_1;
Z0_Love = c * cosh(u_0 + u_pert) .* cos(s) + i * c * sinh(u_0 + u_pert ) .* sin(s);
end
Z0_Love = Z0_Love.';
if i_pert_end_control == 0 % if desired, trim end of perturbation as per section IVC of paper
Z0_start = Z0_Love;
else
Z0_Love_trim =FitLoveEnd(Z0_Love, m_1, pi/8, pi/4, 0);
Z0_start = Z0_Love_trim;
end
if i_calc_pert_area_change == 1; del_area = CalcPertAreaChange(a,b,s,Z0_start); end; % print area diagnostic
% ******************************** PERFORM INTEGRATION ***************************
options = odeset('RelTol',RelTolValue,'AbsTol',AbsTolValue);
switch i_integrator_control
case 0 % do no integration
case 1
int_st = 'Single step integration.';
disp(['Begin integration with end time: ', num2str(sim_end_time)]); tic
[T Y_Love] = ode113(@(t,Z) dZdt(Z,s,zeta),[0:frame_steps:sim_end_time],Z0_start,options); toc
case 2 % integrate in steps, determine safe halt time due to filamentation/fusion
if n_pre_evolve ~= 0
disp(['Begin pre_evolution integration with end time: ', num2str(n_pre_evolve)]); tic
[T Y_Love] = ode113(@(t,Z) dZdt(Z,s,zeta),[0:frame_steps:n_pre_evolve],Z0_start,options); toc
Z0_start = Y_Love(end,:);
end
in_vec = [RelTolValue, AbsTolValue, a, b, zeta, m_1, alpha_m_0, beta_m_0, sim_end_time, sim_steps, frame_steps];
in_vec = [in_vec, safe_max_del_r_mult, safe_max_min_r_mult, safe_max_max_r, safe_mult_time, safe_min_r];
in_vec = [in_vec, i_resample, t_resample, t_resample_start, i_plot_safe_calc ];
in_vec = [in_vec, n_ellipse_points, n_time_points_sim, Z0_start'];
out_vec = PiecewiseIntegration(in_vec);
T = out_vec(:,1) + n_pre_evolve;
Y_Love = out_vec(:,2:size(out_vec,2));
case 3 % create 'fake data' using predictions of linear theory
switch i_pert_control
case 0 % no pert, make non-zero to avoid fitting singularity
A_m = 1.E-13; B_m = 1.E-13;
case 1 % convert alpha_m, beta_m to A_m and B_m
switch i_beta_control
case 0
beta_m_0 = beta_m;
case 1
beta_m_0 = beta_m_plus;
case 2
beta_m_0 = beta_m_neg;
end
A_m = sqrt( alpha_m^2 / ( 1 + tan(m_1 * beta_m_0)^2));
B_m = A_m * tan(m_1 * beta_m_0);
case 2
disp('change i_pert control (linearized perturbation not a valid input for linear theory solution)')
case 3 % use the A_m and B_m specified
end
in_vec = [lambda_m, mu_m_plus, mu_m_neg, m_1, sim_end_time, frame_steps, a, b, zeta, n_ellipse_points, A_m, B_m];
out_vec = CalcLinEvolution(in_vec);
T = out_vec(:,1);
Y_Love = out_vec(:,2:size(out_vec,2));
case 4 % load data from experiment
S1.exp_head = exp_head;
S1.exp_scale_factor = exp_scale_factor;
S1.n_ellipse_points = n_ellipse_points;
tic; out_vec = LoadExpData(S1); toc;
T = out_vec(:,1);
Y_Love = out_vec(:,2:size(out_vec,2));
end % i_integrator_control switch
if i_plot_final_state; plot(real(Y_Love(end,:)),imag(Y_Love(end,:))); axis equal; end;
n_time_points_fit = length(T); % may be different than n_time_points_sim
out_matrix(1,i_vary,1:length(T)) = T;
if i_integrator_control ~= 4; % time_vec is time normalized by rotation period, length may be less than n_time_points_sim
out_matrix(2,i_vary,1:length(T)) = T / T_R;
else;
out_matrix(2,i_vary,1:length(T)) = T / exp_T_R; % with experiment, rotation period separately measured
end
% ******************************** FITTING ******************************
if i_fit_control > 0
switch i_fit_range_control % set range of fit
case 0
i_fit_start = n_time_points_fit - n_fit_param;
i_fit_end = i_fit_start;
case 1
i_fit_start = 1;
size_vec = size(Y_Love);
i_fit_end = size_vec(1);
case 2
i_fit_start = n_time_points_fit - n_fit_param;
i_fit_end = n_time_points_fit;
case 3
i_fit_start = 1;
i_fit_end = n_fit_param;
end
for i_fit = i_fit_start:i_fit_step:i_fit_end % loop over fitting of each point
x_Love = real(Y_Love(i_fit,:).'); % load input vector for fit
y_Love = imag(Y_Love(i_fit,:).');
fit_in_vec(1) = b;
fit_in_vec(2) = m_1;
fit_in_vec(3) = m_2;
fit_in_vec(4) = m_3;
fit_in_vec(5) = a;
fit_in_vec(6) = n_ellipse_points;
fit_in_vec(8) = i_plot_perp_matrix;
fit_in_vec(10) = i_plot_fit;
if i_fit == i_fit_start
fit_in_vec(7) = 1; % i_force_perp_matrix_creation = 1
fit_in_vec(9) = 0; % i_fit_analysis_prev_use = 0
else
fit_in_vec(7) = 0;
fit_in_vec(9) = 1;
end
switch i_fit_control % perform selected type of fitting
case 0 % do nothing
case 1 % do linear fit, up to 3 modes
int_out_vec = LoveLinearFit(x_Love,y_Love, fit_in_vec);
case 2 % do nonlinear fit, up to 2 modes
int_out_vec = LoveNonLinearFit(x_Love,y_Love, fit_in_vec);
end
% ******************************** LOAD AND DISPLAY FIT OUTPUT ******************************
% t, t/T_R, a, b, phi, alpha_1, beta_1, alpha_2, beta_2, alpha_3, beta_3, min_r, max_r, chi out_matrix
% 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14
out_matrix(6,i_vary,i_fit) = int_out_vec(1); % load results into out_matrix
out_matrix(7,i_vary,i_fit) = int_out_vec(2);
out_matrix(8,i_vary,i_fit) = int_out_vec(3);
out_matrix(9,i_vary,i_fit) = int_out_vec(4);
out_matrix(5,i_vary,i_fit) = int_out_vec(7);
out_matrix(12,i_vary,i_fit) = int_out_vec(8);
out_matrix(13,i_vary,i_fit) = int_out_vec(9);
out_matrix(14,i_vary,i_fit) = int_out_vec(10);
if m_3 >= 0
out_matrix(10,i_vary,i_fit) = int_out_vec(5);
out_matrix(11,i_vary,i_fit) = int_out_vec(6);
out_matrix(4,i_vary,i_fit) = b;
out_matrix(3,i_vary,i_fit) = a;
else
out_matrix(3,i_vary,i_fit) = int_out_vec(5); % nonlinear fitter, a and b can be fit parameters
out_matrix(4,i_vary,i_fit) = int_out_vec(6);
end
if i_fit == i_fit_start % keep track of ellipse rotations, assume CCW rotation (positive zeta)
n_half_rotations = 0;
else
if out_matrix(5,i_vary,i_fit) + n_half_rotations * pi < out_matrix(5,i_vary,i_fit - 1)
n_half_rotations = n_half_rotations + 1; % increment half rotation counter
end
out_matrix(5,i_vary,i_fit) = out_matrix(5,i_vary,i_fit) + n_half_rotations * pi;
end
FitString(n_vary, i_fit, i_vary, T, out_matrix, m_1, m_2, m_3, i_fit_control, i_integrator_control) % print fit results
end % i_fits loop over fit
end % i_fit_control > 1 end
end % i_vary loop (not usually used)
% ******************************** WRITE FIT RESULTS ******************************
if i_write_fit_output == 1
for i_write = 1: n_vary
if i_write < 10
file_ext = ['_0',num2str(i_write),'.dat'];
else
file_ext = ['_',num2str(i_write),'.dat'];
end
write_vec = [alpha_m_0, beta_m_0, a, b, zeta, m_1, m_2, m_3];
dir_string = [drive_string,file_head];
fileout = [dir_string,'out_vector',file_ext];
save(fileout, 'write_vec', '-ascii', '-tabs');
out_matrix_b = (squeeze(out_matrix(:,i_write,:)))';
dir_string = [drive_string,file_head];
fileout = [dir_string,'out_matrix',file_ext];
save(fileout, 'out_matrix_b', '-ascii', '-tabs');
end
end
