diff --git a/.gitignore b/.gitignore
index c5f2c33740a1f552d746ddc3d50e7608e54fa272..d318f8caa52a07cd619561d1a79be474c5d9f2e1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
.ipynb_checkpoints/
+*jitfailure*
.DS_Store
__pycache__/
*~
diff --git a/@Ternary_model/create_mesh.m b/@Ternary_model/create_mesh.m
index 11c782857bbb3a0ee3e3f73e6f603a66d16eaff3..fed4a54abef0480a0436096233585d4c9b2ac1bd 100644
--- a/@Ternary_model/create_mesh.m
+++ b/@Ternary_model/create_mesh.m
@@ -3,32 +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
+% % 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.02, 2000), linspace(1.021, 2, 100),...
+ linspace(2.01, 10, 100), linspace(10.01, 300, 200)];
elseif strcmp(T.mode, 'Client')
T.x = linspace(0, 20, 24000);
end
diff --git a/@Ternary_model/plot_sim.m b/@Ternary_model/plot_sim.m
index 218d1752d719a57136c36eb0d1b6076fe828a3c4..0039d5e32fbf72ab9db219b1a09fa1cfc8578339 100644
--- a/@Ternary_model/plot_sim.m
+++ b/@Ternary_model/plot_sim.m
@@ -25,13 +25,13 @@ if strcmp(mode, 'mov')
end
elseif strcmp(mode, 'plot')
figure(20);
- cla;
+% cla;
hold on;
- xlim([-inf, 6.0]); ylim([-inf, max(T.sol(:))]);
+ xlim([0.6, 1.4]); ylim([-inf, 1]);
ax = gca;
ax.FontSize = 12;
xlabel('x [\mum]'); ylabel('volume fraction');
- plot(T.x, T.sol([1, 2:nth:si(1)], :), 'LineWidth', 1, 'Color', color);
+ plot(T.x, 3.47*T.sol([1, 2:nth:si(1)], :), 'LineWidth', 1, 'Color', color);
% plot(T.x, T.sol(nth, :), 'LineWidth', 1, 'Color', color);
plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0, 0),...
'LineWidth', 1, 'LineStyle', '--', 'Color', [0.97, 0.55, 0.03]);
diff --git a/@Ternary_model/solve_tern_frap.m b/@Ternary_model/solve_tern_frap.m
index 7d519c8dbb086f45985e72ebb5c38ea2a72a3618..837fdbcc4c57f31a5eb233ce1dbe59d3129da102 100644
--- a/@Ternary_model/solve_tern_frap.m
+++ b/@Ternary_model/solve_tern_frap.m
@@ -10,25 +10,7 @@ T.sol = pdepe(T.geometry, fh_pde, fh_ic, fh_bc, T.x, T.t);
end
%% Function definitions for pde solver
-function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
- e_g0, u_g0, v, mode)
-% Solve with full ternary model. Analytical derivatives.
-% pt ... phi_tot
-% gra_a ... analytical gradient of phi_tot
-pt = Ternary_model.phi_tot(x, a, b, e, u0, v*t);
-gra_a = Ternary_model.gradient_analytical(x, a, b, e, v*t);
-s = 0;
-c = 1;
-if strcmp(mode, 'Constituent')
- g0 = Ternary_model.gamma0(x, a+t*v, b, e_g0, u_g0, v*t);
- f = g0*(1-pt)/pt*(pt*dudx-u*gra_a);
-elseif strcmp(mode, 'Client')
- g0 = 0.5;
- chi_phi = -4.530864768482371;
- f = g0*(dudx+chi_phi*u*gra_a);
-end
-end
function u = flory_ic(x, a, b, u0, e, ic, x0)
if strcmp(ic, 'FRAP')
diff --git a/flory_hugg_pde.m b/flory_hugg_pde.m
new file mode 100644
index 0000000000000000000000000000000000000000..978587aab2455e1a677e63360a04c85afd1fa331
--- /dev/null
+++ b/flory_hugg_pde.m
@@ -0,0 +1,19 @@
+function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
+ e_g0, u_g0, v, mode)
+% Solve with full ternary model. Analytical derivatives.
+% pt ... phi_tot
+% gra_a ... analytical gradient of phi_tot
+
+pt = Ternary_model.phi_tot(x, a, b, e, u0, v*t);
+gra_a = Ternary_model.gradient_analytical(x, a, b, e, v*t);
+s = 0;
+c = 1;
+if strcmp(mode, 'Constituent')
+ g0 = Ternary_model.gamma0(x, a+t*v, b, e_g0, u_g0, v*t);
+ f = g0.*(1-pt)./pt.*(pt.*dudx-u.*gra_a);
+elseif strcmp(mode, 'Client')
+ g0 = 0.5;
+ chi_phi = -4.530864768482371;
+ f = g0*(dudx+chi_phi*u*gra_a);
+end
+end
\ No newline at end of file
diff --git a/prob_laplace.m b/prob_laplace.m
index 321695f318e696291379a5995078fafed5f4b845..2c4a3c0e9fbe9f6d58f55c32c5dcd08a2a60a193 100644
--- a/prob_laplace.m
+++ b/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.01, 0.05, 50];
+t = [0, 0.0001, 0.05, 1, 50];
%% Try out different precisions
% prec = [0.5, 0.5, 1, 2, 3.5, 5];
fac = [1, 5, 10];%, 10, 100];
@@ -30,16 +30,28 @@ end
Tpt1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-18), 0.5, e(7/3),...
0, 1, 300, 7, 0, 'Constituent'},...
- t, 1);
+ t, 0.3);
Tpt1.solve_tern_frap();
+%% Calculate flux for Tpt1 at t=0
+[dudx, x, sol, flux] = calc_flux(Tga1, 4);
+[dudx1, x1, sol1, flux1] = calc_flux(Tpt1, 4);
+figure; hold on; title('flux');
+plot(x, flux);
+plot(x1, 3.47*flux1);
+axis([0.90, 1.11, 0, inf]);
+figure; hold on; title('dudx');
+plot(x, dudx);
+plot(x1, 3.47*dudx1);
+axis([0.99, 1.01, 0, inf]);
%%
% Rescale mobility, such that phi_tot=1 everywhere achieves same flux
-Gi = (1-Tpt1.phi_t(1))/(1-Tpt1.phi_t(1));
-Go = Tpt1.phi_t(1)/Tpt1.phi_t(end)*(1-Tpt1.phi_t(end));
+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);
%%
Tga1 = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-18), Tpt1.phi_t(1), 0,...
- (Go-Gi), Gi, 300, 7, 0, 'Constituent'},...
- t, 1);
+ -(Gi-Go)/2, Gi, 300, 7, 0, 'Constituent'},...
+ t, 5.5);
Tga1.solve_tern_frap();
%%
Tpt1.plot_sim('plot', 1, 'g')
@@ -313,4 +325,14 @@ for i = 1:length(delta_x0)
end
end
p = sum(p_i);
+end
+
+function [dudx, x, sol, flux] = calc_flux(Temp, t_ind)
+
+dudx = diff(Temp.sol(t_ind, :))./diff(Temp.x);
+x = (Temp.x(1:end-1)+Temp.x(2:end))/2;
+sol = (Temp.sol(t_ind, 1:end-1)+Temp.sol(t_ind, 2:end))/2;
+
+[~, flux ,~] = flory_hugg_pde(x, 0, sol, dudx, Temp.a, Temp.b, Temp.e, ...
+ Temp.u0, Temp.e_g0, Temp.u_g0, Temp.v, Temp.mode);
end
\ No newline at end of file