Skip to content

Commit c9cf27b

Browse files
mrodrig6mrodrig6
committed
pruned several features, both codes are now mostly consistent, need to centralize the unpacking of variables into a script and, perhaps, modularize the pre_process and SVBDODE routine
1 parent 32f6a14 commit c9cf27b

File tree

4 files changed

+118
-134
lines changed

4 files changed

+118
-134
lines changed

‎graphics/s_plot_fd.m‎

Whitespace-only changes.

‎graphics/s_plot_spectral.m‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
% display result
2+
if plotresult == 1
3+
subplot(3,1,1);
4+
plot(t,R,'k','LineWidth',2);
5+
hold('on');
6+
ylabel('$R$','Interpreter','Latex','FontSize',12);
7+
box on;
8+
if bubtherm == 1
9+
if medtherm == 0
10+
b = zeros(size(a));
11+
end
12+
subplot(3,1,2);
13+
hold on;
14+
box on;
15+
semilogy(t,abs(a(end,:)),'k-','LineWidth',2);
16+
semilogy(t,abs(b(end,:)),'b-','LineWidth',2);
17+
ylabel('$a_N$, $b_M$, $e_N$','Interpreter','Latex','FontSize',12);
18+
axis([0 t(end) 1e-20 1]);
19+
set(gca, 'YScale', 'log');
20+
if masstrans == 0
21+
leg1 = legend('$a_N$','$b_M$','Location','NorthEast','FontSize',12);
22+
else
23+
leg1 = legend('$a_N$','$b_M$','$e_N$','Location','NorthEast','FontSize',12);
24+
end
25+
set(leg1,'Interpreter','latex');
26+
set(gca,'TickLabelInterpreter','latex','FontSize',16);
27+
set(gcf,'color','w');
28+
end
29+
if spectral == 1
30+
subplot(3,1,3);
31+
hold on;
32+
box on;
33+
semilogy(t,abs(c(end,:)),'k-','LineWidth',2);
34+
semilogy(t,abs(d(end,:)),'b-','LineWidth',2);
35+
xlabel('$t$','Interpreter','Latex','FontSize',12);
36+
ylabel('$c_P$, $d_P$','Interpreter','Latex','FontSize',12);
37+
set(gca, 'YScale', 'log');
38+
axis([0 t(end) 1e-20 1]);
39+
leg1 = legend('$c_P$','$d_P$','Location','NorthEast','FontSize',12);
40+
set(leg1,'Interpreter','latex');
41+
set(gca,'TickLabelInterpreter','latex','FontSize',16);
42+
set(gcf,'color','w');
43+
end
44+
end

‎src/spectral/m_imrv2_finitediff.m‎

Lines changed: 51 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,12 @@
1-
function varargout = m_imrv2_finitediff(varargin)
1+
% file m_imrv2_finitediff.m
2+
% brief contains module m_imrv2_finitediff
23

3-
% Description: This code the IMRv2 code expanding beyond IMR from
4-
% Estrada et al. (2018) JMPS. Additional physics have been added including
5-
% the Keller-Miksis with enthalpy and non-Newtonian viscosity.
4+
% brief This module features a fourth- and sixth-order accurate finite
5+
% difference solver of the PDEs involving thermal transport and
6+
% viscoelasticity to solve Rayleigh-Plesset equations
7+
function varargout = m_imrv2_finitediff(varargin)
68

7-
% Problem Initialization
9+
% problem Initialization
810
[eqns_opts, solve_opts, init_opts, tspan_opts, out_opts, acos_opts,...
911
wave_opts, sigma_opts, thermal_opts, mass_opts]...
1012
= f_call_params(varargin{:});
@@ -16,39 +18,39 @@
1618
masstrans = eqns_opts(7); perturbed = eqns_opts(8);
1719
nl = eqns_opts(9);
1820
if (stress == 4); ptt = 1; else; ptt = 0; end
21+
1922
% solver options
2023
method = solve_opts(1); spectral = solve_opts(2);
2124
divisions = solve_opts(3); Nv = solve_opts(4);
2225
Nt = solve_opts(5); Mt = solve_opts(6);
2326
Lv = solve_opts(7); Lt = solve_opts(8);
27+
2428
% dimensionless initial conditions
2529
Rzero = init_opts(1); Uzero = init_opts(2);
2630
p0star = init_opts(3); P8 = init_opts(4);
2731
T8 = init_opts(5); Pv_star = init_opts(6);
2832
Req = init_opts(7); S0 = init_opts(8);
2933
alphax = init_opts(9);
30-
if perturbed == 1
31-
azero = init_opts(10:10+nl-1);
32-
adot_zero = init_opts(10+nl:10+2*nl-1);
33-
end
34+
3435
% time span options
3536
tspan = tspan_opts;
3637
tfin = tspan(end);
3738
% output options
3839
dimensionalout = out_opts(1); progdisplay = out_opts(2);
3940
plotresult = out_opts(3);
40-
% physical parameters%
41+
42+
% physical parameters
43+
4144
% acoustic parameters
4245
Cstar = acos_opts(1); GAMa = acos_opts(2);
4346
kappa = acos_opts(3); nstate = acos_opts(4);
4447
% dimensionless waveform parameters
4548
om = wave_opts(1); ee = wave_opts(2);
4649
tw = wave_opts(3); dt = wave_opts(4);
4750
mn = wave_opts(5); wave_type = wave_opts(6);
48-
if perturbed == 1
49-
l = wave_opts(7:7+nl-1)';
50-
end
51+
5152
pvarargin = [om,ee,tw,dt,mn,wave_type];
53+
5254
% dimensionless viscoelastic
5355
We = sigma_opts(1); Re8 = sigma_opts(2);
5456
DRe = sigma_opts(3); v_a = sigma_opts(4);
@@ -57,17 +59,19 @@
5759
JdotA = sigma_opts(9); v_lambda_star = sigma_opts(10);
5860
iWe = 1/We;
5961
if Ca==-1; Ca=Inf; end
62+
6063
% dimensionless thermal
6164
Foh = thermal_opts(1); Br = thermal_opts(2);
6265
alpha = thermal_opts(3); beta = thermal_opts(4);
6366
chi = thermal_opts(5); iota = thermal_opts(6);
67+
6468
% dimensionaless mass transfer
6569
Fom = mass_opts(1); C0 = mass_opts(2);
6670
Rv_star = mass_opts(3); Ra_star = mass_opts(4);
6771
L_heat_star = mass_opts(5); mv0 = mass_opts(6);
6872
ma0 = mass_opts(7);
6973

70-
% pre_process code
74+
% pre_process
7175

7276
% creates finite difference matrices
7377
D_Matrix_T_C = f_finite_diff_mat(Nt,1,0);
@@ -296,54 +300,43 @@
296300
dXdt = [U; Udot; pdot; Taudot; Tmdot; Cdot];
297301

298302
end
299-
300-
301-
% display result
302-
if plotresult == 1
303-
subplot(3,1,1);
304-
plot(t,R,'k','LineWidth',2);
305-
hold('on');
306-
ylabel('$R$','Interpreter','Latex','FontSize',12);
307-
box on;
308-
if isreal(R)
309-
% oldb sims for certain values were imaginary
310-
% this spits out imaginary solution (BI will take care of it)
311-
axis([0 t(end) 0 (max(R) + min(R))]);
312-
set(gca,'TickLabelInterpreter','latex','FontSize',16)
313-
set(gcf,'color','w');
314-
end
315-
if bubtherm == 1
316-
if medtherm == 0
317-
b = zeros(size(a));
318-
end
319-
subplot(3,1,2);
320-
hold on;
321-
box on;
322-
semilogy(t,abs(a(end,:)),'k-','LineWidth',2);
323-
semilogy(t,abs(b(end,:)),'b-','LineWidth',2);
324-
ylabel('$a_N$, $b_M$, $e_N$','Interpreter','Latex','FontSize',12);
325-
axis([0 t(end) 1e-20 1]);
326-
set(gca, 'YScale', 'log');
327-
if masstrans == 0
328-
leg1 = legend('$a_N$','$b_M$','Location','NorthEast','FontSize',12);
329-
else
330-
leg1 = legend('$a_N$','$b_M$','$e_N$','Location','NorthEast','FontSize',12);
331-
end
332-
set(leg1,'Interpreter','latex');
333-
set(gca,'TickLabelInterpreter','latex','FontSize',16);
334-
set(gcf,'color','w');
335-
end
336-
end
337303
disp('--- COMPLETED SIMULATION ---');
338304

339305
% functions called by solver
340306

307+
function Tw= TW(Tauw)
308+
%calculates the temperature at the bubble wall as a fuction of \tau
309+
Tw = (A_star -1 + sqrt(1+2*Tauw*A_star)) / A_star;
310+
end
341311

342-
% function Cw= CW(Tw,P)
343-
% % Calculates the concentration at the bubble wall
344-
% %Function of P and temp at the wall
345-
% thetha = Rv_star/Ra_star*(P./(f_pvsat(Tw*T8)/P8) -1);
346-
% Cw = 1./(1+thetha);
347-
% end
312+
function Cw= CW(Tw,P)
313+
% calculates the temperature and concentration at the bubble wall
314+
thetha = Rv_star/Ra_star*(P./(f_pvsat(Tw*T_inf)/P_inf) -1);
315+
Cw = 1./(1+thetha);
316+
end
348317

318+
function Tauw= Boundary(prelim)
319+
% solves temperature boundary conditions at the bubble wall
320+
% create finite diff. coeffs.
321+
% coefficients in terms of forward difference
322+
323+
% second order
324+
coeff = [-3/2 , 2 ,-1/2 ];
325+
Tm_trans = Tm(2:3);
326+
T_trans = flipud(Tau(end-2:end-1));
327+
C_trans = flipud(C(end-2:end-1));
328+
329+
% sixth order
330+
% coeff= [-49/20 ,6 ,-15/2 ,20/3 ,-15/4 ,6/5 ,-1/6]; %Sixth order coeff
331+
% Tm_trans = Tm(2:7);
332+
% T_trans = flipud(Tau(end-6:end-1));
333+
% C_trans = flipud(C(end-6:end-1));
334+
335+
Tauw =chi*(2*Km_star/L*(coeff*[TW(prelim); Tm_trans] )/deltaYm) +...
336+
chi*(-coeff*[prelim ;T_trans] )/deltaY + Cgrad*...
337+
fom*L_heat_star*P*( (CW(TW(prelim),P)*(Rv_star-Ra_star)+Ra_star))^-1 *...
338+
(TW(prelim) * (1-CW(TW(prelim),P)) ).^(-1).*...
339+
(-coeff*[CW(TW(prelim),P); C_trans] )/deltaY;
349340
end
341+
342+
end

‎src/spectral/m_imrv2_spectral.m‎

Lines changed: 23 additions & 76 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,10 @@
1-
function varargout = m_imrv2_spectral(varargin)
1+
% file m_imrv2_spectral.m
2+
% brief contains module m_imrv2_spectral
23

3-
% Description: This code the IMRv2 code expanding beyond IMR from
4-
% Estrada et al. (2018) JMPS. Additional physics have been added including
5-
% the Keller-Miksis with enthalpy and non-Newtonian viscosity.
4+
% brief This module features a Chebyshev spectral collocation solver of the
5+
% PDEs involving thermal transport and viscoelasticity to solve
6+
% Rayleigh-Plesset equations
7+
function varargout = m_imrv2_spectral(varargin)
68

79
% Problem Initialization
810
[eqns_opts, solve_opts, init_opts, tspan_opts, out_opts, acos_opts,...
@@ -16,38 +18,37 @@
1618
masstrans = eqns_opts(7); perturbed = eqns_opts(8);
1719
nl = eqns_opts(9);
1820
if (stress == 4); ptt = 1; else; ptt = 0; end
21+
1922
% solver options
2023
method = solve_opts(1); spectral = solve_opts(2);
2124
divisions = solve_opts(3); Nv = solve_opts(4);
2225
Nt = solve_opts(5); Mt = solve_opts(6);
2326
Lv = solve_opts(7); Lt = solve_opts(8);
27+
2428
% dimensionless initial conditions
2529
Rzero = init_opts(1); Uzero = init_opts(2);
2630
p0star = init_opts(3); P8 = init_opts(4);
2731
T8 = init_opts(5); Pv_star = init_opts(6);
2832
Req = init_opts(7); S0 = init_opts(8);
2933
alphax = init_opts(9);
30-
if perturbed == 1
31-
azero = init_opts(10:10+nl-1);
32-
adot_zero = init_opts(10+nl:10+2*nl-1);
33-
end
34+
3435
% time span options
3536
tspan = tspan_opts;
3637
tfin = tspan(end);
3738
% output options
3839
dimensionalout = out_opts(1); progdisplay = out_opts(2);
3940
plotresult = out_opts(3);
40-
% physical parameters%
41+
42+
% physical parameters
43+
4144
% acoustic parameters
4245
Cstar = acos_opts(1); GAMa = acos_opts(2);
4346
kappa = acos_opts(3); nstate = acos_opts(4);
4447
% dimensionless waveform parameters
4548
om = wave_opts(1); ee = wave_opts(2);
4649
tw = wave_opts(3); dt = wave_opts(4);
4750
mn = wave_opts(5); wave_type = wave_opts(6);
48-
if perturbed == 1
49-
l = wave_opts(7:7+nl-1)';
50-
end
51+
5152
pvarargin = [om,ee,tw,dt,mn,wave_type];
5253
% dimensionless viscoelastic
5354
We = sigma_opts(1); Re8 = sigma_opts(2);
@@ -170,10 +171,6 @@
170171
end
171172
Z1 = X(:,ic);
172173
Z2 = X(:,id);
173-
if perturbed == 1
174-
aa = X(:,7:7+nl-1);
175-
adot = X(:,7+nl:7+2*nl-1);
176-
end
177174
a = X(:,ia)';
178175
b = X(:,ib)';
179176
c = X(:,ic)';
@@ -317,11 +314,14 @@
317314
J = 0;
318315
JdotX = 0;
319316

320-
% Kelvin-Voigt with neo-Hookean elasticity
317+
% compute stress integral
321318
elseif stress == 1
322-
% compute stress integral
323-
% J = (4*(Req/R) + (Req/R)^4 - 5)/(2*Ca) - 4/Re8*U/R;
324-
% JdotX = -2*U*(Req*(1/R)^2 + Req^4*(1/R)^5)/Ca + 4/Re8*U^2/R^2;
319+
% Kelvin-Voigt with neo-Hookean elasticity
320+
J = (4*(Req/R) + (Req/R)^4 - 5)/(2*Ca) - 4/Re8*U/R;
321+
JdotX = -2*U*(Req*(1/R)^2 + Req^4*(1/R)^5)/Ca + 4/Re8*U^2/R^2;
322+
323+
elseif stress == 2
324+
% quadratic Kelvin-Voigt with neo-Hookean elasticity
325325
J = (3*alphax-1)*(5 - (Req/R)^4 - 4*(Req/R))/(2*Ca) - 4/Re8*U/R + ...
326326
(2*alphax/Ca)*(27/40 + (1/8)*(Req/R)^8 + (1/5)*(Req/R)^5 + (1/2)*(Req/R)^2 - ...
327327
2*R/Req);
@@ -352,7 +352,7 @@
352352
JdotX = 2*sum(cdd.*(Z1dot - Z2dot));
353353

354354
% linear Maxwell, linear Jeffreys, linear Zener
355-
elseif stress == 2
355+
elseif stress == 3
356356
% extract
357357
Z1 = X(ic);
358358
J = Z1/R^3 - 4*LAM/Re8*U/R;
@@ -372,7 +372,7 @@
372372
JdotX = Z1dot/R^3 - 3*U/R^4*Z1 + 4*LAM/Re8*U^2/R^2;
373373

374374
% upper-convected Maxwell, OldRoyd-B
375-
elseif stress == 3
375+
elseif stress == 4
376376
% extract stress sub-integrals
377377
Z1 = X(ic); Z2 = X(id);
378378
% compute new derivatives
@@ -398,63 +398,10 @@
398398
dXdt = [U; Udot; pdot; qdot; Z1dot; Z2dot; Jdot];
399399

400400
end
401-
402-
403-
% display result
404-
if plotresult == 1
405-
subplot(3,1,1);
406-
plot(t,R,'k','LineWidth',2);
407-
hold('on');
408-
ylabel('$R$','Interpreter','Latex','FontSize',12);
409-
box on;
410-
if isreal(R)
411-
% oldb sims for certain values were imaginary
412-
% this spits out imaginary solution (BI will take care of it)
413-
axis([0 t(end) 0 (max(R) + min(R))]);
414-
set(gca,'TickLabelInterpreter','latex','FontSize',16)
415-
set(gcf,'color','w');
416-
end
417-
if bubtherm == 1
418-
if medtherm == 0
419-
b = zeros(size(a));
420-
end
421-
subplot(3,1,2);
422-
hold on;
423-
box on;
424-
semilogy(t,abs(a(end,:)),'k-','LineWidth',2);
425-
semilogy(t,abs(b(end,:)),'b-','LineWidth',2);
426-
ylabel('$a_N$, $b_M$, $e_N$','Interpreter','Latex','FontSize',12);
427-
axis([0 t(end) 1e-20 1]);
428-
set(gca, 'YScale', 'log');
429-
if masstrans == 0
430-
leg1 = legend('$a_N$','$b_M$','Location','NorthEast','FontSize',12);
431-
else
432-
leg1 = legend('$a_N$','$b_M$','$e_N$','Location','NorthEast','FontSize',12);
433-
end
434-
set(leg1,'Interpreter','latex');
435-
set(gca,'TickLabelInterpreter','latex','FontSize',16);
436-
set(gcf,'color','w');
437-
end
438-
if spectral == 1
439-
subplot(3,1,3);
440-
hold on;
441-
box on;
442-
semilogy(t,abs(c(end,:)),'k-','LineWidth',2);
443-
semilogy(t,abs(d(end,:)),'b-','LineWidth',2);
444-
xlabel('$t$','Interpreter','Latex','FontSize',12);
445-
ylabel('$c_P$, $d_P$','Interpreter','Latex','FontSize',12);
446-
set(gca, 'YScale', 'log');
447-
axis([0 t(end) 1e-20 1]);
448-
leg1 = legend('$c_P$','$d_P$','Location','NorthEast','FontSize',12);
449-
set(leg1,'Interpreter','latex');
450-
set(gca,'TickLabelInterpreter','latex','FontSize',16);
451-
set(gcf,'color','w');
452-
end
453-
end
454401
disp('--- COMPLETED SIMULATION ---');
455402

456-
457403
% functions called by solver
404+
458405
% stress differentiator
459406
function [trr,dtrr,t00,dt00] = stressdiff(c,d)
460407
if Nv < 650

0 commit comments

Comments
 (0)