Works nicely. Still missing proper boundary conditions that take into account fluxes.

ternary_frap.m

@@ -4,74 +4,35 @@
 % conditions are introduced. Integration of model via pdepe.
 
 %% Solve pde
-% 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);
+x = linspace(49, 80, 700);
+t = linspace(0, 6, 1000);
 sol = pdepe(0, @flory_hugg_a, @flory_ic, @flory_bc, x, t);
-
 %% Plotting
 figure(1); hold on;
-for i = 1:50 
-    cla; plot(x, phi_tot(x, -50, 1)); plot(x, sol(i, :)); pause(0.1);
+for i = 1:500 
+    cla; xlim([49, 52]); ylim([0, 1.5]);
+    plot(x, phi_tot(x, -50, 0.25)); plot(x, sol(i, :)); pause(0.01);
 end
 
 %% Plot and check derivatives of pt
 figure; hold on;
-% 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)
+x = linspace(40, 60, 100);
+plot(x, phi_tot(x, -50, 0.25));
+plot(x, gradient_analytical(x, -50, 0.25));
+
 %% Function definitions for pde solver
 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/(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);
+pt = phi_tot(x, -50, 0.25);
+gra_a = gradient_analytical(x, -50, 0.25);
+c = 1/100;
+f = (1.3-pt)/pt*(pt*dudx-u*gra_a);
+s = 0;
 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<50
         u0 = 0.0;
@@ -79,6 +40,7 @@ function u0 = flory_ic(x)
         u0 = 0.3;
     end
 %     u0 = 0.3;
+%  u0 = phi_tot(x, -50, 1);
 end
 
 function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
@@ -92,30 +54,6 @@ function p = phi_tot(x, a, b)
     p = (tanh(-(x+a)/b)+1)/2+0.3;
 end
 
-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 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 = gradient_analytical(x, a, b)
     grad = -(1-tanh(-(x+a)/b).^2)*1/b*0.5;
-end
-
-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