From: ibidyouadu Date: Fri, 5 Jun 2020 15:26:59 +0000 (-0400) Subject: Initial commit X-Git-Url: http://git.angelumana.com/?a=commitdiff_plain;h=271b0a028db4ba1920dd012cd09073e14a2181f2;p=viscous-gravity-currents%2F.git Initial commit --- 271b0a028db4ba1920dd012cd09073e14a2181f2 diff --git a/README.md b/README.md new file mode 100644 index 0000000..213b098 --- /dev/null +++ b/README.md @@ -0,0 +1 @@ +"# viscous-gravity-currents" diff --git a/documentation/gravity_current_umana_anderson.pdf b/documentation/gravity_current_umana_anderson.pdf new file mode 100644 index 0000000..5453915 Binary files /dev/null and b/documentation/gravity_current_umana_anderson.pdf differ diff --git a/model_code/TH_sol.m b/model_code/TH_sol.m new file mode 100644 index 0000000..5615f97 --- /dev/null +++ b/model_code/TH_sol.m @@ -0,0 +1,19 @@ +% the last ~40 lines of gc_molND_nonuniform_s.m +H0 = 0.1; +E = 5; +E1 = 1; +E2 = 1; +E3 = 2; +E4 = 3; +qbarUA = 1; % 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); +num_eta = n; +eta = s*eta_N; +phi = ( (E2*E4)/(E*E3) * (eta_N^E3 - eta.^E3) ).^(1/E4); + +t_TH_f = t(end); +xN_TH_f = eta_N* ( (1/3)*qbarUA^E4 * t_TH_f).^(1/E); +h_TH_f = (qbarUA^2 * 3./t_TH_f).^(1/E) * phi; \ No newline at end of file diff --git a/model_code/gc_molND_nonuniform_s.m b/model_code/gc_molND_nonuniform_s.m new file mode 100644 index 0000000..5ac0632 --- /dev/null +++ b/model_code/gc_molND_nonuniform_s.m @@ -0,0 +1,143 @@ +% 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 diff --git a/model_code/gc_rhsND_nonuniform_s.m b/model_code/gc_rhsND_nonuniform_s.m new file mode 100644 index 0000000..623dc20 --- /dev/null +++ b/model_code/gc_rhsND_nonuniform_s.m @@ -0,0 +1,58 @@ +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 diff --git a/model_code/tester.m b/model_code/tester.m new file mode 100644 index 0000000..68366ad --- /dev/null +++ b/model_code/tester.m @@ -0,0 +1,49 @@ +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')