DuDx is 4 x 1. It is used to calculate s3, so s3 is 4 x 1. s3 is used as the third component of s, so s ends up length 1 + 1 + 3 + 1 = 7 when it is expected to be length 4.
pde2fshear_v4Perturbed_randomWalk()
function pde2fshear_v4Perturbed_randomWalk
global chi0; % declare global variables
global D0;
global chi1;
global D1;
global alpha_chi;
global alpha_D;
global H0;
global S0;
global H;
global S;
global data;
global track;
global xstep;
global tstep;
global count;
global chi_growth;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global sigma_turb;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
%define values for constants
data.constant.critgradpressure = 1.2;
data.constant.critgraddensity = 1.0;
count = 1;
chi0 = 0.5;
D0 = 0.5;
chi1 = 5.0;
D1 = 5.0;
alpha_chi = 0.1;
alpha_D = 0.1;
H0 = 27;
S0 = 21;
track = 1;
chi_growth=20;
lambda_suppress=0.5 ;
sigma_turb=0.5;
chi_ano=10;
D_ano=10;
% For random walk model
gamma_nonlin = 5;
alpha_nonlin = 1;
%position and time grids information
xstep = 100;
tstep = 100;
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 10;
%tmax = 1;
m = 0; %define type of equation to solve
%Preallocate vectors for speed improvement
grad_u1 = zeros(tstep,xstep);
grad_u2 = zeros(tstep,xstep);
grad_u3 = zeros(tstep,xstep);
grad_u4 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
curve_u4 = zeros(tstep,xstep);
flowshear_p = zeros(tstep,xstep);
flowshear_n = zeros(tstep,xstep);
Q = zeros(tstep,xstep);
Q0 = zeros(tstep,xstep);
Q1 = zeros(tstep,xstep);
neo_p = zeros(tstep,xstep);
ano_p = zeros(tstep,xstep);
Gam = zeros(tstep,xstep);
Gam0 = zeros(tstep,xstep);
Gam1 = zeros(tstep,xstep);
neo_n = zeros(tstep,xstep);
ano_n = zeros(tstep,xstep);
%x = linspace(xmin,xmax/2,xstep/5);
x = linspace(xmin,xmax,xstep);
t = linspace(tmin,tmax,tstep);
data.variable.x = x;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% Extract the first solution component as u1 = pressure
% second solution component as u2 = density
% third solution component as u3 = turbulence intensity
% fourth solution component as u4 = intensity_diff
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
u4 = sol(:,:,4);
%grad_u = gradient(u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
grad_u1(j,i) = (u1(j,2)-u1(j,1))/(x(2)-x(1));
grad_u2(j,i) = (u2(j,2)-u2(j,1))/(x(2)-x(1));
grad_u3(j,i) = (u3(j,2)-u3(j,1))/(x(2)-x(1));
grad_u4(j,i) = (u4(j,2)-u4(j,1))/(x(2)-x(1));
elseif i == xstep/5
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
grad_u1(j,i) = (u1(j,i+1)-u1(j,i))/(x(i+1)-x(i));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i))/(x(i+1)-x(i));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i))/(x(i+1)-x(i));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i))/(x(i+1)-x(i));
elseif i == xstep
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
end
end
end
for i=1:tstep
for j=1:xstep
v_e = -grad_u1(i,j)*grad_u2(i,j)/u2(i,j)^2; % -g_p*g_n/n^2
flowshear_p(i,j) = 1+ alpha_chi*v_e^2;
flowshear_n(i,j) = 1+ alpha_D*v_e^2;
if abs(grad_u1(i,j)) < abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
Q(i,j) = -grad_u1(i,j)*chi0;
Q0(i,j) = Q(i,j);
Q1(i,j) = 0;
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = 0;
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
elseif abs(grad_u1(i,j)) >= abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
else
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
Gam(i,j) = -grad_u2(i,j)*(D0 + D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j));
Gam0(i,j) = -grad_u2(i,j)*D0;
Gam1(i,j) = -grad_u2(i,j)*D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j);
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = D_ano*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j);
end
end
end
%curve_u = gradient(grad_u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
curve_u1(j,i) = (grad_u1(j,2)-grad_u1(j,1))/(x(2)-x(1));
curve_u2(j,i) = (grad_u2(j,2)-grad_u2(j,1))/(x(2)-x(1));
curve_u3(j,i) = (grad_u3(j,2)-grad_u3(j,1))/(x(2)-x(1));
curve_u4(j,i) = (grad_u4(j,2)-grad_u4(j,1))/(x(2)-x(1));
elseif i == xstep/5
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i))/(x(i+1)-x(i));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i))/(x(i+1)-x(i));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i))/(x(i+1)-x(i));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i))/(x(i+1)-x(i));
elseif i == xstep
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
end
end
end
%to save parameters and variables
data.constant.chi0 = chi0;
data.constant.D0 = D0;
data.constant.chi1 = chi1;
data.constant.D1 = D1;
data.constant.alphachi = alpha_chi;
data.constant.alphaD = alpha_D;
data.constant.H0 = H0;
data.constant.S0 = S0;
data.variable.pressure = u1;
data.variable.density = u2;
data.variable.turbulence = u3;
data.variable.intensity_diff = u4;
data.variable.drift_velFluct = drift_velFluct;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradintensity = grad_u3;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curveturbulence = curve_u3;
data.variable.x = x;
data.variable.t = t;
data.variable.Q = Q;
data.variable.Gamma = Gam;
data.variable.Q0 = Q0;
data.variable.Gamma0 = Gam0;
data.variable.neo_P = neo_p;
data.variable.neo_n = neo_n;
data.variable.Q1 = Q1;
data.variable.Gam1 = Gam1;
data.variable.ano_p = ano_p;
data.variable.ano_n = ano_n;
data.variable.heatsource = H;
data.variable.particlesource = S;
data.variable.wexb_p = flowshear_p;
data.variable.wexb_n = flowshear_n;
data.control.xgrid = xstep;
data.control.tgrid = tstep;
data.control.xmin = xmin;
data.control.xmax = xmax;
data.control.tmin = tmin;
data.control.tmax = tmax;
end
% --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global track;
global xstep;
global alpha_chi;
global alpha_D;
global chi_growth; % total growth rate
global length;
global theta_heaviside1;
global lambda_suppress;
global v_e;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
%lf FFf F Fb vfDdength = 0.01 or 1
length=1;
group_vel = 1;
D_0 = 1;
c = [1;1;1;1];
v_e = -DuDx(1)*DuDx(2)/u(2)^2; % -g_p*g_n/n^2
flowshear_p = 1+ alpha_chi*v_e^2;
flowshear_n = 1+ alpha_D*v_e^2;
term1 = abs(DuDx(1))-data.constant.critgradpressure;
intensity_diff = D_0*u(3)^alpha_nonlin;
drift_velFluct = DuDx(4);
drift_vel = group_vel + drift_velFluct;
% Implementing Heaviside function for H-mode in p, n and I equations
if term1 > 0
theta_heaviside1=1;
else
theta_heaviside1=0;
end
%Turbulence intensity Equation for random walk
s3 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(3)^alpha_nonlin)*u(3)-drift_vel*DuDx ;
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2; s3;0];
f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;intensity_diff*u(3); intensity_diff].*DuDx; % flux term for random walk model
whos c f s DuDx H0 length S0 x s3
disp(f);
disp(s);
end
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01; 0.1*exp(-100*(x-1)^2)];
u0 = [0.4*(1-x^2); 0.4*(1-x^2); 0.4*(1-x^2);0.5*exp(-100*(x-1)^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = [0; 0; 0;0];
ql = [1; 1; 1;1];
pr = [ur(1); ur(2);0; ur(4)];
qr = [eps; 0.1; 1;1];
% qr = [0; 0.1; 1];
%qr = [0; eps; 1];
%---------------------------------------------------------------
end
