Redefining various bits. Most importantly: oscillations gone. Sol doesn't converge yet.

ternary_frap.m

@@ -4,55 +4,81 @@
 % conditions are introduced. Integration of model via pdepe.
 
 %% Solve pde
-x = linspace(0,500, 6000);
-t = linspace(0, 500, 50);
+% x = [linspace(0, 40, 100), linspace(40.4, 60, 50), linspace(60.4, 100, 100)];
+x = linspace(40, 60, 300);
+% x = linspace(0.5, 100.5, 1001);
+
+t = linspace(0, 500000, 50);
+% y = linspace(0.3000001, 1.2999999, 500);
+% x = 0.5*atanh(2*y-1.6)+50;
+a = -50;
+b = 0.5;
+% x_t = linspace(0.5, 101.5, 1011);
+% global phi_tot_g
+% phi_tot_g = phi_tot(x_t, b, a);
+% global der_phi 
+% der_phi = gradient_analytical(x_t, a, b);
+% global lapl_phi
+% lapl_phi = laplacian_analytical(x_t, a, b);
+
 % opt = odeset('RelTol',1e-14, 'AbsTol', 1e-16,'MaxStep',1e-2); 
 % sol = pdepe(0,@flory_pde, @flory_ic, @flory_bc,x,t, opt);
-sol = pdepe(0,@flory_hugg_pde, @flory_ic, @flory_bc,x,t);
+sol = pdepe(0, @flory_hugg_a, @flory_ic, @flory_bc, x, t);
 
 %% Plotting
 figure(1); hold on;
-for i = 1:50
-    cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.04);
+for i = 1:50 
+    cla; plot(x, phi_tot(x, -50, 1)); plot(x, sol(i, :)); pause(0.1);
 end
 
 %% Plot and check derivatives of pt
 figure; hold on;
-x = linspace(140, 170, 500);
-plot(x, pt(x)); plot(x, gra_pt(x, 0.001)); plot(x, lap_pt(x, 0.001));
-plot(x, gralap_pt(x, 0.001)); plot(x, laplap_pt(x, 0.001));
-plot(x, gra_a(x, -155, 0.5));
-plot(x, lap_a(x, -155, 0.5));
-% pt(0)
-% gra_pt(0, 0.00001)
-% lap_pt(0, 0.00001)
-% lap_a(0, -155, 0.5)
-% gralap_pt(0, 0.00001)
-% laplap_pt(0, 0.00001)
+% x = linspace(40, 60, 100);
+plot(x, phi_tot(x, -50, 0.5));% plot(x, gra_pt(x, -50, 0.5, 0.001)); plot(x, lap_pt(x, -50, 0.5, 0.001));
+% plot(x, gralap_pt(x, -50, 0.5, 0.001)); plot(x, laplap_pt(x, -50, 0.5, 0.001));
+plot(x, gradient_analytical(x, -50, 0.5));
+plot(x, laplacian_analytical(x, -50, 0.5));
+%%
+phi_tot(x, -50, 1)
 %% Function definitions for pde solver
-function [c, f ,s] = flory_hugg_pde(x, t, u, dudx)  
-% Solve with full ternary model.
-X = 3.2;
-K = 1.5;
-c = 1/(pt(x)+1);
+function [c, f ,s] = flory_hugg_a(x, t, u, dudx)  
+% Solve with full ternary model. Analytical derivatives.
+% pt ... phi_tot
+% gra_a ... analytical gradient of phi_tot
+% lap_a ... analytical laplacian of phi_tot
+
+pt = @(x) phi_tot(x, -50, 1);
+gra_a = @(x) gradient_analytical(x, -50, 1);
+lap_a = @(x) laplacian_analytical(x, -50, 1);
+% gra_a = @(x) gra_pt(x, -50, 0.5, 0.0001);
+% lap_a = @(x) lap_pt(x, -50, 0.5, 0.0001);
+c = 1/(1.3-pt(x));
 f = dudx;
-s = u*(lap_pt(x, 0.0001)-2*X*(lap_pt(x, 0.0001)-...
-        gra_pt(x, 0.0001)^2-pt(x)*lap_pt(x, 0.0001))+...
-        K*(-laplap_pt(x, 0.0001)+...
-            gra_pt(x, 0.0001)*gralap_pt(x, 0.0001)+...
-            pt(x)*laplap_pt(x, 0.0001)))/(pt(x)+1)+...
-    dudx*(-2*X*(gra_pt(x, 0.0001)-pt(x)*gra_pt(x, 0.0001))-...
-          K*gralap_pt(x, 0.0001)+K*pt(x)*gralap_pt(x, 0.0001))/...
-          (pt(x)+1);
+s = u/(1.3-pt(x))*(lap_a(x)+(gra_a(x)/pt(x))^2-lap_a(x)/pt(x))...
+    -dudx/(1.3-pt(x))/pt(x)*gra_a(x);
 end
 
+% function [c, f ,s] = flory_hugg_a(x, t, u, dudx)  
+% % Solve with full ternary model. 
+% global phi_tot_g
+% global der_phi 
+% global lapl_phi
+% 
+% i = round(10*x);
+% % disp(i)
+% c = 1/(1-phi_tot_g(i));
+% f = dudx;
+% s = u/(1-phi_tot_g(i))*(lapl_phi(i)+(der_phi(i)/phi_tot_g(i))^2-lapl_phi(i)/phi_tot_g(i))...
+%     -dudx/(1-phi_tot_g(i))/phi_tot_g(i)*der_phi(i);
+% end
+
 function u0 = flory_ic(x)
-    if x<5
+    if x<50
         u0 = 0.0;
     else
-        u0 = 0.5;
+        u0 = 0.3;
     end
-%     u0 = pt(x)*0.5;
+%     u0 = 0.3;
 end
 
 function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
@@ -62,42 +88,34 @@ function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
     qr = 1;
 end
 
-function [c, f ,s] = flory_pde(x, t, u, dudx) 
-% Tentatively solve system where binary mixture is assumed, with variable
-% '1' (now called pt). Basically pt is fixed and now phi_u and
-% phi_b can swap places but nothing else.
-c = 1/pt(x);
-f = dudx;
-s = -u*lap_pt(x)/pt(x);
-end
-
-function p = pt(x)
-    p = (tanh(-(x-155)*2)+1)/2+0.3;
-end
-
-function gpt = gra_pt(x, delta)
-    gpt = (pt(x+delta)-pt(x-delta))/(2*delta);
-end
-
-function lpt = lap_pt(x, delta)
-    lpt = (gra_pt(x+delta, delta)-...
-           gra_pt(x-delta, delta))/(2*delta);
+function p = phi_tot(x, a, b)
+    p = (tanh(-(x+a)/b)+1)/2+0.3;
 end
 
-function glpt = gralap_pt(x, delta)
-    glpt = (lap_pt(x+delta, delta)-...
-            lap_pt(x-delta, delta))/(2*delta);
+function gpt = gra_pt(x, a, b, delta)
+    gpt = (phi_tot(x+delta, a, b)-...
+           phi_tot(x-delta, a, b))/(2*delta);
 end
 
-function llpt = laplap_pt(x, delta)
-    llpt = (gralap_pt(x+delta, delta)-...
-            gralap_pt(x-delta, delta))/(2*delta);
+function lpt = lap_pt(x, a, b, delta)
+    lpt = (gra_pt(x+delta, a, b, delta)-...
+           gra_pt(x-delta, a, b, delta))/(2*delta);
 end
+% 
+% function glpt = gralap_pt(x, a, b, delta)
+%     glpt = (lap_pt(x+delta, a, b, delta)-...
+%             lap_pt(x-delta, a, b, delta))/(2*delta);
+% end
+% 
+% function llpt = laplap_pt(x, a, b, delta)
+%     llpt = (gralap_pt(x+delta, a, b, delta)-...
+%             gralap_pt(x-delta, a, b, delta))/(2*delta);
+% end
 
-function grad = gra_a(x, a, b)
+function grad = gradient_analytical(x, a, b)
     grad = -(1-tanh(-(x+a)/b).^2)*1/b*0.5;
 end
 
-function lap = lap_a(x, a, b)
+function lap = laplacian_analytical(x, a, b)
     lap = -2*tanh(-(x+a)/b).*(1-tanh(-(x+a)/b).^2)*1/b^2*0.5;
 end
\ No newline at end of file