debugged, in part, the bubtherm, need the Rayleigh solution to work

Author: mrodrig6

src/f_call_params.m

@@ -33,7 +33,8 @@
 % remaining parameters
 
 % viscosity variables
-[mu8,Dmu,v_a,v_nc,v_lambda] = f_nonNewtonian_Re(vmaterial); % non-Newtonian viscosity
+[mu8,Dmu,v_a,v_nc,v_lambda] = ...
+    f_nonNewtonian_Re(vmaterial); 
 
 % pressure variables
 Pv              = f_pvsat(T8);  
@@ -61,7 +62,6 @@
         case 'eps3',        eps3 = varargin{n+1};
         case 'vapor',       vapor = varargin{n+1};
         case 'masstrans',   masstrans = varargin{n+1};
-        case 'perturbed',   perturbed = varargin{n+1};
 
         % solver options
         case 'method',  method = varargin{n+1};
@@ -81,14 +81,10 @@
         case 'u0',      U0 = varargin{n+1};
         case 'req',     Req = varargin{n+1};
                 P0 = (P8 + 2*S/Req - Pv*vapor)*((Req/R0)^(3));
-        case 'a0',      a0 = varargin{n+1};
-        case 'adot0',   adot0 = varargin{n+1}; 
-        case 's0',      S0 = varargin{n+1};
 
         % output options
         case 'dimout',  dimensionalout = varargin{n+1};
         case 'pdisp',   progdisplay = varargin{n+1};
-        case 'plot',    plotresult = varargin{n+1};
 
         % acoustic options
         case 'rho8',    rho8 = varargin{n+1};
@@ -111,7 +107,7 @@
         case 'g',       G = varargin{n+1};
         case 'lambda1', lambda1 = varargin{n+1};
         case 'lambda2', lambda2 = varargin{n+1};
-        case 'alphax', alphax = varargin{n+1};
+        case 'alphax',  alphax = varargin{n+1};
         case 'surft',    S = varargin{n+1};
                 P0 = (P8 + 2*S/Req - Pv*vapor)*((Req/R0)^(3));
         case 'vmaterial', vmaterial = varargin{n+1}; 
@@ -151,10 +147,6 @@
     C8 = sqrt(nstate*(P8 + GAM)/rho8); 
 end
 
-if ~medtherm && Mt >= 0
-    Mt = -1;
-    disp('RESETTING Mt = -1 for medtherm == 0');
-end
 check = 1-isnumeric(radial);
 if check || radial > 4 || radial <= 0 
     error('INPUT ERROR: radial must be 1, 2, 3, or 4');
@@ -366,7 +358,7 @@
 % solver options
 solve_opts = [method spectral divisions Nv Nt Mt Lv Lt];
 % dimensionless initial conditions
-init_opts = [Rzero Uzero pzero P8 T8 Pv_star Req_zero S0 alphax];
+init_opts = [Rzero Uzero pzero P8 T8 Pv_star Req_zero alphax];
 % time span options
 tspan_opts = tvector;
 % output options

src/f_odesolve.m

@@ -11,14 +11,14 @@
     if divisions == 0
         options = odeset();
     else
-        options = odeset('MaxStep',tfin/divisions,'RelTol',1e-8,'AbsTol',1e-8);
+        options = odeset('MaxStep',tfin/divisions);%'RelTol',1e-8,'AbsTol',1e-8);
     end
     [t,X] = ode23tb(bubble,tspan,init,options);
 elseif method == 45
     if divisions == 0
         options = odeset('NonNegative',1,'AbsTol',1e-8,'RelTol',1e-8);
     else
-        options = odeset('NonNegative',1,'MaxStep',tfin/divisions,'RelTol',1e-8);
+        options = odeset('NonNegative',1,'MaxStep',tfin/divisions);%'RelTol',1e-8);
     end
     [t,X] = ode45(bubble,tspan,init,options);
 else

src/m_imrv2_finitediff.m

@@ -28,8 +28,7 @@
 Rzero           = init_opts(1);  Uzero           = init_opts(2); 
 p0star          = init_opts(3);  P8              = init_opts(4); 
 T8              = init_opts(5);  Pv_star         = init_opts(6); 
-Req             = init_opts(7);  S0              = init_opts(8);
-alphax          = init_opts(9);
+Req             = init_opts(7);  alphax          = init_opts(8);
 
 % time span options
 tspan = tspan_opts;
@@ -104,9 +103,12 @@
 T = zeros(-1,1);
 if bubtherm
     Tau0 = zeros(Nt,1);
-    Tm0 = ones(Mt ~= -1);
 else
     Tau0 = zeros(-1,1);
+end
+if medtherm
+    Tm0 = ones(Mt ~= -1);
+else
     Tm0 = zeros(-1,1);
 end
 if masstrans
@@ -239,8 +241,8 @@
 
         % temperature inside the bubble
         U_vel = (chi/R*(kappa-1)*DTau-y*R*pdot/3)/(kappa*p);
-        first_term = (DDTau*chi/R^2+pdot)*( K_star*T/p*kapover);
-        second_term = -DTau*((1/R)*(U_vel-y*U));
+        first_term = (DDTau*chi./R^2+pdot)*(K_star*T./p*kapover);
+        second_term = -DTau.*((1./R).*(U_vel-y.*U));
 
         Taudot= first_term+second_term; 
         Taudot(end) = 0;
@@ -262,9 +264,9 @@
         % internal pressure equation
         pdot = 3/R*(chi*(kappa-1)*DTau(end)/R-kappa*p*U);
         % temperature inside the bubble
-        U_vel = (chi/R*(kappa-1).*DTau-y*R*pdot/3)/(kappa*p);
-        first_term = (DDTau.*chi./R^2+pdot).*( K_star.*T/p*kapover );
-        second_term = -DTau.*(1/(R).*(U_vel-y*U));
+        U_vel = (chi/R*(kappa-1)*DTau-y*R*pdot/3)/(kappa*p);
+        first_term = (DDTau.*chi./R^2+pdot).*(K_star.*T./p*kapover);
+        second_term = -DTau.*((1/R).*(U_vel-y*U));
 
         Taudot= first_term+second_term; 
         Taudot(end) = 0;

src/m_imrv2_spectral.m

@@ -28,8 +28,7 @@
 Rzero           = init_opts(1);  Uzero           = init_opts(2); 
 p0star          = init_opts(3);  P8              = init_opts(4); 
 T8              = init_opts(5);  Pv_star         = init_opts(6); 
-Req             = init_opts(7);  S0              = init_opts(8);
-alphax          = init_opts(9);
+Req             = init_opts(7);  alphax          = init_opts(8);
 
 % time span options
 tspan = tspan_opts;

toolchain/f_generate_test.m

@@ -21,20 +21,20 @@
     end
 end
 
-testb = testb + 1;
-teste = teste + 1;
 % bubtherm tests
-for id = testb:teste
-    filename = strcat('../unit_tests/',ids{id},'.dat'); 
-    [ts,Rs,Us] = m_imrv2_spectral('bubtherm',1);
-    [t,R,U] = m_imrv2_finitediff('bubtherm',1);
+tvector = linspace(0,20E-6,100);
+for id = testb+3
+    filename = strcat('../unit_tests/',ids{id+4},'.dat'); 
+    [ts,Rs,Us] = m_imrv2_spectral('radial',id,'bubtherm',1,'tvector',tvector);
+    [tf,Rf,Uf] = m_imrv2_finitediff('radial',id,'bubtherm',1,'Nt',140,'tvector',tvector);
     display(filename)
-    if (norm(Rs-Rf,2) < 1E-15)
+    if (norm(Rs-Rf,2) < 1E-2)
         disp('----> SUCCESS! <------')
         save(filename,"ts","Rs","Us")
     else
         disp('error radial not working')
     end
+    plot(tf,Rf,'--r'); hold on; plot(ts,Rs,'ks');
 end
 %%
 testb = testb + 1;