% MATLAB script to calculate Ra_cr using aspect.
clear;
close all;

% all user-specified stuff goes here:
addpath(fullfile('../Dropbox/Curtis_Williams/steinberger/export_fig'))
base_input_file = 'convection-box-base-nd.prm';
path_to_aspect = '/opt/aspect/build/aspect';              % normal aspect
%path_to_aspect = '/opt/aspect_ba/aspect/build/aspect';      % Boussinesq aspect
initial_global_refinement = 4;

aspect_min = pi/3;  % minimum aspect ratio
aspect_max = 5*pi;  % maximum aspect ratio
naspect = 10;        % number of aspect ratios

% material properties
thermal_conductivity = 2.2;
reference_density = 3300;
heat_capacity = 1250;
alpha = 3e-5;
Tbot = 25;
g=10;
p=.01;                        %perturbation amplitude
b = 3e6; % layer depth

% these parameters must be consistent with defs in input file
kappa = thermal_conductivity/(reference_density*heat_capacity);
%aspect_list = logspace(log10(aspect_min),log10(aspect_max),naspect);
aspect_list = 1./(linspace(1/aspect_max,1/aspect_min,naspect));


figure(1);
hold on;
as = linspace(aspect_min/2,aspect_max,100);
lams = as*b; % lambda is width of box
racr1 = (pi^2 + 4*pi^2*b^2./(lams).^2).^3./(4*pi^2*b^2./(lams).^2); % analytic solution from T&S

figure(1);
plot(2*pi*b./lams,racr1,'k--');
xlabel('Aspect ratio (2{\pi}b/\lambda)');
ylabel('Ra');
set(gca,'YLim',[0 6000]);
set(gca,'XLim',[0 7]);

for ias=1:length(aspect_list)
    aspect_ratio = aspect_list(ias);
    lam = aspect_ratio*b;
    %% this section defines a dictionary of strings to replace in the 'base' input file
    dictionary = {'WIDTH',num2str(lam);
        'HEIGHT',num2str(b);
        'IREF',num2str(initial_global_refinement)
        'PMAG',num2str(p);
        'TBOT',num2str(Tbot);
        'RHO',num2str(reference_density);
        'CP',num2str(heat_capacity);
        'KT',num2str(thermal_conductivity);
        'ALPHA',num2str(alpha);
        'GRAV',num2str(g);
        };
    
    %
    racr_high = 1e5;
    racr_low  = 250;
    tol = 1;
    
    while( racr_high-racr_low > tol ) % outer loop for bisection algorithm
        this_ra = 10^((log10(racr_high)+log10(racr_low))/2);
        % generate new input file
        unix('rm -f input.prm');
        unix(['cp ' base_input_file ' input.prm']);
        
        this_eta = reference_density*g*alpha*Tbot*(b)^3/this_ra/kappa;
        
        % update dictionary
        this_dictionary = dictionary;
        this_dictionary(end+1,1:2) = {'VISCOSITY',sprintf('%e',this_eta)};
        this_dictionary(end+1,1:2) = {'END_TIME',sprintf('%e',1e10)};
%         this_dictionary(end+1,1:2) = {'DTMAX',sprintf('%e',dtmax)};
        
        % use dictionary to replace new input file using keys using 'sed'
        % command
        nkeys = size(this_dictionary,1);
        for i=1:nkeys
            unix( ['sed -i ''s/' this_dictionary{i,1} '/' this_dictionary{i,2} '/g'' input.prm'] );
        end
        
        % run aspect
        unix(['rm -rf output']);
        unix(['unset DYLD_LIBRARY_PATH; unset LD_LIBRARY_PATH; source ~/.bashrc; mpirun -n 1 ' path_to_aspect ' input.prm']);
        
        % read statistics file
        [t,vrms] = read_statistics_ba('output/statistics');

        dvdt = (vrms(end)-vrms(1))/(t(end)-t(1));
        
        if( dvdt > 0 )
            racr_high = this_ra;
            figure(1);
            plot(2*pi*b/lam,this_ra,'ro');
        else
            racr_low = this_ra;
            figure(1);
            plot(2*pi*b/lam,this_ra,'b.');
        end
        drawnow;
    end
    
    racr = (pi^2 + 4*pi^2*b^2/lam^2)^3/(4*pi^2*b^2/lam^2); % formula from T&S
    
    figure(1);
    plot(2*pi*b/lam,this_ra,'ro','MarkerFaceColor','r');
    
    err = norm(this_ra-racr)/norm(racr);
    errors(ias) = err;
    results(1,ias) = 2*pi*b/lam;
    results(2,ias) = err;
    
    disp(['L2 norm of error ' num2str(err)]);
end


print('racr_plot.eps','-depsc');

results
