Comparison ML Python works. Spherical still off. Slab fine.

a9c0d6b5 · Lars Hubatsch · 9b0346c3 · a9c0d6b5 · a9c0d6b5 · a9c0d6b5
Commit a9c0d6b5 authored 4 years ago by Lars Hubatsch
--- a/@Ternary_model/create_mesh.m
+++ b/@Ternary_model/create_mesh.m
@@ -4,7 +4,7 @@ function create_mesh(T)
 if strcmp(T.mode, 'Constituent') | T.v == 0
 % Start with left side of mesh, at r=0/x=0, with a given step size
-    x = linspace(0, 0.75, 40);
+    x = linspace(0, -0.75*T.a, 40);
    x = [x, x(end)+0.001/T.precision];
    while x(end)<-T.a
        x_temp = x(end)-x(end-1);

--- a/@Ternary_model/plot_sim.m
+++ b/@Ternary_model/plot_sim.m
@@ -12,26 +12,26 @@ if strcmp(mode, 'mov')
    for i = 1:nth:length(T.t)
        cla;
        axis([0, T.system_size, 0, max(T.sol(:))]);
-        ax = gca;
+        ax = gca; ax.FontSize = 12;
-        ax.FontSize = 12;
        xlabel('position'); ylabel('probability [not normalized]');
        plot(T.x, T.phi_tot(T.x, T.a, T.b, T.e, T.u0, T.t(i)*T.v), 'LineWidth', 2, ...
            'LineStyle', '--');
        plot(T.x, T.sol(i, :), 'LineWidth', 2);
        text(abs(-T.a+2), 0.7, num2str(T.t(i)));
-        shg
+        shg; pause();
-        pause();
 %         print([num2str(i),'.png'],'-dpng')
    end
 elseif strcmp(mode, 'plot')
    figure(20);
-    cla;
+%     cla;
    hold on;
-%     xlim([-inf, 6.0]); ylim([-inf, max(T.sol(:))]);
+    xlim([-inf, -2*T.a]); %ylim([-inf, max(T.sol(:))]);
    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);
+%     N = max(max(T.sol([1, 2:nth:si(1)], :)));
+    N = 1;
+    plot(T.x, T.sol([1, 1:nth:si(1)], :)/N, '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]);

--- a/@Ternary_model/solve_tern_frap.m
+++ b/@Ternary_model/solve_tern_frap.m
 function solve_tern_frap(T)
 %SOLVE_TERN_FRAP solves Fokker Planck equation for ternary FRAP
 % Assumes fixed profile of phi_tot
+tic
 fh_ic = @(x) flory_ic(x, T.a, T.b, T.u0, T.e, T.ic, T.x0);
 fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, T.u0, T.e);
 fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, T.a, T.b, T.e,...
                                          T.u0, T.e_g0, T.u_g0, T.v,...
                                          T.mode);
 T.sol = pdepe(T.geometry, fh_pde, fh_ic, fh_bc, T.x, T.t);
+toc
 end
 %% Function definitions for pde solver

--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -94,8 +94,9 @@
   "outputs": [],
   "source": [
    "bp = 0.1 # Boundary position at IC\n",
-    "sym = 2 # Symmetry of the problem\n",
+    "sym = 0 # Symmetry of the problem\n",
-    "dt = 0.0001\n",
+    "dt = 0.001\n",
+    "nt = 10\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",
@@ -121,12 +122,12 @@
    "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_1 = calc_sim(c0_1, c_tot_1, Ga0_1, dt, nt, sym)\n",
-    "c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, dt, 10, sym)\n",
+    "# c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, dt, nt, sym)\n",
-    "c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, dt, 10, sym)\n",
+    "# c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, dt, nt, sym)\n",
-    "c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, dt, 10, sym)\n",
+    "# c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, dt, nt, sym)\n",
-    "c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, dt, 10, sym)\n",
+    "# c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, dt, nt, sym)\n",
-    "c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, dt, 10, sym)"
+    "c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, dt, nt, sym)"
   ]
  },
  {
@@ -270,6 +271,32 @@
    "plt.ylim(0, 1.1)"
   ]
  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Comparison with Matlab"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ML_9 = np.genfromtxt('Matlab_workspaces/ML_9_1.csv', delimiter=',')\n",
+    "ML_1 = np.genfromtxt('Matlab_workspaces/ML_1_9.csv', delimiter=',')\n",
+    "plt.plot(ML_9[0, :], ML_9[-1, :])\n",
+    "plt.plot(x, eval_func(c0_9, x))\n",
+    "plt.xlim(0, 0.5)\n",
+    "# plt.show()\n",
+    "\n",
+    "plt.plot(ML_1[0, :], ML_1[-1, :])\n",
+    "plt.plot(x, eval_func(c0_1, x))\n",
+    "plt.xlim(0, 0.5)\n",
+    "plt.show()"
+   ]
+  },
  {
   "cell_type": "markdown",
   "metadata": {},

 %% Cell type:code id: tags:
 ``` python
 import dolfin as df
 import matplotlib.pyplot as plt
 import mshr as ms
 import numpy as np
 import time
 df.set_log_level(40)
 # domain = ms.Sphere(df.Point(0, 0, 0), 1.0)
 # mesh = ms.generate_mesh(domain, 50)
 mesh = df.UnitIntervalMesh(10000)
 F = df.FunctionSpace(mesh, 'CG', 1)
 ```
 %% Cell type:code id: tags:
 ``` python
 def calc_sim(c0, c_tot, Ga0, dt, n_t, sym):
    tc = df.TestFunction(F)
    c = df.Function(F)
    X = df.SpatialCoordinate(mesh)
 #     X.interpolate(df.Expression('x[0]', degree=1))
    if sym == 0:
    # Weak form 1D:
 #         form = (df.inner((c-c0)/dt, tc) +
 #              df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc)) -
 #              df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc))-
 #              tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0))+
 #              tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0))) * df.dx
 #     # Weak form 1D short:
        form = ((c-c0)/dt*tc + (1-c_tot)*Ga0*
                df.inner((df.grad(c)-c/c_tot*df.grad(c_tot)),
                          df.grad(tc))) * df.dx
    elif sym == 2:
    # Weak form radial symmetry:
 #         form = (df.inner((c-c0)/dt, tc*X[0]*X[0]) +
 #              df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc*X[0]*X[0])) -
 #              df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc*X[0]*X[0]))-
 #              tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*X[0]*X[0]))+
 #              tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0*X[0]*X[0]))-
 #              (1-c_tot)*Ga0*2*X[0]*c.dx(0)*tc+
 #              (1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx
 #     # Weak form radial symmetry:
 #         form = ((c-c0)/dt*tc*X[0]**2 +
 #                 (1-c_tot)*Ga0*(c.dx(0)-c/c_tot*c_tot.dx(0))*
 #                          tc.dx(0)*X[0]**2) * df.dx
        form = ((c-c0)/dt*tc + (1-c_tot)*Ga0*
                df.inner((df.grad(c)-c/c_tot*df.grad(c_tot)),
                          df.grad(tc)))*X[0]*X[0]*df.dx
    t = 0
    # Solve in time
    ti = time.time()
    for i in range(n_t):
 #         print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)]))
        df.solve(form == 0, c)
        df.assign(c0, c)
        t += dt
    print(time.time() - ti)
    return c0
 def p_tot(p_i, p_o):
    return str(p_i-p_o)+'*(-0.5*tanh(3500*(x[0]-0.1))+0.5)+'+str(p_o)
 #     return '(x[0]<0.1 ? '+str(p_i)+':'+ str(p_o)+')'
 def create_func(f_space, expr_str, deg):
    f = df.Function(f_space)
    f.interpolate(df.Expression(expr_str, degree=deg))
    return f
 def eval_func(func, x):
    return np.array([func([x]) for x in x])
 def eval_P(func, x):
    return func(x[0])/func(x[1])
 def eval_D(func_ga, func_tot, x):
    return(func_ga(x)*(1-func_tot(x)))
 ```
 %% Cell type:code id: tags:
 ``` python
 bp = 0.1 # Boundary position at IC
-sym = 2 # Symmetry of the problem
+sym = 0 # Symmetry of the problem
-dt = 0.0001
+dt = 0.001
+nt = 10
 c_tot_1 = create_func(F, p_tot(0.9, 0.9), 1)
 c0_1 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' + p_tot(0.9, 0.9), 1)
 Ga0_1 = create_func(F, p_tot(1, 1/9), 1)
 c_tot_2 = create_func(F, p_tot(0.9, 0.45), 1)
 c0_2 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.45), 1)
 Ga0_2 = create_func(F,p_tot(1, 0.08), 1)
 c_tot_3 = create_func(F, p_tot(0.9, 0.3), 1)
 c0_3 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.3), 1)
 Ga0_3 = create_func(F,p_tot(1, 0.145), 1)
 c_tot_5 = create_func(F, p_tot(0.9, 0.18), 1)
 c0_5 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.18), 1)
 Ga0_5 = create_func(F,p_tot(1, 0.34), 1)
 c_tot_8 = create_func(F, p_tot(0.9, 0.1125), 1)
 c0_8 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1125), 1)
 Ga0_8 = create_func(F,p_tot(1, 0.8), 1)
 c_tot_9 = create_func(F, p_tot(0.9, 0.1), 1)
 c0_9 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1), 1)
 Ga0_9 = create_func(F,p_tot(1, 1), 1)
-c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, dt, 10, sym)
+c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, dt, nt, sym)
-c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, dt, 10, sym)
+# c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, dt, nt, sym)
-c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, dt, 10, sym)
+# c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, dt, nt, sym)
-c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, dt, 10, sym)
+# c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, dt, nt, sym)
-c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, dt, 10, sym)
+# c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, dt, nt, sym)
-c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, dt, 10, sym)
+c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, dt, nt, sym)
 ```
 %% Cell type:code id: tags:
 ``` python
 x = np.linspace(0, 1, 10000)
 plt.plot(x, eval_func(c0_1, x))
 plt.plot(x, eval_func(c0_2, x))
 plt.plot(x, eval_func(c0_3, x))
 plt.plot(x, eval_func(c0_5, x))
 plt.plot(x, eval_func(c0_8, x))
 plt.plot(x, eval_func(c0_9, x))
 # plt.xlim(0.08, 0.125125)
 # plt.ylim(0., 0.325125)
 plt.show()
 # Diffusion versus partitioning
 list_Ga = [Ga0_1, Ga0_2, Ga0_3, Ga0_5, Ga0_8, Ga0_9]
 list_Tot = [c_tot_1, c_tot_2, c_tot_3, c_tot_5, c_tot_8, c_tot_9]
 D = [eval_D(Ga, Tot, 1) for (Ga, Tot) in zip (list_Ga, list_Tot)]
 P = [eval_P(Tot, [0, 1]) for Tot in list_Tot]
 plt.plot(D, P)
 plt.ylim(0, 11); plt.xlim(0, 1)
 plt.ylabel('P'); plt.xlabel('Diffusion coefficient')
 plt.plot(np.linspace(0, 1, 100), 9.5*(np.linspace(0, 1, 100))**(1/2))
 ```
 %% Cell type:code id: tags:
 ``` python
 # Plot outside, check invariance
 x = np.linspace(bp, 1, 1000)
 y1 = eval_func(c0_1, x)
 y2 = eval_func(c0_2, x)
 y3 = eval_func(c0_3, x)
 y9 = eval_func(c0_9, x)
 plt.plot((x-bp), y1)
 plt.plot((x-bp)/2, 2*y2)
 plt.plot((x-bp)/3, 3*y3)
 plt.plot((x-bp)/9, 9*y9)
 plt.xlim(0.0, 0.02)
 plt.ylim(0.2, 1)
 plt.show()
 ```
 %% Cell type:code id: tags:
 ``` python
 # Plot inside, check invariance
 x = np.linspace(0, 0.1, 1000)
 y1 = eval_func(c0_1, x)-np.min(eval_func(c0_1, x))
 y2 = eval_func(c0_2, x)-np.min(eval_func(c0_2, x))
 y3 = eval_func(c0_3, x)-np.min(eval_func(c0_3, x))
 y9 = eval_func(c0_9, x)-np.min(eval_func(c0_9, x))
 plt.plot((x-np.min(x)), y1/eval_func(c0_1, [0.08]))
 plt.plot((x-np.min(x)), y2/eval_func(c0_2, [0.08]))
 plt.plot((x-np.min(x)), y3/eval_func(c0_3, [0.08]))
 plt.plot((x-np.min(x)), y9/eval_func(c0_9, [0.08]))
 plt.xlim(0.08, bp)
 plt.show()
 ```
 %% Cell type:code id: tags:
 ``` python
 # Check partitioning
 cp = eval_func(c0_9, x)
 print('Partitioning: ' + str(np.max(cp[1:1050])/np.min(cp[999:1050])))
 ```
 %% Cell type:code id: tags:
 ``` python
 # Set B=1
 pi = 0.8
 po = 0.5 - np.sqrt(0.25 + (pi**2-pi))
 ```
 %% Cell type:code id: tags:
 ``` python
 p1_i = 0.9; p1_o = 0.1
 p2_i = 0.8; p2_o = 0.2
 p3_i = 0.7; p3_o = 0.3
 p4_i = 0.9; p4_o = 0.9
 ct_1 = create_func(F, p_tot(p1_i, p1_o), 1)
 ct_2 = create_func(F, p_tot(p2_i, p2_o), 1)
 ct_3 = create_func(F, p_tot(p3_i, p3_o), 1)
 ct_4 = create_func(F, p_tot(p4_i, p4_o), 1)
 g_1= create_func(F, '1', 1)
 g_2= create_func(F, '0.5', 1)
 g_3= create_func(F, '0.33333333333333333333', 1)
 g_4= create_func(F, '1', 1)
 c0_1 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p1_i, p1_o), 1)
 c0_2 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p2_i, p2_o), 1)
 c0_3 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p3_i, p3_o), 1)
 c0_4 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p4_i, p4_o), 1)
 for i in range(10):
    c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 2, 2)
 #     c0_2 = calc_sim(c0_2, ct_2, g_2, 0.01, 2, 2)
 #     c0_3 = calc_sim(c0_3, ct_3, g_3, 0.01, 2, 2)
    c0_4 = calc_sim(c0_4, ct_4, g_4, 0.001, 2, 2)
    plt.plot(x, eval_func(c0_1, x), 'r')#/eval_func(c0_1, [0.099])
 #     plt.plot(x, eval_func(c0_2, x), 'g')
 #     plt.plot(x, eval_func(c0_3, x)/eval_func(c0_3, [0.099]), 'b')
    plt.plot(x, eval_func(c0_4, x), 'k')
    plt.xlim(0.0, 0.625)
    plt.ylim(0, 1.1)
 ```
 %% Cell type:code id: tags:
 ``` python
 plt.plot(x, eval_func(c0_1, x)/0.9)
 plt.plot(x, eval_func(c0_2, x)/0.8)
 plt.plot(x, eval_func(c0_3, x)/0.7)
 plt.xlim(0.0, 0.125)
 plt.ylim(0, 1.1)
 ```
 %% Cell type:markdown id: tags:
+### Comparison with Matlab
+%% Cell type:code id: tags:
+``` python
+ML_9 = np.genfromtxt('Matlab_workspaces/ML_9_1.csv', delimiter=',')
+ML_1 = np.genfromtxt('Matlab_workspaces/ML_1_9.csv', delimiter=',')
+plt.plot(ML_9[0, :], ML_9[-1, :])
+plt.plot(x, eval_func(c0_9, x))
+plt.xlim(0, 0.5)
+# plt.show()
+plt.plot(ML_1[0, :], ML_1[-1, :])
+plt.plot(x, eval_func(c0_1, x))
+plt.xlim(0, 0.5)
+plt.show()
+```
+%% Cell type:markdown id: tags:
 ## Radial diffusion equation
 %% Cell type:code id: tags:
 ``` python
 mesh = df.UnitIntervalMesh(1000)
 dt = 0.001
 F = df.FunctionSpace(mesh, 'CG', 1)
 c0 = df.Function(F)
 c0.interpolate(df.Expression('x[0]<0.5 && x[0]>0.2 ? 1:0', degree=1))
 q = df.TestFunction(F)
 c = df.Function(F)
 X = df.SpatialCoordinate(mesh)
 g = df.Expression('.00', degree=1)
 u_D = df.Expression('1', degree=1)
 def boundary(x, on_boundary):
    return on_boundary
 bc = df.DirichletBC(F, u_D, boundary)
 # Weak form spherical symmetry
 form = ((c-c0)/dt*q*X[0]*X[0] +
        X[0]*X[0]*df.inner(df.grad(c), df.grad(q))-
        c.dx(0)*2*X[0]*q) * df.dx
 # Weak form 1D
 # form = ((c-c0)/dt* q + df.inner(df.grad(c), df.grad(q))) * df.dx
 # Weak form 1D with .dx(0) notation for derivative in 1st direction.
 # form = ((c-c0)/dt*q + c.dx(0)*q.dx(0)) * df.dx
 t = 0
 # Solve in time
 for i in range(60):
    print(np.sum([x*x*c0([x]) for x in np.linspace(0, 1, 1000)]))
    df.solve(form == 0, c)
    df.assign(c0, c)
    t += dt
    plt.plot(np.linspace(0, 1, 1000), [c0([x]) for x in np.linspace(0, 1, 1000)])
 ```
 %% Cell type:code id: tags:
 ``` python
 ```

--- a/prob_laplace.m
+++ b/prob_laplace.m
@@ -2,7 +2,24 @@
 % 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.02, 0.05, 0.1, 1, 50];
+% t = [0, 0.0001, 0.02, 0.05, 0.1, 1, 50];
+t = linspace(0.0, 0.01, 11);
+% t=0.01;
+% t = logspace(-4, -2, 20);
+%% Comparison with Python
+T1 = Ternary_model(0, 'FRAP', {-0.1, b(7/3, 10^-12), 0.5, 0.4,...
+                                     0, 1, 300, 7, 0, 'Constituent'},...
+                                     t, 1);
+T1.solve_tern_frap();
+csvwrite('ML_9_1.csv', [T1.x; T1.sol]);
+T2 = Ternary_model(0, 'FRAP', {-0.1, b(7/3, 10^-12), 0.9, 0,...
+                                     -0.895, 1, 300, 7, 0, 'Constituent'},...
+                                     t, 1);
+T2.x = T1.x;
+T2.solve_tern_frap();
+csvwrite('ML_1_9.csv', [T2.x; T2.sol]);
 %% Try out different precisions
 % prec = [0.5, 0.5, 1, 2, 3.5, 5];
 fac = [1, 5, 10];%, 10, 100];