Boundary conditions, initial conditions, phi_tot are now all parameterised.

ternary_frap.m

@@ -2,16 +2,22 @@
 % unbleached species. Model is assumed to be equilibrated
 % (bleached+unbleached=const.=pt). Then bleached species initial
 % conditions are introduced. Integration of model via pdepe.
-
+a = -50;
+b = 0.05;
+c = 1/1000000000;
+u0 = 0.1;
 %% Solve pde
-x = linspace(49, 80, 700);
-t = linspace(0, 100, 1000);
-sol = pdepe(0, @flory_hugg_a, @flory_ic, @flory_bc, x, t);
+x = [linspace(49, 55, 600), linspace(55.01, 200, 300)];
+t = linspace(0, 100, 100);
+fh_ic = @(x) flory_ic(x, u0);
+fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, u0);
+fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, a, b, c, u0);
+sol = pdepe(0, fh_pde, fh_ic, fh_bc, x, t);
 %% Plotting
 figure(1); hold on;
-for i = 1:500 
+for i = 1:100 
     cla; xlim([49, 53]); ylim([0, 1.5]);
-    plot(x, phi_tot(x, -50, 0.25)); plot(x, sol(i, :)); pause(0.01);
+    plot(x, phi_tot(x, a, b, u0)); plot(x, sol(i, :)); pause(0.01);
 end
 
 %% Plot and check derivatives of pt
@@ -21,41 +27,41 @@ 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)  
+function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, c_p, u0)  
 % Solve with full ternary model. Analytical derivatives.
 % pt ... phi_tot
 % gra_a ... analytical gradient of phi_tot
 
-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);
+pt = phi_tot(x, a, b, u0);
+gra_a = gradient_analytical(x, a, b);
+c = c_p;
+f = ((u0+1)-pt)/pt*(pt*dudx-u*gra_a);
 s = 0;
 end
 
-function u0 = flory_ic(x)
+function u = flory_ic(x, u0)
     if x<50
-        u0 = 0.0;
+        u = 0.0;
     else
-        u0 = 0.3;
+        u = u0;
     end
 %     u0 = 0.3;
 %  u0 = phi_tot(x, -50, 1);
 end
 
-function [pl,ql,pr,qr] = flory_bc(xl,ul,xr,ur,t)
+function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0)
     pl = 0;
     ql = 1;
 %     % No flux
 %     pr = 0;%ur - 0.01;
 %     qr = 1;
     % Dirichlet BC
-    pr = ur- 0.3;
+    pr = ur - u0;
     qr = 0;
 end
 
-function p = phi_tot(x, a, b)
-    p = (tanh(-(x+a)/b)+1)/2+0.3;
+function p = phi_tot(x, a, b, u0)
+    p = (tanh(-(x+a)/b)+1)/2+u0;
 end
 
 function grad = gradient_analytical(x, a, b)