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

17 views (last 30 days)
Widitha Samrakoon on 18 Feb 2016
Commented: Torsten on 22 Feb 2016
I get this error in writing the code.... The code I wrote is as follows
function pdeop global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a; m = 0; x = linspace(0,1.5,31); t = linspace(0,3600,3601); sol = pdepe(m,@pdeoppde,@pdeopic,@pdeopbc,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2); u3=sol(:,:,3); u4=sol(:,:,4);
% -------------------------------------------------------------- function [c,f,s] = pdeoppde(x,t,u,DuDx) c = [A*void*(1-hl)*Rhog;A*void*hl*Rhol;A*void*(1-hl)*Rhog;A*void*hl*Rhol]; f = [-G;L;-G;L] .*u; Cs=Rhog*Kos*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff); Ct=Rhog*Kot*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff); ysstar=-4.2913*u(2)^4+11.791*u(2)^3-12.01*u(2)^2+5.508*u(2); ytstar=0.0023*exp(6.0058*u(4)); M=Cs*(ysstar-u(1)); N=Ct*(ytstar-u(3)); s = [M;-M;N;-N]; end % -------------------------------------------------------------- function u0 = pdeopic(x) u0 = [0;0.5;0;0.5]; end % -------------------------------------------------------------- function[pl,ql,pr,qr] = pdeopbc(xl,ul,xr,ur,t) pl = [ul(1);0;ul(3);0]; ql = [0;0;0;0]; pr = [ur(1)-0.5;0;ur(3)-0.5;0]; qr = [0;0;0;0]; end end
Can anyone please suggest me a solution?
Widitha Samrakoon on 19 Feb 2016
Edited: Widitha Samrakoon on 19 Feb 2016
My code contains two functions named paracal.m which calculated parameters to be used in the pde system and then the main function named pdeop.m. First paracal.m is run and then pdeop.m is run. First I will paste the code of paracal.m and then pdeop.m. The error occurs when running pdeop.m.
function paracalc
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a;
A=pi*0.0254^2*0.25;
a=1710;
L=7.2888*10^-5;
G=2.4296*10^-3;
Rhog=842.77;
Rhol=875;
void=0.86;
Usg=G/(Rhog*A);
Usl=L/(Rhol*A);
hl=(1.5688*Usl^(1/3))*(1+4.1435*Usg^0.05)*(1-exp(-88.9560*Usl^(0.4)))^(2/3);
Kos=(1278102.274+6643.8785*Usg^-0.5)^-1;
Kot=(1620456.645+10029.7854*Usg^-0.5)^-1;
aeff=a*(1-exp(-88.9560*Usl^0.4));
Ugeff=1.6444*Usg/(1-hl);
Uleff=1.6444*Usl/hl;
end
Now I will paste the code of pdeop.m which gives the error
function pdeop
global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a;
m = 0;
x = linspace(0,1.5,31);
t = linspace(0,3600,61);
sol = pdepe(m,@pdeoppde,@pdeopic,@pdeopbc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3=sol(:,:,3);
u4=sol(:,:,4);
function [c,f,s] = pdeoppde(x,t,u,DuDx)
c = [A*void*(1-hl)*Rhog;A*void*hl*Rhol;A*void*(1-hl)*Rhog;A*void*hl*Rhol];
f =[-G;L;-G;L].*u+[-0.001;0.001;-0.001;0.001].*DuDx;
Cs=Rhog*Kos*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff);
Ct=Rhog*Kot*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff);
ysstar=-4.2913*u(2)^4+11.791*u(2)^3-12.01*u(2)^2+5.508*u(2);
ytstar=0.0023*exp(6.0058*u(4));
M=Cs*(ysstar-u(1));
N=Ct*(ytstar-u(3));
s = [M;-M;N;-N];
function u0 = pdeopic(x)
u0 = [0;0.5;0;0.5];
end
function[pl,ql,pr,qr] = pdeopbc(xl,ul,xr,ur,t)
pl = [ul(1);0;ul(3);0];
ql = [0;0;0;0];
pr = [0;ur(2)-0.5;0;ur(4)-0.5];
qr = [0;0;0;0];
end
end
Please note that my system of PDEs has four elliptical equations but the problem is all four of them do not have second order spatial derivatives. Therefore in my system, the flux vector doesn't have a DuDx term but U term only. By a search in the internet I got to know that pdepe function cannot work with such equations.However in order to get rid of this, as suggested in an answer to a similar question, I included a second order spatial derivative to all four equations which you can see as a separate column vector in the flux vector as [-0.001;0.001,-0.001,0.001]. I hoped this would solve the problem. But the error message keeps popping out. I am really helpless. Please help.
I am attaching a pdf file of my PDE system with the initial conditions and boundary conditions.

Bill Greene on 20 Feb 2016
I see two problems.
First, your boundary conditions are incorrectly defined. They should be:
pl = [ul(1);ul(2);ul(3);ul(4)];
pr = [ur(1);ur(2)-0.5;ur(3);ur(4)-0.5];
(ql and qr are correct). All four of your equations are hyperbolic which pdepe is not designed to handle. But, you have added some artificial diffusion (the small second derivative terms). That is a good approach but because of the negative signs on some of the coefficients you have de-stabilized the system instead stabilizing it. Try this for the vector of diffusion coefficients, instead:
eps=.001;
[eps eps eps eps]'
##### 2 CommentsShowHide 1 older comment
Bill Greene on 21 Feb 2016
OK, I took a guess at what your intention was with respect to the boundary conditions. pdepe requires boundary conditions at both ends and p=q=0 is not acceptable. I suggest p=0 and q=1 at the ends where you really don't want a BC. This will set your "artificial" flux to zero there but have minimal impact on your "real" PDE.

Torsten on 22 Feb 2016
As Bill already pointed out, pdepe is not designed to solve PDEs without 2nd derivative terms.
You will have to discretize the spatial first derivatives (e.g. by a first-order upwind scheme) and solve the resulting system of ordinary differential equations using ODE15S, e.g. . Look up "method of lines" for more details.
Best wishes
Torsten.
##### 2 CommentsShowHide 1 older comment
Torsten on 22 Feb 2016
If you have the NAG Toolbox, use
Best wishes
Torsten.