Skip to content
Snippets Groups Projects
Commit 5f1b99c3 authored by Lars Hubatsch's avatar Lars Hubatsch
Browse files

WIP: introducing moving boundary. v=0 works as expected.

parent c1e040c0
No related branches found
No related tags found
No related merge requests found
......@@ -5,7 +5,7 @@ classdef Ternary_model < handle
% conditions are introduced. Integration of model via pdepe.
properties
pa = '~/ownCloud/Dropbox_Lars_Christoph/DropletFRAP/Latex/Figures/FigModel/';
pa = '~/Desktop/DropletFRAP/Latex/Figures/FigModel/';
geometry = 2; % 2 ... Radial, 1 ... Cylindrical, 0 ... Slab
ic = 'FRAP'
precision = 5;
......@@ -24,6 +24,7 @@ classdef Ternary_model < handle
system_size = 300;
x0 = 1; % Center of Gauss initial condition
real_params;
v; % velocity of boundary
end
methods
function T = Ternary_model(geometry, ic, params, t, prec)
......@@ -39,12 +40,11 @@ classdef Ternary_model < handle
T.e_g0 = params(5);
T.u_g0 = params(6);
T.system_size = params(7);
if strcmp(ic, 'Gauss')
T.x0 = params(8);
end
T.x0 = params(8);
T.v = params(9);
end
T.create_mesh();
T.phi_t = Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0);
T.phi_t = Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0, 0);
T.calc_real_params;
end
create_mesh(T)
......@@ -54,7 +54,7 @@ classdef Ternary_model < handle
calc_real_params(T);
end
methods(Static)
pt = phi_tot(x, a, b, e, u0)
pt = phi_tot(x, a, b, e, u0, vt)
g = gradient_analytical(x, a, b, e)
g0 = gamma0(x, a, b, e_g0, u_g0)
sp = spacing(x, a, b, e_g0, u_g0)
......
function calc_real_params(T)
% Calculate important real-life parameters, such as D_in, D_out,
% partitioning etc..
T.real_params.phi_in = Ternary_model.phi_tot(-1000, T.a, T.b, T.e, T.u0);
T.real_params.phi_out = Ternary_model.phi_tot(1000, T.a, T.b, T.e, T.u0);
T.real_params.phi_in = Ternary_model.phi_tot(-1000, T.a, T.b, T.e, T.u0, 0);
T.real_params.phi_out = Ternary_model.phi_tot(1000, T.a, T.b, T.e, T.u0, 0);
T.real_params.Gamma_in = Ternary_model.gamma0(-1000, T.a, T.b, T.e_g0,...
T.u_g0);
T.real_params.Gamma_out = Ternary_model.gamma0(1000, T.a, T.b, T.e_g0,...
......
......@@ -8,11 +8,11 @@ x = linspace(0, 0.75, 40);
x = [x, x(end)+0.01/T.precision/10];%/10
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);
phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0);
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);
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
......@@ -26,7 +26,7 @@ T.x = T.x(1:find(T.x-T.system_size>=0, 1)-1);
[~, ind] = min(abs(T.x-T.x0));
if T.x(ind+20)-T.x(ind) > 0.1
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)];
ceil((T.x(ind+20)-T.x(ind-20))*200)),...
T.x(ind+21:end)];
end
function p = phi_tot(x, a, b, e, u0)
p = e*tanh(-(x+a)/b)+u0;
function p = phi_tot(x, a, b, e, u0, vt)
p = e*tanh(-(x+a+vt)/b)+u0;
end
\ No newline at end of file
......@@ -4,18 +4,18 @@ function solve_tern_frap(T)
fh_ic = @(x) flory_ic(x, T.a, 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.u0, T.e_g0, T.u_g0, T.v);
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, T)
e_g0, u_g0, v)
% 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);
pt = Ternary_model.phi_tot(x, a, b, e, u0, t*v);
gra_a = Ternary_model.gradient_analytical(x, a, b, e);
g0 = Ternary_model.gamma0(x, a, b, e_g0, u_g0);
c = 1;
......
......@@ -7,8 +7,8 @@ t = [0, 0.1, 1, 9.9];
prec = [0.25, 0.5, 1, 2, 3.5, 5];
parfor i = 1:length(prec)
tic
T_prec(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
0.16, 0.2, 300, 7], t, prec(i));
T_prec(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, ...
e(7.7/3), 0, 1, 300, 7], t, prec(i));
T_prec(i).solve_tern_frap()
toc
end
......@@ -16,8 +16,8 @@ end
%% Same precisions, different starting positions
x0 = [-0.05, -0.03, -0.01, 0.01, 0.03, 0.05];
parfor i = 1:length(x0)
T_diff_x0(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3),...
0.16, 0.2, 300, 6+x0(i)], t, 0.5);
T_diff_x0(i) = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, ...
e(7.7/3), 0, 1, 300, 6+x0(i)], t, 0.5);
T_diff_x0(i).solve_tern_frap()
end
%%
......@@ -30,7 +30,7 @@ end
t = linspace(0, 10, 1000);
tic
T_mov = Ternary_model(0, 'Gauss', [-6, b(7/3, 10^-15), 0.5, e(7/3),...
0, 1, 300, 7], t, 0.5);
0, 1, 300, 7, 0], t, 0.02);
T_mov.solve_tern_frap()
toc
%%
......@@ -61,7 +61,8 @@ for i = 1:100%length(T_mov.t)
shg; pause();
end
%% Check integral of solution, should be mass conserving and sum to 1.
[~, ind] = min(abs(T_prec(4).x-50)); % integrate up to 50, to avoid right boundary
% integrate to 50, to avoid right boundary
[~, ind] = min(abs(T_prec(4).x-50));
bins = (T_prec(4).sol(:, 1:ind-1)+T_prec(4).sol(:, 2:ind))/2;
bin_size = diff(T_prec(4).x(1:ind));
sum(bins(4, :).*bin_size)
......@@ -87,21 +88,21 @@ chi = 7.7/3;
b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
e = @(chi) sqrt(3/8*(chi-2));
t = [0, 0.05, 0.1];
T1 = Ternary_model(0, 'Gauss', [-5, b(chi, nu), 0.5, e(chi), ...
0.16, 0.2, 10, a], t, 2);
% T1 = Ternary_model(0, 'Gauss', [-5, b(chi, nu), 0.5, e(chi), ...
% 0, 1, 10, a], t, 0.2);
x0 = 5.0:0.002:7.01;
%%
T = {};
tic
parfor i = 1:length(x0)
tic
b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
e = @(chi) sqrt(3/8*(chi-2));
T{i} = Ternary_model(0, 'Gauss', [-a, b(chi, nu), 0.5, e(chi), ...
0.16, 0.2, 10, x0(i)], t, 2);
0, 1, 10, x0(i)], t, 0.2);
T{i}.solve_tern_frap();
end
toc
save prob_laplace_X_7_7_short_2
end
% save prob_laplace_X_7_7_short_2
%%
ls = 0.0005:0.005:2;
p = nan(1, length(ls));
......@@ -114,13 +115,13 @@ plot(ls, p/N);
% plot(lx, xt);
% plot(ls, q/50000);
%%
int_prob(0.3, T, x0, 0.1)
%%
int_prob_simple(0.01, T, x0)
%%
for i = 1:2:length(T)
T{i}.plot_sim('plot', 3, 'red')
end
%%
int_prob(0.3, T, x0, 0.1)
%%
int_prob_simple(0.01, T, x0)
%% Normalization factor
tic
N = normalization(T, x0,0.1);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment