Clear Filters
Clear Filters

Why is the error in the attached code?

3 views (last 30 days)
Hi,
I'm getting an error in the attached code which I'm unable to decode.
Hence need your advice in the same.
Attached herewith the code.
% solve 3-F bifurcaiton model
function pde2fshear_v4Perturbed_nonlinear
global x;
global X;
global t;
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 gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global xmax;
global xmin;
%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.05;
chi_ano=10;
D_ano=10;
% For nonlinear model
gamma_nonlin = 50;
alpha_nonlin = 1;
D_0 = 10;
%position and time grids information
xstep = 100;
tstep = 1000;
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 50;
%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);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = 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);
safetyFact=zeros(tstep,xstep);
elec_diaVel=zeros(tstep,xstep);
grad_length=zeros(tstep,xstep);
mag_shear=zeros(tstep,xstep);
theta_deg=zeros(tstep,xstep);
elec_temp=zeros(tstep,xstep);
T_magFld=zeros(tstep,xstep);
P_magFld=zeros(tstep,xstep);
R0=zeros(tstep,xstep);
I_p=zeros(tstep,xstep);
jb=zeros(tstep,xstep);
B0=zeros(tstep,xstep);
%define first inner half of the plasma
%x = linspace(xmin,xmax/2,xstep/5);
x = linspace(xmin,xmax,xstep);
t = linspace(tmin,tmax,tstep);
%Add smaller grid size near plasma edge
%for i=(xstep/5)+1:xstep
% x(i) = x(i-1)+(xmax/2)/(xstep-(xstep/5));
%end
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
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
%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;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
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;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
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);
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(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=intensity_diff;
data.variable.drift_velFluct=drift_velFluct;
data.variable.drift_vel=drift_vel;
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;
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u1(end,:),'.');
axis ([0 1 0 20]);
ylabel('p')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u1(end,:),'.');
axis([0 1 -60 0]);
ylabel('p\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u1(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('p\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u1(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|p\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u2(end,:),'.');
axis ([0 1 0 20]);
ylabel('n')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u2(end,:),'.');
axis([0 1 -60 0]);
ylabel('n\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u2(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('n\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u2(end,:)),Gam(end,:),'.');
%plot(nn_pre,Gam_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('\Gamma')
xlabel('|n\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u3(end,:),'.');
axis ([0 1 0 20]);
ylabel('I')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u3(end,:),'.');
axis([0 1 -60 0]);
ylabel('I\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u3(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('I\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u3(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|I\prime|')
%title('Bifurcation diagram')
grid on
% --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global x;
global t;
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
global tstep;
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;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global X;
global xmax;
global xmin;
%lf FFf F Fb vfDdength = 0.01 or 1
x_count=1;
length=1;
jb0=1;
jd0=1;
x0=1;
grad_length=DuDx(2);
R0=linspace(1,5,100);
B0=1;
c = [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;
X = linspace(xmin,xmax,xstep-1);
term1 = abs(DuDx(1))-data.constant.critgradpressure;
drift_velFluct = D_0*DuDx(3)^alpha_nonlin;
intensity_diff = D_0*u(3)^alpha_nonlin;
% Implementing Heaviside function for H-mode in p, n and I equations
if term1 > 0
theta_heaviside1=1;
else
theta_heaviside1=0;
end
% Determining group velocity of turbulence intensity
if x_count > xstep-1
x_count=1;
end
if x_count <= xstep-1
X(x_count)=x;
num=15;
if x_count==1
jb=jb0*DuDx(1);
jd=jd0.*((1-X(x_count)-x0).^2).^num;
j_tot=jb(x_count)+jd;
for n=1:5
I_p=trapz(x(1:n),jb(x(1:n))+jd(x(1:n)));
end
end
T_magFld=B0./R0;
P_magFld=B0.*I_p./x;
safetyFact=(x./R0).*T_magFld./P_magFld;
if x_count == 1
grad_safetyFact(x_count) = (safetyFact(2)-safetyFact(1))/(x(2)-x(1));
elseif x_count == xstep/5
grad_safetyFact(x_count) = (safetyFact(x_count)-safetyFact(x_count-1))/(x(x_count)-x(x_count-1));
else
grad_safetyFact(x_count) = (safetyFact(x_count+1)-safetyFact(x_count-1))/(x(xstep+1)-x(x_count-1));
end
mag_shear=(x/safetyFact)*grad_safetyFact;
end
elec_diaVel = -elec_temp/1.6*10^(-19)*T_magFld*grad_length;
group_vel = - (elec_diaVel*(2*grad_length/mag_shear*R0)*sin(theta_deg));
%Turbulence intensity Equation for nonlinear turbulence intensity
s3 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(3)^alpha_nonlin)*u(3)+group_vel*u(3)-(D_0*(1+alpha_nonlin)*alpha_nonlin*u(3)^(alpha_nonlin-1)*u(3)^2);
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2;s3];
f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;-(intensity_diff-D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin)].*DuDx; % flux term for nonlinear model
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01;0.1*exp(-100*(x-1)^2)];
u0 = [0.1*(1-x^2); 0.1*(1-x^2); 0.5*exp(-100*(x-1)^2)];
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = [0; 0; 0];
ql = [1; 1; 1];
pr = [ur(1); ur(2);0];
qr = [0.01; 0.1; 1];
%---------------------------------------------------------------

Accepted Answer

Les Beckham
Les Beckham on 31 Jan 2024
After adding end to all of the functions so that this would even run, you get an dimension error (see below).
You are trying assign a 1x100 array x to a single element of X. You need to figure out what you really meant to do there.
pde2fshear_v4Perturbed_nonlinear
Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.
Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.
Name Size Bytes Class Attributes X 1x99 792 double global x 1x100 800 double global x_count 1x1 8 double
Unable to perform assignment because the left and right sides have a different number of elements.

Error in solution>pdex2pde (line 552)
X(x_count)=x;

Error in pdepe (line 242)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

Error in solution>pde2fshear_v4Perturbed_nonlinear (line 143)
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% solve 3-F bifurcaiton model
function pde2fshear_v4Perturbed_nonlinear
global x;
global X;
global t;
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 gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global xmax;
global xmin;
%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.05;
chi_ano=10;
D_ano=10;
% For nonlinear model
gamma_nonlin = 50;
alpha_nonlin = 1;
D_0 = 10;
%position and time grids information
xstep = 100;
tstep = 1000;
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 50;
%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);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = 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);
safetyFact=zeros(tstep,xstep);
elec_diaVel=zeros(tstep,xstep);
grad_length=zeros(tstep,xstep);
mag_shear=zeros(tstep,xstep);
theta_deg=zeros(tstep,xstep);
elec_temp=zeros(tstep,xstep);
T_magFld=zeros(tstep,xstep);
P_magFld=zeros(tstep,xstep);
R0=zeros(tstep,xstep);
I_p=zeros(tstep,xstep);
jb=zeros(tstep,xstep);
B0=zeros(tstep,xstep);
%define first inner half of the plasma
%x = linspace(xmin,xmax/2,xstep/5);
x = linspace(xmin,xmax,xstep);
t = linspace(tmin,tmax,tstep);
%Add smaller grid size near plasma edge
%for i=(xstep/5)+1:xstep
% x(i) = x(i-1)+(xmax/2)/(xstep-(xstep/5));
%end
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
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
%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;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
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;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
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);
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(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=intensity_diff;
data.variable.drift_velFluct=drift_velFluct;
data.variable.drift_vel=drift_vel;
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;
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u1(end,:),'.');
axis ([0 1 0 20]);
ylabel('p')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u1(end,:),'.');
axis([0 1 -60 0]);
ylabel('p\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u1(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('p\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u1(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|p\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u2(end,:),'.');
axis ([0 1 0 20]);
ylabel('n')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u2(end,:),'.');
axis([0 1 -60 0]);
ylabel('n\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u2(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('n\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u2(end,:)),Gam(end,:),'.');
%plot(nn_pre,Gam_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('\Gamma')
xlabel('|n\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u3(end,:),'.');
axis ([0 1 0 20]);
ylabel('I')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u3(end,:),'.');
axis([0 1 -60 0]);
ylabel('I\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u3(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('I\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u3(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|I\prime|')
%title('Bifurcation diagram')
grid on
end
% --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global x;
global t;
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
global tstep;
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;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global X;
global xmax;
global xmin;
%lf FFf F Fb vfDdength = 0.01 or 1
x_count=1;
length=1;
jb0=1;
jd0=1;
x0=1;
grad_length=DuDx(2);
R0=linspace(1,5,100);
B0=1;
c = [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;
X = linspace(xmin,xmax,xstep-1);
term1 = abs(DuDx(1))-data.constant.critgradpressure;
drift_velFluct = D_0*DuDx(3)^alpha_nonlin;
intensity_diff = D_0*u(3)^alpha_nonlin;
% Implementing Heaviside function for H-mode in p, n and I equations
if term1 > 0
theta_heaviside1=1;
else
theta_heaviside1=0;
end
% Determining group velocity of turbulence intensity
if x_count > xstep-1
x_count=1;
end
if x_count <= xstep-1
whos X x x_count
X(x_count)=x;
num=15;
if x_count==1
jb=jb0*DuDx(1);
jd=jd0.*((1-X(x_count)-x0).^2).^num;
j_tot=jb(x_count)+jd;
for n=1:5
I_p=trapz(x(1:n),jb(x(1:n))+jd(x(1:n)));
end
end
T_magFld=B0./R0;
P_magFld=B0.*I_p./x;
safetyFact=(x./R0).*T_magFld./P_magFld;
if x_count == 1
grad_safetyFact(x_count) = (safetyFact(2)-safetyFact(1))/(x(2)-x(1));
elseif x_count == xstep/5
grad_safetyFact(x_count) = (safetyFact(x_count)-safetyFact(x_count-1))/(x(x_count)-x(x_count-1));
else
grad_safetyFact(x_count) = (safetyFact(x_count+1)-safetyFact(x_count-1))/(x(xstep+1)-x(x_count-1));
end
mag_shear=(x/safetyFact)*grad_safetyFact;
end
elec_diaVel = -elec_temp/1.6*10^(-19)*T_magFld*grad_length;
group_vel = - (elec_diaVel*(2*grad_length/mag_shear*R0)*sin(theta_deg));
%Turbulence intensity Equation for nonlinear turbulence intensity
s3 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(3)^alpha_nonlin)*u(3)+group_vel*u(3)-(D_0*(1+alpha_nonlin)*alpha_nonlin*u(3)^(alpha_nonlin-1)*u(3)^2);
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2;s3];
f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;-(intensity_diff-D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin)].*DuDx; % flux term for nonlinear model
end
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01;0.1*exp(-100*(x-1)^2)];
u0 = [0.1*(1-x^2); 0.1*(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];
ql = [1; 1; 1];
pr = [ur(1); ur(2);0];
qr = [0.01; 0.1; 1];
end
%---------------------------------------------------------------

More Answers (0)

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Tags

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!