From ae8fe0e3e5e3482dbbe056a7680fcacc2dbc9669 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Sun, 14 Mar 2021 17:11:36 +0100
Subject: [PATCH] Introduce constant flux case with Dirichlet BC. Works in
 principle.

---
 @Ternary_model/Ternary_model.m   |  6 ++++--
 @Ternary_model/solve_tern_frap.m | 30 +++++++++++++++++++-----------
 flory_hugg_pde.m                 |  5 ++++-
 prob_laplace.m                   |  5 +++++
 4 files changed, 32 insertions(+), 14 deletions(-)

diff --git a/@Ternary_model/Ternary_model.m b/@Ternary_model/Ternary_model.m
index 48c6438..a8eafaf 100644
--- a/@Ternary_model/Ternary_model.m
+++ b/@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);
diff --git a/@Ternary_model/solve_tern_frap.m b/@Ternary_model/solve_tern_frap.m
index 2b17c88..5c20dff 100644
--- a/@Ternary_model/solve_tern_frap.m
+++ b/@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
diff --git a/flory_hugg_pde.m b/flory_hugg_pde.m
index eb37e12..b32428a 100644
--- a/flory_hugg_pde.m
+++ b/flory_hugg_pde.m
@@ -1,5 +1,5 @@
 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
diff --git a/prob_laplace.m b/prob_laplace.m
index a9b795b..e69edc7 100644
--- a/prob_laplace.m
+++ b/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
-- 
GitLab