Introduce constant flux case with Dirichlet BC. Works in principle.

@Ternary_model/Ternary_model.m

@@ -23,6 +23,7 @@ classdef Ternary_model < handle
         e = 0.4;
         e_g0 = 0.16; % mobility spread. Also used in square mobility ansatz.
         ic_c % initial concentration inside droplet
+        j % parameter characterizing the flux in Fokker Planck eq.
         u_g0 = 0.2;
         system_size = 300;
         x0 = 1; % Center of Gauss initial condition
@@ -45,8 +46,9 @@ classdef Ternary_model < handle
                 T.system_size = params{7};
                 T.x0 = params{8};
                 T.v = params{9};
-                T.mode = params{10};
-                T.ic_c = params{11};
+                T.j = params{10};
+                T.mode = params{11};
+                T.ic_c = params{12};
             end
             T.create_mesh();
             T.phi_t = Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0, 0);

@Ternary_model/solve_tern_frap.m

@@ -3,9 +3,9 @@ function solve_tern_frap(T)
 % 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, T.ic_c);
-fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, T.u0, T.e);
+fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, T.u0, T.e, T.j);
 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.u0, T.e_g0, T.u_g0, T.v, T.j,...
                                           T.mode);
 opts = odeset('RelTol',1e-3,'AbsTol',1e-6);
 T.sol = pdepe(T.geometry, fh_pde, fh_ic, fh_bc, T.x, T.t, opts);
@@ -63,13 +63,21 @@ else
 end
 end
 
-function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0, e)
-    pl = 0;
-    ql = 1;
-%     % No flux
-    pr = 0;%ur - 0.01;
-    qr = 1;
-    % Dirichlet BC
-%     pr = ur-u0+e;
-%     qr = 0;
+function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0, e, j)
+    if j == 0
+        pl = 0;
+        ql = 1;
+        % No flux
+        pr = 0;%ur - 0.01;
+        qr = 1;
+        % Dirichlet BC
+%         pr = ur-u0+e;
+%         qr = 0;
+    % Dirichlet BC constant flux case
+    else
+        pr = ur;
+        qr = 0;
+        pl = ul;
+        ql = 0;
+    end
 end
\ No newline at end of file

flory_hugg_pde.m

 function [c, f, s, g0, pt] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
-                                    e_g0, u_g0, v, mode)
+                                    e_g0, u_g0, v, j, mode)
 % Solve with full ternary model. Analytical derivatives.
 % pt ... phi_tot
 % gra_a ... analytical gradient of phi_tot
@@ -19,5 +19,8 @@ elseif strcmp(mode, 'Client')
     g0 = 0.5;
     chi_phi = -4.530864768482371;
     f = g0*(dudx+chi_phi*u*gra_a);
+elseif strcmp(mode, 'Const_flux')
+    g0 = Ternary_model.gamma0(x, a+t*v, b, e_g0, u_g0, v*t, e, u0);
+    f = g0.*(1-pt)./pt.*(pt.*dudx-u.*gra_a)-j*u/pt;
 end
 end
\ No newline at end of file

prob_laplace.m

@@ -337,6 +337,11 @@ run_jump_lengths(-5, -0.00756985349, 7.7/3, 1, 4, '773_30_01', 10^-6, 3, 0.5);
 run_jump_lengths(-5, -0.00756985349, 7.7/3, 1, 4, '773_30_10', 10^-6, 4, 0.5);
 
 % output is saved manually from displayed values to data/mean_cross_jl_ML.csv
+%% Constant flux time course
+t = linspace(0, 5, 300);
+T_mov = Ternary_model(0, 'Gauss', {-5, b(7/3, 10^-6), 0.5, e(7/3),...
+                   0, 1, 10, 7, 0, 1, 'Const_flux', 0}, t, 0.2);
+T_mov.solve_tern_frap();
 %% %%%%%%%%%%%%%%%%%%% MOVING BOUNDARY TIME COURSE %%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % t = linspace(0, 10, 9000); % Entropy