Works in 1D, use different mobility! Doesn't work in 2D/3D.

prob_laplace.m

@@ -2,7 +2,7 @@
 % Define parameters for tanh, according to Weber, Zwicker, Julicher, Lee.
 b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
 e = @(chi) sqrt(3/8*(chi-2));
-t = [0, 0.0001, 0.05, 1, 50];
+t = [0, 0.0001, 0.02, 0.05, 0.1, 1, 50];
 %% Try out different precisions
 % prec = [0.5, 0.5, 1, 2, 3.5, 5];
 fac = [1, 5, 10];%, 10, 100];
@@ -27,35 +27,40 @@ end
 
 %% Test partitioning == 1 compared to natural no mobility case.
 % Tg1: g0 = 1 everywhere. 
-
-Tpt1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-9), 0.5, e(7/3),...
+prec = linspace(0.2, 3, 6);
+parfor i = 1:length(prec)
+Tpt1(i) = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-15), 0.5, e(7/3),...
                                 0, 1, 300, 7, 0, 'Constituent'},...
-                     t, 0.5);
+                     t, prec(i));
 % Tpt1.x = Tga1.x;
-Tpt1.solve_tern_frap();
+Tpt1(i).solve_tern_frap();
+end
 %%
 % Rescale mobility, such that phi_tot=1 everywhere achieves same flux
 Gi = 1;%(1-Tpt1.phi_t(1))/(1-Tpt1.phi_t(1));
 % Go = Tpt1.phi_t(1)/Tpt1.phi_t(end)*(1-Tpt1.phi_t(end));
 Go = Tpt1.phi_t(end)/Tpt1.phi_t(1)*2;
 %%
-Tga1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-9), Tpt1.phi_t(1), 0,...
-                                 0, Gi, 300, 7, 0, 'Constituent'},...
-                                 t, 0.5);
-% Tga1.x = Tpt1.x;
-Tga1.solve_tern_frap();
+parfor i = 1:length(prec)
+    Tga1(i) = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-15), Tpt1(i).phi_t(1), 0,...
+                                     -0.83, Gi, 300, 7, 0, 'Constituent'},...
+                                     t, prec(i));
+    Tga1(i).x = Tpt1(i).x;
+    Tga1(i).solve_tern_frap();
+end
 %%
-Tpt1.plot_sim('plot', 1, 'r', 1.72)
-Tga1.plot_sim('plot', 1, 'k', 1)
-
+for i = 6:length(prec)
+Tpt1(i).plot_sim('plot', 1, 'g', 1)
+Tga1(i).plot_sim('plot', 1, 'k', 1.0)
+end
 %% Calculate flux for Tpt1 at t=0
-fac = 1;
-[dudx, x, sol, flux, g0, pt] = calc_flux(Tpt1, 1);
-[dudx1, x1, sol1, flux1, g01, pt1] = calc_flux(Tga1, 1);
+fac = 1;%0.585;
+[dudx, x, sol, flux, g0, pt] = calc_flux(Tpt1(1), 2);
+[dudx1, x1, sol1, flux1, g01, pt1] = calc_flux(Tga1(1), 2);
 figure; hold on; title('flux');
 plot(x, flux); 
 plot(x1, fac*flux1);
-axis([0.90, 1.11, 0, inf]);
+axis([0.70, 1.51, 0, inf]);
 % figure; hold on; title('flux stretched');
 % plot(fac*(x-1), flux); 
 % plot(x1-1, fac*flux1);