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/Ternary_model.m b/@Ternary_model/Ternary_model.m
index aa39745ecf0c188ebae5a6ae5f10b0da8ead9e3d..0d53ba46c6213bb0a1b3531b34d4af39158f35cf 100644
--- a/@Ternary_model/Ternary_model.m
+++ b/@Ternary_model/Ternary_model.m
@@ -5,7 +5,7 @@ classdef Ternary_model < handle
% conditions are introduced. Integration of model via pdepe.
properties
- pa = '~/Desktop/DropletFRAP/Latex/Figures/FigModel/';
+ pa = '~/Desktop/DropletFRAP/Latex/Figures/Fig2/';
geometry = 2; % 2 ... Radial, 1 ... Cylindrical, 0 ... Slab
ic %'FRAP', 'Gauss', 'Uniform'
mode % 'Constituent': LLPS molecules move. 'Client': client moves
@@ -52,7 +52,7 @@ classdef Ternary_model < handle
T.calc_real_params;
end
create_mesh(T)
- plot_sim(T, mode, nth, color)
+ plot_sim(T, mode, nth, color, N, i)
solve_tern_frap(T)
calc_time_scales(T);
calc_real_params(T);
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/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index 6aa45f9771ab894030b8bbfd3941e2c8dfa77a08..c87251276bca4db12ab6ed1f416f557d9c8822bb 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -13,49 +13,306 @@
"import time\n",
"df.set_log_level(40)\n",
"\n",
- "domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
- "# mesh = ms.generate_mesh(domain, 35)\n",
- "mesh = df.UnitSquareMesh(50,50)\n",
- "dt = 0.01\n",
+ "# domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
+ "# mesh = ms.generate_mesh(domain, 50)\n",
+ "mesh = df.UnitIntervalMesh(10000)\n",
+ "F = df.FunctionSpace(mesh, 'CG', 1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def calc_sim(c0, c_tot, Ga0, dt, n_t, sym):\n",
+ " tc = df.TestFunction(F)\n",
+ " c = df.Function(F)\n",
+ " X = df.SpatialCoordinate(mesh)\n",
+ "# X.interpolate(df.Expression('x[0]', degree=1))\n",
+ " if sym == 0:\n",
+ " # Weak form 1D:\n",
+ "# form = (df.inner((c-c0)/dt, tc) +\n",
+ "# df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc)) -\n",
+ "# df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc))-\n",
+ "# tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0))+\n",
+ "# tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0))) * df.dx\n",
+ "# # Weak form 1D short:\n",
+ " form = ((c-c0)/dt*tc + (1-c_tot)*Ga0*\n",
+ " df.inner((df.grad(c)-c/c_tot*df.grad(c_tot)),\n",
+ " df.grad(tc))) * df.dx\n",
+ " elif sym == 2:\n",
+ " # Weak form radial symmetry:\n",
+ "# form = (df.inner((c-c0)/dt, tc*X[0]*X[0]) +\n",
+ "# df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc*X[0]*X[0])) -\n",
+ "# df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc*X[0]*X[0]))-\n",
+ "# tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*X[0]*X[0]))+\n",
+ "# tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0*X[0]*X[0]))-\n",
+ "# (1-c_tot)*Ga0*2*X[0]*c.dx(0)*tc+\n",
+ "# (1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx\n",
+ "# # Weak form radial symmetry:\n",
+ "# form = ((c-c0)/dt*tc*X[0]**2 +\n",
+ "# (1-c_tot)*Ga0*(c.dx(0)-c/c_tot*c_tot.dx(0))*\n",
+ "# tc.dx(0)*X[0]**2) * df.dx\n",
+ " form = ((c-c0)/dt*tc + (1-c_tot)*Ga0*\n",
+ " df.inner((df.grad(c)-c/c_tot*df.grad(c_tot)),\n",
+ " df.grad(tc)))*X[0]*X[0]*df.dx\n",
+ " t = 0\n",
+ " # Solve in time\n",
+ " ti = time.time()\n",
+ " for i in range(n_t):\n",
+ "# print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)]))\n",
+ " df.solve(form == 0, c)\n",
+ " df.assign(c0, c)\n",
+ " t += dt\n",
+ " print(time.time() - ti)\n",
+ " return c0\n",
"\n",
- "F = df.FunctionSpace(mesh, 'CG', 1)\n",
+ "def p_tot(p_i, p_o):\n",
+ " return str(p_i-p_o)+'*(-0.5*tanh(3500*(x[0]-0.1))+0.5)+'+str(p_o)\n",
+ "# return '(x[0]<0.1 ? '+str(p_i)+':'+ str(p_o)+')'\n",
+ "\n",
+ "def create_func(f_space, expr_str, deg):\n",
+ " f = df.Function(f_space)\n",
+ " f.interpolate(df.Expression(expr_str, degree=deg))\n",
+ " return f\n",
+ "\n",
+ "def eval_func(func, x):\n",
+ " return np.array([func([x]) for x in x])\n",
+ "\n",
+ "def eval_P(func, x):\n",
+ " return func(x[0])/func(x[1])\n",
+ "\n",
+ "def eval_D(func_ga, func_tot, x):\n",
+ " return(func_ga(x)*(1-func_tot(x)))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "bp = 0.1 # Boundary position at IC\n",
+ "sym = 2 # Symmetry of the problem\n",
+ "dt = 0.0001\n",
+ "\n",
+ "c_tot_1 = create_func(F, p_tot(0.9, 0.9), 1)\n",
+ "c0_1 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' + p_tot(0.9, 0.9), 1)\n",
+ "Ga0_1 = create_func(F, p_tot(1, 1/9), 1)\n",
+ "\n",
+ "c_tot_2 = create_func(F, p_tot(0.9, 0.45), 1)\n",
+ "c0_2 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.45), 1)\n",
+ "Ga0_2 = create_func(F,p_tot(1, 0.08), 1)\n",
+ "\n",
+ "c_tot_3 = create_func(F, p_tot(0.9, 0.3), 1)\n",
+ "c0_3 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.3), 1)\n",
+ "Ga0_3 = create_func(F,p_tot(1, 0.145), 1)\n",
+ "\n",
+ "c_tot_5 = create_func(F, p_tot(0.9, 0.18), 1)\n",
+ "c0_5 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.18), 1)\n",
+ "Ga0_5 = create_func(F,p_tot(1, 0.34), 1)\n",
+ "\n",
+ "c_tot_8 = create_func(F, p_tot(0.9, 0.1125), 1)\n",
+ "c0_8 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1125), 1)\n",
+ "Ga0_8 = create_func(F,p_tot(1, 0.8), 1)\n",
+ "\n",
+ "c_tot_9 = create_func(F, p_tot(0.9, 0.1), 1)\n",
+ "c0_9 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1), 1)\n",
+ "Ga0_9 = create_func(F,p_tot(1, 1), 1)\n",
+ "\n",
+ "c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, dt, 10, sym)\n",
+ "c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, dt, 10, sym)\n",
+ "c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, dt, 10, sym)\n",
+ "c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, dt, 10, sym)\n",
+ "c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, dt, 10, sym)\n",
+ "c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, dt, 10, sym)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x = np.linspace(0, 1, 10000)\n",
+ "plt.plot(x, eval_func(c0_1, x))\n",
+ "plt.plot(x, eval_func(c0_2, x))\n",
+ "plt.plot(x, eval_func(c0_3, x))\n",
+ "plt.plot(x, eval_func(c0_5, x))\n",
+ "plt.plot(x, eval_func(c0_8, x))\n",
+ "plt.plot(x, eval_func(c0_9, x))\n",
+ "# plt.xlim(0.08, 0.125125)\n",
+ "# plt.ylim(0., 0.325125)\n",
+ "plt.show()\n",
+ "# Diffusion versus partitioning\n",
+ "list_Ga = [Ga0_1, Ga0_2, Ga0_3, Ga0_5, Ga0_8, Ga0_9]\n",
+ "list_Tot = [c_tot_1, c_tot_2, c_tot_3, c_tot_5, c_tot_8, c_tot_9]\n",
+ "D = [eval_D(Ga, Tot, 1) for (Ga, Tot) in zip (list_Ga, list_Tot)]\n",
+ "P = [eval_P(Tot, [0, 1]) for Tot in list_Tot]\n",
+ "plt.plot(D, P)\n",
+ "plt.ylim(0, 11); plt.xlim(0, 1)\n",
+ "plt.ylabel('P'); plt.xlabel('Diffusion coefficient')\n",
+ "plt.plot(np.linspace(0, 1, 100), 9.5*(np.linspace(0, 1, 100))**(1/2))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plot outside, check invariance\n",
+ "x = np.linspace(bp, 1, 1000)\n",
+ "y1 = eval_func(c0_1, x)\n",
+ "y2 = eval_func(c0_2, x)\n",
+ "y3 = eval_func(c0_3, x)\n",
+ "y9 = eval_func(c0_9, x)\n",
+ "plt.plot((x-bp), y1)\n",
+ "plt.plot((x-bp)/2, 2*y2)\n",
+ "plt.plot((x-bp)/3, 3*y3)\n",
+ "plt.plot((x-bp)/9, 9*y9)\n",
+ "plt.xlim(0.0, 0.02)\n",
+ "plt.ylim(0.2, 1)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Plot inside, check invariance\n",
+ "x = np.linspace(0, 0.1, 1000)\n",
+ "y1 = eval_func(c0_1, x)-np.min(eval_func(c0_1, x))\n",
+ "y2 = eval_func(c0_2, x)-np.min(eval_func(c0_2, x))\n",
+ "y3 = eval_func(c0_3, x)-np.min(eval_func(c0_3, x))\n",
+ "y9 = eval_func(c0_9, x)-np.min(eval_func(c0_9, x))\n",
+ "plt.plot((x-np.min(x)), y1/eval_func(c0_1, [0.08]))\n",
+ "plt.plot((x-np.min(x)), y2/eval_func(c0_2, [0.08]))\n",
+ "plt.plot((x-np.min(x)), y3/eval_func(c0_3, [0.08]))\n",
+ "plt.plot((x-np.min(x)), y9/eval_func(c0_9, [0.08]))\n",
+ "plt.xlim(0.08, bp)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Check partitioning\n",
+ "cp = eval_func(c0_9, x)\n",
+ "print('Partitioning: ' + str(np.max(cp[1:1050])/np.min(cp[999:1050])))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Set B=1\n",
+ "pi = 0.8\n",
+ "po = 0.5 - np.sqrt(0.25 + (pi**2-pi))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p1_i = 0.9; p1_o = 0.1\n",
+ "p2_i = 0.8; p2_o = 0.2\n",
+ "p3_i = 0.7; p3_o = 0.3\n",
+ "p4_i = 0.9; p4_o = 0.9\n",
+ "\n",
+ "ct_1 = create_func(F, p_tot(p1_i, p1_o), 1)\n",
+ "ct_2 = create_func(F, p_tot(p2_i, p2_o), 1)\n",
+ "ct_3 = create_func(F, p_tot(p3_i, p3_o), 1)\n",
+ "ct_4 = create_func(F, p_tot(p4_i, p4_o), 1)\n",
"\n",
+ "g_1= create_func(F, '1', 1)\n",
+ "g_2= create_func(F, '0.5', 1)\n",
+ "g_3= create_func(F, '0.33333333333333333333', 1)\n",
+ "g_4= create_func(F, '1', 1)\n",
+ "\n",
+ "c0_1 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p1_i, p1_o), 1)\n",
+ "c0_2 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p2_i, p2_o), 1)\n",
+ "c0_3 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p3_i, p3_o), 1)\n",
+ "c0_4 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p4_i, p4_o), 1)\n",
+ "for i in range(10):\n",
+ " c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 2, 2)\n",
+ "# c0_2 = calc_sim(c0_2, ct_2, g_2, 0.01, 2, 2)\n",
+ "# c0_3 = calc_sim(c0_3, ct_3, g_3, 0.01, 2, 2)\n",
+ " c0_4 = calc_sim(c0_4, ct_4, g_4, 0.001, 2, 2)\n",
+ " plt.plot(x, eval_func(c0_1, x), 'r')#/eval_func(c0_1, [0.099])\n",
+ "# plt.plot(x, eval_func(c0_2, x), 'g')\n",
+ "# plt.plot(x, eval_func(c0_3, x)/eval_func(c0_3, [0.099]), 'b')\n",
+ " plt.plot(x, eval_func(c0_4, x), 'k')\n",
+ " plt.xlim(0.0, 0.625)\n",
+ " plt.ylim(0, 1.1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(x, eval_func(c0_1, x)/0.9)\n",
+ "plt.plot(x, eval_func(c0_2, x)/0.8)\n",
+ "plt.plot(x, eval_func(c0_3, x)/0.7)\n",
+ "plt.xlim(0.0, 0.125)\n",
+ "plt.ylim(0, 1.1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Radial diffusion equation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mesh = df.UnitIntervalMesh(1000)\n",
+ "dt = 0.001\n",
+ "F = df.FunctionSpace(mesh, 'CG', 1)\n",
"c0 = df.Function(F)\n",
- "c_tot = df.Function(F)\n",
+ "c0.interpolate(df.Expression('x[0]<0.5 && x[0]>0.2 ? 1:0', degree=1))\n",
+ "q = df.TestFunction(F)\n",
"c = df.Function(F)\n",
- "tc = df.TestFunction(F)\n",
- "\n",
- "# Interpolate c_tot and initial conditions\n",
- "# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.5)) ', degree=1))\n",
- "c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\n",
- "# tanh from left to right initial conditions\n",
- "# c_tot.interpolate(df.Expression('0.4*tanh(5*(x[0]-0.5)) - 1', degree=1))\n",
- "# Step function left to right initial condition\n",
- "# c0.interpolate(df.Expression('(x[0]<0.5) ? 0 : 1.0', degree=1))\n",
- "# Somewhat realistic half FRAP\n",
- "c0.interpolate(df.Expression(('(x[0]<0.5) && sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))<0.2 ? 0 :'\n",
- " '0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5'),\n",
- " degree=1))\n",
- "\n",
- "# Weak form\n",
- "form = ((df.inner((c-c0)/dt, tc) + (1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -\n",
- " (1-c_tot)/c_tot* c*df.inner(df.grad(c_tot), df.grad(tc))) * df.dx\n",
- "\n",
- "# Open file for writing\n",
- "cFile = df.XDMFFile('conc.xdmf')\n",
+ "X = df.SpatialCoordinate(mesh)\n",
+ "g = df.Expression('.00', degree=1)\n",
+ "u_D = df.Expression('1', degree=1)\n",
+ "def boundary(x, on_boundary):\n",
+ " return on_boundary\n",
+ "bc = df.DirichletBC(F, u_D, boundary)\n",
+ "# Weak form spherical symmetry\n",
+ "form = ((c-c0)/dt*q*X[0]*X[0] +\n",
+ " X[0]*X[0]*df.inner(df.grad(c), df.grad(q))-\n",
+ " c.dx(0)*2*X[0]*q) * df.dx\n",
+ "# Weak form 1D\n",
+ "# form = ((c-c0)/dt* q + df.inner(df.grad(c), df.grad(q))) * df.dx\n",
+ "# Weak form 1D with .dx(0) notation for derivative in 1st direction.\n",
+ "# form = ((c-c0)/dt*q + c.dx(0)*q.dx(0)) * df.dx\n",
"t = 0\n",
- "cFile.write(c0, t)\n",
"\n",
"# Solve in time\n",
- "ti = time.time()\n",
- "for i in range(80):\n",
- " print(t)\n",
+ "for i in range(60):\n",
+ " print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)]))\n",
" df.solve(form == 0, c)\n",
" df.assign(c0, c)\n",
" t += dt\n",
- " cFile.write(c0, t)\n",
- "cFile.close()\n",
- "print(time.time() - ti)"
+ " plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])"
]
},
{
@@ -63,11 +320,7 @@
"execution_count": null,
"metadata": {},
"outputs": [],
- "source": [
- "# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.1)) ', degree=1))\n",
- "c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\n",
- "plt.plot(np.linspace(0, 1, 100), [c_tot([x, 0.5]) for x in np.linspace(0, 1, 100)])"
- ]
+ "source": []
}
],
"metadata": {
diff --git a/Plots_Droplet_FRAP.ipynb b/Plots_Droplet_FRAP.ipynb
index 434198bad2cd687014854c82aeb14bc89c6a9844..0e860b07cd52a0996a3f98cf0461347e39bb1d45 100644
--- a/Plots_Droplet_FRAP.ipynb
+++ b/Plots_Droplet_FRAP.ipynb
@@ -72,7 +72,7 @@
"G_in = [1, 1, 1, .1, .1, .1]\n",
"G_out = [1, 1, 1, 0.99/0.9, 0.99/0.9, 0.99/0.9]\n",
"\n",
- "f_i_01 = []\n",
+ "f_i = []\n",
"\n",
"for p, m, p_i, p_e, G_i, G_o in zip(point_lists, me, phi_tot_int,\n",
" phi_tot_ext, G_in, G_out):\n",
diff --git a/flory_hugg_pde.m b/flory_hugg_pde.m
new file mode 100644
index 0000000000000000000000000000000000000000..981342d5b0103ffc060535b5ab61394ad33ef77b
--- /dev/null
+++ b/flory_hugg_pde.m
@@ -0,0 +1,19 @@
+function [c, f ,s, g0, pt] = 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 9bcb17baf3e19133edf416468de63a901199991a..be323698344900ac22ef6f980cf3eae7de26f770 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.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,24 +27,56 @@ 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),...
+prec = linspace(0.2, 3, 6);
+parfor i = 1:length(prec)
+Tpt1(i) = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-18), 0.5, e(7/3),...
0, 1, 300, 7, 0, 'Constituent'},...
- t, 1);
-Tpt1.solve_tern_frap();
+ t, prec(i));
+Tpt1(i).solve_tern_frap();
+end
%%
% 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)*2;
%%
-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);
-Tga1.solve_tern_frap();
+parfor i = 1:length(prec)
+ Tga1(i) = Ternary_model(0, 'FRAP', {-1, b(7/3, 10^-18), Tpt1(i).phi_t(1), 0,...
+ -0.83, Gi, 300, 7, 0, 'Constituent'},...
+ t, prec(i));
+ Tga1(i).solve_tern_frap();
+end
%%
-Tpt1.plot_sim('plot', 1, 'g')
+for i = 1:length(prec)
+Tpt1(i).plot_sim('plot', 1, 'g', 1, i)
+Tga1(i).plot_sim('plot', 1, 'k', 1.0, i)
+end
+%% Calculate flux for Tpt1 at t=0
+i = 6; ti=3;
+fac = 1;%0.585;
+[dudx, x, sol, flux, g0, pt] = calc_flux(Tpt1(i), ti);
+[dudx1, x1, sol1, flux1, g01, pt1] = calc_flux(Tga1(i), ti);
+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(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]);
+figure; hold on; title('Diffusion coefficients');
+plot(x, g0.*(1-pt));
+plot(x1, g01.*(1-pt1));
+axis([0.99, 1.01, 0, inf]);
%%
-Tga1.plot_sim('plot', 1, 'k')
+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)*
@@ -353,4 +385,14 @@ for i = 1:length(delta_x0)
end
end
p = sum(p_i);
+end
+
+function [dudx, x, sol, flux, g0, pt] = 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 ,~, g0, pt] = 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
diff --git a/ternary_frap.m b/ternary_frap.m
index b05f8825fd45b0849e00aa3e48f3820af2fdb5a5..0a8cb286b2437f34acef58ad2eac2487563d7dd1 100644
--- a/ternary_frap.m
+++ b/ternary_frap.m
@@ -10,12 +10,12 @@ make_graph_pretty('time', 'ratio outside/inside', 'bleached component');
%% Figures
% Time course
figure; hold on;
-xlim([0, 3]); ylim([0, 0.5]);
+xlim([0, 3]); ylim([0, 1]);
ax = gca;
ax.FontSize = 16;
xlabel('x [\mum]'); ylabel('volume fraction');
-plot(T.x, T.sol(1:2:200, :), 'LineWidth', 2, 'Color', [135/255 204/255 250/255]);
-plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 4,...
+plot(T.x, T.sol(1:1:100, :), 'LineWidth', 2, 'Color', [135/255 204/255 250/255]);
+plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0, 0), 'LineWidth', 4,...
'LineStyle', '--', 'Color', [247/255, 139/255, 7/255]);
print([T.pa, 'ternary_time_course'], '-depsc');
%% Plot and check derivatives of pt
@@ -32,14 +32,16 @@ plot(Ternary_model.gamma0(T.x, T.a, T.b, T.e_g0, T.o_g0),...
Ternary_model.spacing(Ternary_model.gamma0(T.x, T.a, T.b, T.e_g0, T.o_g0),...
T.a, T.b, T.e_g0, T.o_g0));
%% Fit with experimental BCs
-T = Ternary_model(2, 'FRAP', [], linspace(0.0, 400, 5000));
+T = Ternary_model(2, 'FRAP', {-1, b(7/3, 10^-6), 0.5, e(7/3),...
+ 0, 1, 100, 0, 0, 'Constituent'},...
+ linspace(0.0, 4000, 5000), 1);
T.solve_tern_frap();
-x_seed = 60;
-s = T.sol(1:end, 1:30)./max(T.sol(1:end, 1:30))/1.022;
-ts = T.t(1:end)';
-f_min = @(x) to_minimize(s, ts, x, 'simple_drop', s(:, end));
-opt = optimset('MaxFunEvals', 2000, 'PlotFcns',@optimplotfval);
-x_seed = fminsearch(f_min, x_seed, opt);
+% x_seed = 60;
+% s = T.sol(1:end, 1:30)./max(T.sol(1:end, 1:30))/1.022;
+% ts = T.t(1:end)';
+% f_min = @(x) to_minimize(s, ts, x, 'simple_drop', s(:, end));
+% opt = optimset('MaxFunEvals', 2000, 'PlotFcns',@optimplotfval);
+% x_seed = fminsearch(f_min, x_seed, opt);
%%
[cost, v_fit, r_n, v_fit_end] = to_minimize(s, ts, x_seed, 'simple_drop', s(:, end));
%% Plot for