pdepe with cross elements in the differentiation "c"

Question

Francesca Turco 2023-8-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2013547-pdepe-with-cross-elements-in-the-differentiation-c

评论： carlos Hernando 2024-9-24

x=linspace(0,1,151); t=linspace(0,10,100);
[sol] = pdepe(m,@(x,t,U,dUdx)pdefun(x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPi,fPd,fPrad,Snbi,Spump, ...
                                    xx,cz0psm,neped,T0p,n0p,Pin0p,Tout,Vol,alfa,tshift), ...
                @(x)pdeic(x,Te,Lz,xx,cz0psm,T0p,n0p,Pin0p,alfa,tshift,Vol), ...
                @pdebc, ...
                x,t);
function [c,f,s] = pdefun (x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPi,fPd,fPrad,Snbi,Spump, ...
                           xx,cz0psm,neped,T0p,n0p,Pin0p,Tout,Vol,alfa,tshift)
intT0=1; intPin0=1; intn0=1; 
intnep=1; intcz=1; intV=1;
intDT=1; intDN=1; intTout=1;
intfPp=1; intfPi=1; intfPd=0.01;
Pr=L*1e34*interp1(Te,Lz,U(1),'pchip','extrap')*U(3)^2*intcz+ brem*0.00534*2*U(3)^2*sqrt(U(1))+ liner*0.01*2*U(3)^2;
Pa=0.01*U(3);
berror=U(1)*U(3)-intT0*intn0; inter=trapz(berror,xx);
%c=[1 1 1 1 1]';
%f=[intDT*dUdx(1) 0 intDN*dUdx(3) 0 0]'; 
%s=[fTpin*U(2)-fTpr*Pr-intTout*U(1) + alfa*Pa, ...  % Te
%   -intfPp*(berror)-intfPi*inter+fPrad*Pr, ...        % Pin   
%   Snbi*U(2)+intnep-Spump*U(3), ...                            % ne
%   1e34*interp1(Te,Lz,U(1),'pchip','extrap')*intcz, ...    % coefficient x Prad(Lz)
%   alfa*Pa]';  % Palpha 
c=[   1      0     0      0 0      % U1 =T
   -fPd*U(3) 1  -fPd*U(1) 0 0      % U2 =Pin     
      0      0     1      0 0      % U3 =ne
      0      0     0      1 0      % U4 =Lz
      0      0     0      0 1];    % U5 =Palpha
f=[intDT*dUdx(1) 0 0 0 0 
    0 0 0 0 0 
    0 0 intDN*dUdx(3) 0 0 
    0 0 0 0 0 
    0 0 0 0 0]; 
s=[fTpin*U(2)-fTpr*Pr-intTout*U(1)+alfa*Pa 0 0 0 0 
   0 -intfPp*(berror)-intfPi*inter+fPrad*Pr 0 0 0    
   0 0 Snbi*U(2)+intnep-Spump*U(3) 0 0 
   0 0 0 1e34*interp1(Te,Lz,U(1),'pchip','extrap')*intcz 0 
   0 0 0 0 alfa*Pa];   
end
function u0 = pdeic(x,Te,Lz,xx,cz0psm,T0p,n0p,Pin0p,alfa,tshift,Vol)
intT0=pchip(xx,T0p,x); intPin0=pchip(xx,Pin0p,x); intn0=pchip(xx,n0p,x); 
intcz=pchip(xx,cz0psm,x); intLz0=interp1(Te,Lz,intT0,'pchip','extrap'); intV=pchip(xx,Vol,x);
%u0=[intT0, intPin0, intn0, 1e34*intLz0*intcz, alfa*intT0^2*(0.78*intn0)^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(intT0+tshift)-0.0089152*(intT0+tshift).^2+0.000103742*(intT0+tshift).^3)]';
u0=[intT0 0 0 0 0
    0 intPin0 0 0 0 
    0 0 intn0 0 0 
    0 0 0 1e34*intLz0*intcz 0
    0 0 0 0 alfa*intT0^2*(0.78*intn0)^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(intT0+tshift)-0.0089152*(intT0+tshift).^2+0.000103742*(intT0+tshift).^3)];
u0=[intT0, intPin0, intn0, 1e34*intLz0*intcz, alfa*intT0^2*(0.78*intn0)^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(intT0+tshift)-0.0089152*(intT0+tshift).^2+0.000103742*(intT0+tshift).^3)]';
end
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) 
%pl = [0 0 0 0 0]';                   % Te, Pin, ne, coeff, Palpha
%ql = [1 1 1 1 1]'; 
%pr = [ur(1)-0.01 ur(2)-0.02 ur(3)-0.06 0 0]';  % Te, Pin, ne, coeff, Palpha
%qr = [1 1 1 1 1]';
pl = zeros(5,5);                   % Te, Pin, ne, coeff, Palpha
ql = ones(5,5); 
pr = [ur(1)-0.01 0 0 0 0
      0 ur(2)-0.02 0 0 0 
      0 0 ur(3)-0.06 0 0
      0 0 0 0 0
      0 0 0 0 0];  % Te, Pin, ne, coeff, Palpha
qr = ones(5,5);
end

Hello, I have used the pdepe to solve a system of 5 coupled PDEs, that works fine. I want to add two terms in the LHS dUdt that makes "c" be a combination of dUdt(1)+dUdt(2)+dUdt(3). Is this possible? Or does the LHS of pdepe "c" can only have non-zero diagonal values?

More basically: can pdepe accept matrices or just vector inputs in c, s, and f?

The code below is the main script that I made. It works fine with the c, f, s coefficients that are commented out. With the c, f, s expressed as matrices, it gives the error "Error using pdepe: Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 5."

I am not sure why this is. Am I expressing the matrices in the wrong way, or is this just not possible with pdepe?

2 个评论
显示无隐藏无

Walter Roberson 2023-8-28

Please edit your post to format the code.

Delete the current code, then press the > button in the CODE toolbar over the edit window. That will create a code insertion window. Copy the code from your editor and paste it in at that point, and it will be stored the same as is in your editor.

Francesca Turco 2023-8-28

Thanks, it reads much better now.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bill Greene 2023-8-29

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2013547-pdepe-with-cross-elements-in-the-differentiation-c#answer_1296291

在 MATLAB Online 中打开

If you are willing to consider an alternative to pdepe, I have written a PDE solver that mostly supports the same input syntax as pdepe but has several enhancements including the capability for a non-diagonal c-matrix. It can be downloaded here.

Here is a very simple two-equation PDE system that shows how you can define a non-diagonal c-matrix where the terms also depend on the dependent variables in the PDE.

function coupledMassExample
L=1;
n=21;
n2 = ceil(n/2);
x = linspace(0,L,n);
tf=.1;
nt=20;
t = linspace(0,tf,nt);
pdeFunc = @(x,t,u,DuDx) pdeDefinition(x,t,u,DuDx);
icFunc = @(x) icDefinition(x, L);
bcFunc = @(xl,ul,xr,ur,t) dirichletBC(xl,ul,xr,ur,t);
m=0;
sol = pde1dm(m, pdeFunc,icFunc,bcFunc,x,t);
u=sol(:,:,1);
v=sol(:,:,2);
figure; plot(t, u(:, n2), t, v(:, n2)); grid on;
xlabel 'time'
ylabel 'u,v'
title 'solution at center';
legend('u1', 'u2', 'Location', 'west')
figure; plot(x, u(end, :), x, v(end, :)); grid on;
xlabel 'x'
ylabel 'u,v'
title('solution at final time');
legend('u1', 'u2', 'Location', 'south')
fprintf('Solution: Time=%g, u_center=%g, v_center=%g\n', ...
  t(end), u(end,n2), v(end,n2));
end
function [c,f,s] = pdeDefinition(x,t,u,DuDx)
c=[3-.1*u(2) 1+.2*u(1)
   1+u(1) 2-.3*u(2)];
f = DuDx;
s=[0 0]';
end
function u0 = icDefinition(x,L)
u0 = [sin(pi*x/L); 3*sin(pi*x/L)];
end
function [pl,ql,pr,qr] = dirichletBC(xl,ul,xr,ur,t)
pl = ul;
ql = [0 0]';
pr = ur;
qr = [0 0]';
end

5 个评论
显示 3更早的评论隐藏 3更早的评论

Bill Greene 2023-8-30

Referencing the github page is fine.

carlos Hernando 2024-9-24

Hi @Bill Greene. Thank you for your code. I do have a quick question, is it possible to use vectorize calculation in your coupledMassExample.

Thank you!

请先登录，再进行评论。

Answer 2

Torsten 2023-8-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2013547-pdepe-with-cross-elements-in-the-differentiation-c#answer_1295251

移动：Torsten 2023-8-28

More basically: can pdepe accept matrices or just vector inputs in c, s, and f?

Why not reading the documentation of "pdepe" ?

pl, ql, pr and qr must be vectors.

s must be a vector.

c must be a diagonal matrix with elements either zero or positive.

f must be a vector.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Francesca Turco 2023-8-28

编辑：Francesca Turco 2023-8-29

在 MATLAB Online 中打开

Ok, this shoudl work, but it gives an error...

I added the d/dx(dUdx*U) terms to f, and subtracted the 2*dUdx(1)*dUdx(3) term in s, which should give the right UdUdx2 in total. But it gives me the error of wrong size vectors:

"Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 5: Error in burn_TnP_piD (line 8)

[sol] = pdepe(m,@(x,t,U,dUdx)pdefun(x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPd,fPrad,Snbi,Spump, ..."

The full code is attached, the coefficients are passed from burn_inputs and Lzinputs, etc.

All the f, s and c entries seem to have the right format...

function [x,t,T,ne,Prad,PrZ,Pbrem,Pline,Pin,Te,Lz,Lzb,Palpha]=burn_TnP_piD;
x=linspace(0,1,151); t=linspace(0,2,100); makeplots=0; 
burn_inputs; load Lzinputs; makeplots=0; Te=newTe; Lz=newLzKr; m=1; 
Vol=2*pi^2*1.75*(x*0.56).^2*1.85*(1-(1-8/3/pi)*0.59/1.75); 
[sol] = pdepe(m,@(x,t,U,dUdx)pdefun(x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPd,fPrad,Snbi,Spump, ...
                                    xx,cz0psm,neped,T0p,n0p,Pin0p,Tout,Vol,alfa,tshift), ...
                @(x)pdeic(x,Te,Lz,xx,cz0psm,T0p,n0p,Pin0p,alfa,tshift,Vol), ...
                @pdebc, ...
                x,t);
T=sol(:,:,1); Pin=sol(:,:,2); ne=sol(:,:,3);
clear nutty; nutty=sol(:,:,4);   % b_eff
for i=1:size(nutty,2), bE(:,i)=gradient(nutty(:,i),t); end
for i=1:length(t), Lzb(i,:)=bE(i,:)/1e34./cz0psm; end 
PrZ=bE.*ne.^2; Pbrem=brem*0.00534*2*ne.^2.*sqrt(T); Pline=liner*0.01*2*ne.^2; Prad=PrZ+Pbrem+Pline; 
clear nutty; nutty=sol(:,:,5);   % Palpha
for i=1:size(nutty,2), Palpha(:,i)=gradient(nutty(:,i),t); end
disp(['  T0=' num2str(T(end,1))]); 
disp(['  Palpha0=' num2str(Palpha(end,1))])
if makeplots == 1, plotburn_delay, end
end
function [c,f,s] = pdefun (x,t,U,dUdx,Te,Lz,DT,DN,fTpin,fTpr,fTout,L,brem,liner,fPp,fPd,fPrad,Snbi,Spump, ...
                           xx,cz0psm,neped,T0p,n0p,Pin0p,Tout,Vol,alfa,tshift)
intT0=pchip(xx,T0p,x); intPin0=pchip(xx,Pin0p,x); intn0=pchip(xx,n0p,x); 
intnep=pchip(xx,neped,x); intcz=pchip(xx,cz0psm,x); intV=pchip(xx,Vol,x);
intDT=pchip(xx,DT,x); intDN=pchip(xx,DN,x); intTout=pchip(xx,Tout,x);
intfPp=pchip(xx,fPp,x); 
Pr=L*1e34*interp1(Te,Lz,U(1),'pchip','extrap')*U(3)^2*intcz+ brem*0.00534*2*U(3)^2*sqrt(U(1))+ liner*0.01*2*U(3)^2;
Pa=U(1)^2*(0.78*U(3))^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(U(1)+tshift)-0.0089152*(U(1)+tshift).^2+0.000103742*(U(1)+tshift).^3);
berror=U(1)*U(3)-intT0*intn0; 
c=[1 1 1 1 1]';
f=[intDT*dUdx(1) -fPd*(U(3)*dUdx(1)+U(1)*dUdx(3)) intDN*dUdx(3) 0 0]'; 
s=[fTpin*U(2)-fTpr*Pr-intTout*U(1)+alfa*Pa, ...  % Te
   -intfPp*(berror)+fPrad*Pr -fPd*U(3)*(fTpin*U(2)-fTpr*Pr-intTout*U(1)+alfa*Pa) -fPd*U(1)*(Snbi*U(2)+intnep-Spump*U(3)) +2*fPd*dUdx(1)*dUdx(3), ...    % Pin   
   Snbi*U(2)+intnep-Spump*U(3), ...                            % ne
   1e34*interp1(Te,Lz,U(1),'pchip','extrap')*intcz, ...    % coefficient x Prad(Lz)
   alfa*Pa]';  % Palpha 
end
function u0 = pdeic(x,Te,Lz,xx,cz0psm,T0p,n0p,Pin0p,alfa,tshift,Vol)
intT0=pchip(xx,T0p,x); intPin0=pchip(xx,Pin0p,x); intn0=pchip(xx,n0p,x); 
intcz=pchip(xx,cz0psm,x); intLz0=interp1(Te,Lz,intT0,'pchip','extrap'); intV=pchip(xx,Vol,x);
u0=[intT0, intPin0, intn0, 1e34*intLz0*intcz, alfa*intT0^2*(0.78*intn0)^2*4.79e+16*1.64e-19*(-0.455383+0.218998*(intT0+tshift)-0.0089152*(intT0+tshift).^2+0.000103742*(intT0+tshift).^3)]';
end
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) 
pl = [0 0 0 0 0]';                   % Te, Pin, ne, coeff, Palpha
ql = [1 1 1 1 1]'; 
%pr = [ur(1)-0.01 ur(2)-0.02 ur(3)-0.06 0 0]';  % Te, Pin, ne, coeff, Palpha
pr = [ur(1)-0.01 0 ur(3)-0.06 0 0]';
qr = [1 1 1 1 1]';
end

Francesca Turco 2023-8-29

编辑：Francesca Turco 2023-8-29

It works!! The error was because... I had 2 spaces between sum terms in the s(2) expression... which told s to have extra elements *eyeroll goes here* Even the BC seems to be ok: I am now using pl(2)=0, ql(2)=0, pr(2)=ur(2)-0.02 and qr(2)=0, which seems to be keeping the Pin (U(2)) down to the low value I need at the outer edge.

I had to put a differentiable function in fPd instead of a scalar, because I have profiles with a "pedestal", ie a relatively sharp drop near the outer edge, and the derivative of that was giving edge spikes in Pin. The effect of these new terms should be constrained near the centre of the spatial domain anyway, so getting a coefficient that goes to zero smothly was necessary anyway.

I tested this script with fPd=0 against the one that doesn't have all the new terms, just as a sanity check, and it gives exactly the same results! fPd (the derivative gain, essentially) needs to be kept in check because it can prevent convergence if it gets too big, but this also is a typical behaviour for this physics topic, I believe.

I can't thank you enough. What is the best way to reference your help and work? If I can make this work generally, it is supposed to go into a Phys. Rev Letters and/or a Nucl. Fusion paper late rin the year, I would like to acknowledge your participation there.

Torsten 2023-8-30

编辑：Torsten 2023-8-30

I'd check the results with a second PDE solver (e.g. @Bill Greene 's code). Your model contains at least two ODEs without spatial derivatives (U(4) and U(5)), and the conditions pl(2)=0, ql(2)=0 practically mean that no boundary condition is set. Strictly speaking, "pdepe" is not designed for this kind of system of differential equations. Despite your enthusiasm, I'd classify the results as at least necessary to be evaluated.

Francesca Turco 2023-8-30

Agreed, I will try to implement Bill's solver tomorrow. The two ODE without spatial derivative here are just to extract two quantities that otherwise are difficult to check, ie they are conponents of the first 3 equations that I wanted to independently evaluate. I can actually remove both of them, I believe. The behaviour that I see when I scan the fPd vector of coefficients seems consistent with the origin profiles, so it is another indication that this method should be correct, besides recovering the results of the original system when these terms are set to zero (but exist in the system). I'll keep checking...

请先登录，再进行评论。

pdepe with cross elements in the differentiation "c"

2 个评论
显示无隐藏无

回答（2 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

pdepe with cross elements in the differentiation "c"

2 个评论 显示 无隐藏 无

回答（2 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

8 个评论 显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论