--- /dev/null
+% close all;
+% clear all;
+% n=50; %800; % >=4
+
+% boundary interval
+dr = 1/(n-1);
+r = linspace(0,1,n);
+% lambda = 5.5;
+s = (1 - exp(-lambda*r))/(1 - exp(-lambda));
+
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% ORIGINAL ATTEMPT AT AN INITIAL CONDITION
+%
+% % % H0 = 0.1; % dimensionless curly H (see 8-14-19 notes, p.9 ... can relate to volume if we want)
+% % % theta = 0.0;
+% % % hIn = (H0 + tan(theta))*(1-s(1:n-1).^2) + ...
+% % % (s(1:n-1)-1)*tan(theta);
+% % % yIn = [hIn 1]; %IC for h and IC for xN
+% % % tInit = 0;
+% % % tFinal = 100000.0; %100000000.0;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% SIMILARITY INITIAL CONDITION
+%
+theta = 0;
+tInit = 1;
+% tFinal = 10.0;
+%
+% match initial condition length with length of s(1:n-1)
+%
+E = 5;
+E1 = 1;
+E2 = 1;
+E3 = 2;
+E4 = 3;
+qbarUA = 1; %(2*H0/3); % initial volume for unit width channel (see integral of
+ % initial condition for the case theta = 0 which gives
+ % qbarUA = 2 * H0 / 3
+B_TH = beta(E1/E3,1 + E2/E4);
+eta_N = (E3 * ( (E*E3)/(E2*E4) )^(E2/E4) /B_TH )^(E4/E);
+eta_INIT = s*eta_N;
+phi_INIT = ( (E2*E4)/(E*E3) * (eta_N^E3 - eta_INIT(1:n-1).^E3) ).^(1/E4);
+t_TH_INIT = tInit;
+x_N_TH_INIT = eta_N * ( (1/3) * qbarUA^E4 * t_TH_INIT).^(1/E);
+x_TH_INIT_plot = s*x_N_TH_INIT;
+h_TH_INIT = (qbarUA^2 * 3./t_TH_INIT ).^(1/E) * phi_INIT;
+
+yIn = [h_TH_INIT x_N_TH_INIT]; %IC for h and IC for xN
+
+% plot the initial condition to help debug ... (DMA, 11-8-2019)
+%%figure(88);plot(x_TH_INIT_plot,[h_TH_INIT 0]);hold on;
+%%return
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% time domain, ode23 config
+tspan = [tInit tFinal];
+RelTolVal=10^(-3);
+AbsTolVal=10^(-3);
+options = odeset ('RelTol',RelTolVal,'AbsTol',AbsTolVal,'InitialStep',1e-3); %,'MaxStep',1e-8);
+
+[t,y] = ode23 (@gc_rhsND_nonuniform_s, tspan, yIn', options, n, r, dr, lambda, bta, theta, a);
+
+
+%%%%
+%%%% plotting
+%%%%
+% figure(1);plot(t,y(:,end))
+% xlabel('time','FontSize',16);
+% ylabel('xN','FontSize',16);
+% %
+% figure(4);loglog(t,y(:,end),'*');hold on;
+% xlabel('time','FontSize',16);
+% ylabel('xN','FontSize',16);
+% %
+% % estimate for power law xN = A t^B
+% %
+% B_est = log(y(end,end)/y(end-1,end))/log(t(end)/t(end-1));
+% Ashift = 1;
+% xNpowerlaw = Ashift * t.^B_est;
+% loglog(t,xNpowerlaw,'r--');
+% %
+% %
+% figure(2);
+% plot(s,[y(1,1:n-1) 0]);hold on;
+% %plot(s,[y(100,1:n-1) 0]);
+% plot(s,[y(end,1:n-1) 0]);
+% xlabel('s','FontSize',16);
+% ylabel('h','FontSize',16);
+% %
+%
+%
+% figure(3);
+% plot(s*y(1,end),[y(1,1:n-1) 0],'r--','LineWidth',3);hold on; % initial condition
+% ncurves = 6;
+% for jj=1:ncurves;
+% curve_time_index = floor(jj/(ncurves+1) * length(t));
+% plot(s*y(curve_time_index,end),[y(curve_time_index,1:n-1) 0],'b-.','LineWidth',1);hold on;
+% end
+% plot(s*y(end,end),[y(end,1:n-1) 0],'k-','LineWidth',3);
+% xlabel('x','FontSize',16);
+% ylabel('h','FontSize',16);
+%
+%
+% %%%%% Takagi/Huppert 2008 similarity solution
+% %%%%% constant volume (a=0)
+% %%%%% horizontal flow (H=1)
+% %
+%
+% %
+% num_eta = n;
+% eta = linspace(0,eta_N,num_eta);
+% phi = ( (E2*E4)/(E*E3) * (eta_N^E3 - eta.^E3) ).^(1/E4);
+% %
+% % use the same choices for time as in Figure 3 to plot this similarity solution
+% %
+% figure(3)
+% ncurves = 6;
+% for jj=1:ncurves;
+% curve_time_index = floor(jj/(ncurves+1) * length(t));
+% t_TH_plot = t(curve_time_index);
+% x_N_TH = eta_N * ( (1/3) * qbarUA^E4 * t_TH_plot).^(1/E);
+% x_TH = linspace(0,x_N_TH,num_eta);
+% h_TH = (qbarUA^2 * 3./t_TH_plot ).^(1/E) * phi;
+% %% plot(eta,h_TH,'c--','LineWidth',3);hold on;
+% plot(x_TH,h_TH,'c--','LineWidth',3);hold on;
+% end
+% t_TH_f = t(end);
+% xN_TH_f = eta_N* ( (1/3)*qbarUA^E4 * t_TH_f).^(1/E);
+% x_TH_f = linspace(0,xN_TH_f,num_eta);
+% h_TH_f = (qbarUA^2 * 3./t_TH_f).^(1/E) * phi;
+% plot(x_TH_f, h_TH_f, 'm-', 'LineWidth',3);
+%
+%
+% figure(101)
+% ncurves = 6;
+% for jj=1:ncurves;
+% curve_time_index = floor(jj/(ncurves+1) * length(t));
+% t_TH_plot = t(curve_time_index);
+% x_N_TH = eta_N * ( (1/3) * qbarUA^E4 * t_TH_plot).^(1/E);
+% loglog(t_TH_plot,x_N_TH,'c*');hold on;
+% end
\ No newline at end of file
--- /dev/null
+function rhs = gc_rhsND_nonuniform_s(t, yIn, n, r, dr, lambda, bta, theta, a);
+%
+% define constants
+%a=0;
+q=1; % make this equal to qbarUA!
+%W=1;
+%rho=1;
+%nu=1;
+%g=1;
+%%% theta=0.1; %pi/6; % defined in gc_mol.m
+b0=bta;
+b1=bta;
+b2=bta; % dimensionless betas (CURRENTLY USING FOR PARAMETER TESTING)
+
+% y = [h_0, h(s=0), ... h(s=1-ds), h(s=1), xN]
+y = [0 yIn(1:n-1)' 0 yIn(n)]; % h0=y(1) will be set shortly; h_n = y(n+1) = 0 (BC)
+h = y(1:n+1); % h0, h1, ..., hn
+
+% construct s from r; compute drds
+% using g2 from Dalwadi et. al
+% r = g2(s; lambda) so r->g2, s->eta
+s = (1 - exp(-lambda*r))/(1 - exp(-lambda));
+drds = (1 - exp(-lambda))./(lambda * (1 - s*(1 - exp(-lambda))));
+
+% define terms for the rhs equations,
+f = 1/3*h.^3 + b0*h.^2 + (b1+b2)*h;
+if a==0
+ y(1) = y(3);
+else
+ y(1) = y(3) + 2*dr*y(end) * (a*q*t^(a-1)/(cos(theta)*f(2)) - tan(theta))/drds(1);
+end
+
+h_j = y(2:n); % h1, h2, ..., h_(n-1)
+h_jp = y(3:n+1); % h2, h3, ..., hn
+h_jm = y(1:n-1); % h0, h1, ..., h_(n-2)
+
+f(1) = 1/3*y(1).^3 + b0*y(1).^2 + (b1+b2)*y(1);
+
+
+dhdr = (h_jp - h_jm)/(2*dr); % defined on h_1,...,h_(n-1)
+dhdr_xN = (3*y(n+1) - 4*y(n) + y(n-1))/(2*dr); % defined on h_n; for dxNdt eq
+
+fp = (f(3:n+1) + f(2:n))/2; % don't need fp(h_0)
+fm = (f(2:n) + f(1:n-1))/2; % don't need fm(h_n)
+
+drds_p = (drds(2:n) + drds(1:n-1))/2;
+drds_m = (drds(1:n-1) + [drds(1) drds(1:n-2)])/2;
+
+Gp = fp .* (drds_p/y(end).*(h_jp - h_j)/dr - tan(theta)); % dont need Gp(h_0)
+Gm = fm .* (drds_m/y(end).*(h_j - h_jm)/dr - tan(theta)); % dont need Gm(h_n)
+dGdr = (Gp - Gm)/dr;
+
+% define rhs
+% this is for all n-1 dhdt and one dxNdt eq'ns
+% indexing changes here; idx=1:n-1 is for h_1 to h_(n-1), idx=n is for xN
+DY(n) = - cos(theta)*(b1 + b2)*(drds(end)*dhdr_xN/y(end) - tan(theta));
+DY(1:n-1) = s(1:n-1)/y(end)*DY(n).*dhdr.*drds(1:n-1) + cos(theta)*(dGdr/y(end).*drds(1:n-1));
+rhs = DY';
\ No newline at end of file
--- /dev/null
+n = 50; %linspace(25,100,4);
+lambdavec = linspace(7,8,2);
+bta = 1e-4; %10.^linspace(-6,-4,3)
+tFinalvec = 10.^linspace(1,3,2);
+avec = [0.1, 0.5, 1, 2, 10];
+% increase domain of lambda, and startng from large betas, incorporate a
+% beta vector
+
+% the following will be length(nvec) x length(lambdavec) x length(betavec)
+xN = [];
+err = [];
+runTimes = [];
+params = [];
+total_runs = num2str(length(avec)*length(lambdavec)*length(tFinalvec));
+done = 0;
+for i=1:length(avec)
+ a = avec(i);
+ for j=1:length(lambdavec)
+ lambda = lambdavec(j);
+ for k=1:length(tFinalvec)
+ tFinal = tFinalvec(k);
+ tic
+ gc_molND_nonuniform_s; % plotting is off for these purposes
+ runTime = toc;
+ TH_sol;
+ UA_sol = [y(end,1:n-1) 0];
+ astr = num2str(a);
+ lambdastr=num2str(lambda);
+ tfstr=num2str(tFinal);
+ display(['Solving a=',astr,', lambda=',lambdastr,', tFinal=',tfstr]);
+ sol_err = norm(h_TH_f-UA_sol,1)/norm(h_TH_f,1);
+
+ xN(i,j,k) = y(end,end);
+ err(i,j,k) = sol_err;
+ runTimes(i,j,k) = runTime;
+ done = done + 1;
+ done_str = num2str(done);
+ display(['Completed ', done_str, ' out of ', total_runs]);
+% dlmwrite('xN.txt', xN);
+ %type xN.txt;
+% dlmwrite('err.txt.', err);
+ %type err.txt;
+% dlmwrite('runTimes.txt', runTimes);
+ %type runTimes.txt;
+ save('tester_out_05_22.mat', 'xN', 'err', 'runTimes','params');
+ end
+ end
+end
+%save('xN.m', '-ascii', 'xN')
projects / viscous-gravity-currents / .git / commitdiff