Issue with solving pdepe
Mostra commenti meno recenti
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()
% 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
%---------------------------------------------------------------
7 Commenti
Bill Greene
il 25 Dic 2023
I think you are going to have to provide a complete, mathematical description of your problem-- PDE plus boundary conditions.
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".
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
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
il 27 Dic 2023
Risposte (0)
Categorie
Scopri di più su Mathematics in Centro assistenza e File Exchange
il 25 Dic 2023
il 28 Dic 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
