Problem independent from boundary width and accuracy.

@Ternary_model/create_mesh.m

@@ -3,34 +3,34 @@ function create_mesh(T)
 %   Mesh density larger for steeper gradients.
 
 if strcmp(T.mode, 'Constituent')
-% % Start with left side of mesh, at r=0/x=0, with a given step size
-%     x = linspace(0, 0.75, 40);
-%     x = [x, x(end)+0.001/T.precision];
-%     while x(end)<-T.a
-%         x_temp = x(end)-x(end-1);
-%         phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0, 0);
-%         phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
-%         while phi_nm1-phi_n > 0.002/T.precision
-%             x_temp = x_temp/2;
-%             phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
-%         end
-%         x = [x, x(end)+x_temp];
-%     end
-%     x = [x, x(end)+cumsum(fliplr(diff(x)))];
-%     x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 100*T.precision);
-%     x_right = linspace(-10*T.a+0.1, T.system_size, 100*T.precision);
-%     T.x = [x, x_middle, x_right];
-%     T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1);
-% 
-%     % Fine mesh close to x0 if ic = 'Gauss'
-%     [~, ind] = min(abs(T.x-T.x0));
-%     if (T.x(ind+20)-T.x(ind) > 0.1) && strcmp(T.ic, 'Gauss')
-%         T.x = [T.x(1:ind-21), linspace(T.x(ind-20), T.x(ind+20),...
-%                                        ceil((T.x(ind+20)-T.x(ind-20))*200)),...
-%                T.x(ind+21:end)];
-%     end
-T.x = [linspace(0, 0.98, 100), linspace(0.9801, 1.06, 2000), linspace(1.061, 2, 100),...
-       linspace(2.01, 10, 100), linspace(10.01, 300, 200)];
+% Start with left side of mesh, at r=0/x=0, with a given step size
+    x = linspace(0, 0.75, 40);
+    x = [x, x(end)+0.001/T.precision];
+    while x(end)<-T.a
+        x_temp = x(end)-x(end-1);
+        phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0, 0);
+        phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
+        while phi_nm1-phi_n > 0.002/T.precision
+            x_temp = x_temp/2;
+            phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0, 0);
+        end
+        x = [x, x(end)+x_temp];
+    end
+    x = [x, x(end)+cumsum(fliplr(diff(x)))];
+    x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 100*T.precision);
+    x_right = linspace(-10*T.a+0.1, T.system_size, 100*T.precision);
+    T.x = [x, x_middle, x_right];
+    T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1);
+
+    % Fine mesh close to x0 if ic = 'Gauss'
+    [~, ind] = min(abs(T.x-T.x0));
+    if (T.x(ind+20)-T.x(ind) > 0.1) && strcmp(T.ic, 'Gauss')
+        T.x = [T.x(1:ind-21), linspace(T.x(ind-20), T.x(ind+20),...
+                                       ceil((T.x(ind+20)-T.x(ind-20))*200)),...
+               T.x(ind+21:end)];
+    end
+% T.x = [linspace(0, 0.99, 100), linspace(0.9901, 1.01, 2000), linspace(1.0161, 2, 100),...
+%        linspace(2.01, 10, 100), linspace(10.01, 300, 200)];
 elseif strcmp(T.mode, 'Client')
     T.x = linspace(0, 20, 24000);
 end

prob_laplace.m

@@ -28,45 +28,51 @@ end
 %% Test partitioning == 1 compared to natural no mobility case.
 % Tg1: g0 = 1 everywhere. 
 
-Tpt1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-18), 0.5, e(7/3),...
+Tpt1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-9), 0.5, e(7/3),...
                                 0, 1, 300, 7, 0, 'Constituent'},...
-                     t, 0.3);
+                     t, 0.5);
+% Tpt1.x = Tga1.x;
 Tpt1.solve_tern_frap();
+%%
+% 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();
+%%
+Tpt1.plot_sim('plot', 1, 'r', 1.72)
+Tga1.plot_sim('plot', 1, 'k', 1)
+
 %% Calculate flux for Tpt1 at t=0
-fac = 3.58;
-[dudx, x, sol, flux, g0, pt] = calc_flux(Tpt1, 4);
-[dudx1, x1, sol1, flux1, g01, pt1] = calc_flux(Tga1, 4);
+fac = 1;
+[dudx, x, sol, flux, g0, pt] = calc_flux(Tpt1, 1);
+[dudx1, x1, sol1, flux1, g01, pt1] = calc_flux(Tga1, 1);
 figure; hold on; title('flux');
 plot(x, flux); 
 plot(x1, fac*flux1);
 axis([0.90, 1.11, 0, inf]);
-figure; hold on; title('flux stretched');
-plot(0.9*fac*(x-1), flux); 
-plot(x1-1, fac*flux1);
-axis([-0.1, 0.11, 0, inf]);
+% figure; hold on; title('flux stretched');
+% plot(fac*(x-1), flux); 
+% plot(x1-1, fac*flux1);
+% axis([-0.1, 0.11, 0, inf]);
 figure; hold on; title('dudx');
 plot(x, dudx); 
 plot(x1, dudx1);
 axis([0.99, 1.01, 0, inf]);
-fluxratio = flux1./flux;
 figure; hold on; title('Diffusion coefficients');
 plot(x, g0.*(1-pt));
 plot(x1, g01.*(1-pt1));
 axis([0.99, 1.01, 0, inf]);
 %%
-% 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^-18), Tpt1.phi_t(1), 0,...
-                                 0, Gi, 300, 7, 0, 'Constituent'},...
-                                 t, 5.5);
-Tga1.solve_tern_frap();
-%%
-Tpt1.plot_sim('plot', 1, 'g', 4.1)
-%%
-Tga1.plot_sim('plot', 1, 'k', 1)
+figure; hold on; title('sol');
+plot(x, sol); 
+plot(x1, sol1);
+axis([0.90, 1.11, 0, inf]);
 %%
 figure(3); hold on; 
 % fact = [0.9965, 0.94, 1.01, 1.98];fact(i)*