Issue with solving pdepe

Hello,

Wishing you all a Merry Christmas

I'm trying to solve a pdepe consisting of 6 equations which consists pressure, density, turbulence intensity, diffusion intensity, drift velocity, drift fluctuation intensity.

while trying to solve the same I get error as follows:

Error using pdepe

Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term

involving spatial derivative.

Hence, needs your valuable advice please.

Attached herewith the code (let me tell you that it is slightly long code. Kindly bear with me please for this long code.

Sincerely yours,

rc

pde2fshear_v4Perturbed_nonlinear()

Error using pdepe
Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

Error in solution>pde2fshear_v4Perturbed_nonlinear (line 112)
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

% solve 3-F bifurcaiton model 
function pde2fshear_v4Perturbed_nonlinear
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;
global D_0;


  
%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 = 1;
alpha_nonlin =0.9;





%position and time grids information
xstep = 100;
tstep = 1000; 
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);
grad_u5 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
curve_u4 = zeros(tstep,xstep);
curve_u5 = 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);


%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 = intensity_diffusion
% fourth solution component as u4 = Drift velocity intensity
% fifth solution component as u5 = fluct_drift_intensity
% sixth solution component as u6 = turbulence
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3); 
u4 = sol(:,:,4);
u5 = sol(:,:,5);
u6 = sol(:,:,6);
%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));
            grad_u5(j,i) = (u5(j,2)-u5(j,1))/(x(2)-x(1));
            grad_u6(j,i) = (u6(j,2)-u6(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));
            grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
            grad_u6(j,i) = (u6(j,i)-u6(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
            grad_u6(j,i) = (u6(j,i+1)-u6(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(j,i))/(x(i+1)-x(i));
            grad_u6(j,i) = (u6(j,i+1)-u6(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));
            grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
            grad_u6(j,i) = (u6(j,i)-u6(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
            grad_u6(j,i) = (u6(j,i+1)-u6(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*u5(i,j)^alpha_nonlin;
             drift_velFluct = grad_u(3);
             drift_vel = group_vel + drift_velFluct;
         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)*u5(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)*u5(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*u5(i,j)^alpha_nonlin;
            drift_velFluct = grad_u(3);
            drift_vel = group_vel + drift_velFluct;
       else
             Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(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)*u5(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)*u5(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)*u5(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)*u5(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*u5(i,j)^alpha_nonlin;
             drift_velFluct = grad_u(3);
             drift_vel = group_vel + drift_velFluct;
         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));
            curve_u5(j,i) = (grad_u5(j,2)-grad_u5(j,1))/(x(2)-x(1));
            curve_u6(j,i) = (grad_u6(j,2)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i+1)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i))/(x(i+1)-x(i));
            curve_u6(j,i) = (grad_u6(j,i+1)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i+1)-grad_u6(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.intensity_diffusion = u3;
data.variable.drift_velocity= u4;
data.variable.drift_velFluct= u5;
data.variable.turbulence= u6;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradintensity = grad_u5;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curveturbulence = curve_u5;
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 sigma_turb;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;


%lf FFf F  Fb vfDdength = 0.01 or 1
length=1;
D_0 = 1;
group_vel=1;


c = [1;1;0;0;0;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(6)^alpha_nonlin;
drift_velFluct = DuDx(3);
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

s6 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(6)^alpha_nonlin)*u(6)-drift_vel*DuDx(6) ;

 
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2;intensity_diff;0; group_vel;s6];


f = [chi0+chi1*u(6)/flowshear_p ; D0+D1*u(6)/flowshear_n ;0;drift_velFluct;intensity_diff;0].*DuDx;                 % flux term for nonlinear model
end

% --------------------------------------------------------------
function u0 = pdex2ic(x)
  u0 = [eps; eps; eps;eps;eps;eps];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)

pl = [0; 0; 0;0;0;0];
ql = [1; 1; 1;1;1;1];


pr = [ur(1); ur(2); ur(3);ur(4);ur(5);0];
qr = [0.01; 0.1; 1;1;1;1];

end

%---------------------------------------------------------------

[Bill Greene]

I think you are going to have to provide a complete, mathematical description of your problem-- PDE plus boundary conditions.

[Torsten]

Please include your equations together with initial and boundary conditions in a mathematical notation so that we can compare with your implementation for "pdepe".

Looking over your code, there are many reasons to believe that your problem is not suited to be solved with "pdepe".

[Rahul]

Hi,

plz find attached the set of equations

Also i found one mistake in the code that I corrected in the code below. However, it gives the same error when run as before.

% solve 3-F bifurcaiton model 
function pde2fshear_v4Perturbed_nonlinear
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;
global D_0;


  
%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 nonlinear model
gamma_nonlin = 1;
alpha_nonlin =0.5;


%position and time grids information
xstep = 100;
tstep = 1000; 
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);
grad_u5 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
curve_u4 = zeros(tstep,xstep);
curve_u5 = 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);


%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 = intensity_diffusion
% fourth solution component as u4 = Drift velocity intensity
% fifth solution component as u5 = fluct_drift_intensity
% sixth solution component as u6 = turbulence
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3); 
u4 = sol(:,:,4);
u5 = sol(:,:,5);
u6 = sol(:,:,6);
%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));
            grad_u5(j,i) = (u5(j,2)-u5(j,1))/(x(2)-x(1));
            grad_u6(j,i) = (u6(j,2)-u6(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));
            grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
            grad_u6(j,i) = (u6(j,i)-u6(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
            grad_u6(j,i) = (u6(j,i+1)-u6(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(j,i))/(x(i+1)-x(i));
            grad_u6(j,i) = (u6(j,i+1)-u6(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));
            grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
            grad_u6(j,i) = (u6(j,i)-u6(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));
            grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
            grad_u6(j,i) = (u6(j,i+1)-u6(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*u5(i,j)^alpha_nonlin;
             drift_velFluct = grad_u(3);
             drift_vel = group_vel + drift_velFluct;
         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)*u5(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)*u5(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*u5(i,j)^alpha_nonlin;
            drift_velFluct = grad_u(3);
            drift_vel = group_vel + drift_velFluct;
       else
             Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(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)*u5(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)*u5(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)*u5(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)*u5(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*u5(i,j)^alpha_nonlin;
             drift_velFluct = grad_u(3);
             drift_vel = group_vel + drift_velFluct;
         end
     end
 end


%for i=1:xstep
 %   a5(i) = sol_dirac(i) ; 
  %  data.variable.a5_term(i) = a5(i);
%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));
            curve_u5(j,i) = (grad_u5(j,2)-grad_u5(j,1))/(x(2)-x(1));
            curve_u6(j,i) = (grad_u6(j,2)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i+1)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i))/(x(i+1)-x(i));
            curve_u6(j,i) = (grad_u6(j,i+1)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i)-grad_u6(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));
            curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
            curve_u6(j,i) = (grad_u6(j,i+1)-grad_u6(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.intensity_diffusion = u3;
data.variable.drift_velocity= u4;
data.variable.drift_velFluct= u5;
data.variable.turbulence= u6;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradintensity = grad_u5;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curveturbulence = curve_u5;
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 sigma_turb;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;


%lf FFf F  Fb vfDdength = 0.01 or 1
length=1;
D_0 = 1;
group_vel=1;


c = [1;1;0;0;0;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(6)^alpha_nonlin;
drift_velFluct = DuDx(3);
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

s6 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(6)^alpha_nonlin)*u(6)-drift_vel*DuDx(6) ;

 
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2;intensity_diff;0; group_vel;s6];


f = [chi0+chi1*u(6)/flowshear_p ; D0+D1*u(6)/flowshear_n ;0;drift_velFluct;drift_velFluct;intensity_diff*u(6)].*DuDx;                 % flux term for nonlinear model

end


% --------------------------------------------------------------
function u0 = pdex2ic(x)
  u0 = [eps; eps; eps;eps;eps;eps];

% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)

pl = [0; 0; 0;0;0;0];
ql = [1; 1; 1;1;1;1];


pr = [ur(1); ur(2); ur(3);ur(4);ur(5);0];
qr = [0.01; 0.1; 1;1;1;1];


%---------------------------------------------------------------

[Torsten]

Only use the equations for p, n and I to be solved with pdepe - equations (3), (4) and (5) are algebraic relations that can be deduced from p, n and I and their spatial derivatives.

The equation for I doesn't seem to fit in the scheme for equations that are solvable with "pdepe". Or how do you arrive to write it in a form

c*dI/dt = d/dx(f*dI/dx) + s

?

The term -d/dx(V_I * I) will give another second-order derivative of I into your equation, not only the d^2/dx^2(D*I) term.

And the boundary condition setting at x = a look incorrect. Note that the boundary conditions are set via p and q in a relation of the form p + q*f = 0 where f is the flux function you set in "pdex2pde".

E.g. for your first equation pr = ur(1), qr = 0.01 means that you set

p + 0.01*(chi0+chi1*I/flowshear_p) = 0

[Rahul]

Hi,

Thanks a lot for your elaborated remarks which are indeed helpful.

1) Regarding your comment about c*dI/dt = d/dx(f*dI/dx) + s my answer is yes it is with the 1st and 2nd term represents the source while 3rd term representing the flux

2) Regarding your comment "The term -d/dx(V_I * I) will give another second-order derivative of I into your equation, not only the d^2/dx^2(D*I) term" then how should I proceed? Because this is the line which is presently giving me error which shows

Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 3.

3) I'm revisiting by boundary conditions.

4) But I have one query regarding my equation (4) i.e. the equation as below, this has to be solved using pdepe or has to be solved as a normal algegraic? kindly let me know.

V_I_Tilda=del/delX(D_0 * I^zeta)

[Torsten]

V_I_Tilda = D_0*zeta*I^(zeta-1)*dI/dx
d/dx(V_I * I) = 
V_I * dI/dx + d/dx(V_I_Tilda) * I = 
V_I * dI/dx + (D_0*zeta*I^(zeta-1)*d^2I/dx^2 + D_0*zeta*(zeta-1)*I^(zeta-2)*(dI/dx)^2 )*I

I don't see how the term D_0*zeta*I^zeta*d^2I/dx^2 could be accounted for in "pdepe".

Regarding your comment about c*dI/dt = d/dx(f*dI/dx) + s my answer is yes

So my answer is no.

[Rahul]

Hi Torsten,

Thanks a lot for your elaborated reply. I'm doing corrections. If any issue I'll get back to you.

Thank you once again.