Clear Filters
Clear Filters

Why am I receiving the error?

3 views (last 30 days)
Rahul
Rahul on 4 Dec 2023
Commented: Walter Roberson on 4 Dec 2023
Hello,
I'm receiving an error in my code when run. It shows
Error using pdepe
Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 4.
I tried a lot to debug the error but is unsuccessful.
Hence, I request your generous help in this regard.
My code is slightly lengthy. Thus I'm not sure if I should insert it here or attach it along with this.
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;
% --------------------------------------------------------------
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
disp(f);
disp(s);
% --------------------------------------------------------------
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)];
% --------------------------------------------------------------
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];
%---------------------------------------------------------------

Sign in to answer this question.

Accepted Answer

Walter Roberson
Walter Roberson on 4 Dec 2023
Ran in:
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()
Name Size Bytes Class Attributes DuDx 4x1 32 double H0 1x1 8 double global S0 1x1 8 double global c 4x1 32 double f 4x1 32 double length 1x1 8 double global s 7x1 56 double s3 4x1 32 double x 1x1 8 double -0.0101 -0.0101 -0.0006 0.0000 40.4312 10.5000 -0.7959 -0.7959 -0.7959 -0.7999 0
Error using pdepe
Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 4.

Error in solution>pde2fshear_v4Perturbed_randomWalk (line 106)
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
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
  2 Comments
Rahul
Rahul on 4 Dec 2023
Hi Walter Roberson,
Thanks a lot. So, what should be the correction. Can u please suggest the correction for the same?
With regards
Walter Roberson
Walter Roberson on 4 Dec 2023
Are you expecting your s3 to be a 4 x 1 vector? If not then you need to think about why you are computing it using the 4 x 1 DuDx .

Sign in to comment.

More Answers (0)

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!