pdepe: [c f s]: extract f (the flux term)

Hi everyone,

I would need to extract the flux term, f, (not gradiaent) in the pdepe function. Does anyone know a way to extract f at certain times?

Thanks!

[Torsten]

During the solution process or afterwards ?

[Zana Taher]

Thanks for your reply! I actually meant during the process. I found out this way to solve it, if you have time looking at the code:

global exper exper=[-0.4,-112,-121,-128,-131,-131.1,-133.2,-133.1,-134.4,-134.7,-135.12,-135.7,-135.8,-135.1,-135.4,-135.8,-135.9,-134,-135.7,-135.1,-135.3,-135.4,-136,-139,-136]; % p=[0.218,0.52,0.0115,2.03,31.6]; % qr f a n ks p0=[0.218,0.52,0.0115,2.03,31.6]; lb=[0.01,0.4,0,1,15]; %lower bound ub=[Inf,Inf,10,10,100]; %upper bound options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter'); [p,resnorm,residual,exitflag,output] = lsqnonlin(@richards1,p0,lb,ub,options)

uu=richards1(p); global t plot(t,uu+exper') hold on plot(t,exper,'ko') function uu=richards1(p) % Solution of the Richards equation % using MATLAB pdepe % % $Ekkehard Holzbecher $Date: 2006/07/13 $ % % Soil data for Guelph Loam (Hornberger and Wiberg, 2005) % m-file based partially on a first version by Janek Greskowiak %-------------------------------------------------------------------------- L = 200; % length [L] s1 = 0.5; % infiltration velocity [L/T] s2 = 0; % bottom suction head [L] T = 4; % maximum time [T] % qr = 0.218; % residual water content % f = 0.52; % porosity % a = 0.0115; % van Genuchten parameter [1/L] % n = 2.03; % van Genuchten parameter % ks = 31.6; % saturated conductivity [L/T]

x = linspace(0,L,100); global t t = linspace(0,T,25);

options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7) u = pdepe(0,@unsatpde,@unsatic,@unsatbc,x,t,options,s1,s2,p); global exper uu=u(:,1,1)-exper'; % figure; % title('Richards Equation Numerical Solution, computed with 100 mesh points');

% subplot (1,3,1); % plot (x,u(1:length(t),:)); % xlabel('Depth [L]'); % ylabel('Pressure Head [L]');

% subplot (1,3,2); % plot (x,u(1:length(t),:)-(x'*ones(1,length(t)))'); % xlabel('Depth [L]'); % ylabel('Hydraulic Head [L]');

% for j=1:length(t) % for i=1:length(x) % [q(j,i),k(j,i),c(j,i)]=sedprop(u(j,i),p); % end % end % % subplot (1,3,3); % plot (x,q(1:length(t),:)*100) % xlabel('Depth [L]'); % ylabel('Water Content [%]');

% ------------------------------------------------------------------------- function [c,f,s] = unsatpde(x,t,u,DuDx,s1,s2,p) [q,k,c] = sedprop(u,p); f = k.*DuDx-k; s = 0; end % ------------------------------------------------------------------------- function u0 = unsatic(x,s1,s2,p) u0 = -200+x; if x < 10 u0 = -0.5; end end

% ------------------------------------------------------------------------- function [pl,ql,pr,qr] = unsatbc(xl,ul,xr,ur,t,s1,s2,p) pl = s1; ql = 1; pr = ur(1)-s2; qr = 0; end %------------------- soil hydraulic properties ---------------------------- function [q,k,c] = sedprop(u,p) qr = p(1); % residual water content f = p(2); % porosity a = p(3); % van Genuchten parameter [1/L] n = p(4); % van Genuchten parameter ks = p(5); % saturated conductivity [L/T] m = 1-1/n; if u >= 0 c=1e-20; k=ks; q=f; else q=qr+(f-qr)*(1+(-a*u)^n)^-m; c=((f-qr)*n*m*a*(-a*u)^(n-1))/((1+(-a*u)^n)^(m+1))+1.e-20; k=ks*((q-qr)/(f-qr))^0.5*(1-(1-((q-qr)/(f-qr))^(1/m))^m)^2; end end end

[Torsten]

But the flux is "f = k.*DuDx-k" and you don't get it by a call to "sedprop", do you ?

[Zana Taher]

No, because the flux contains the partial derivative that is not included in the sedprop function. So I have to calculate k with the following code and the get f.

% for j=1:length(t) % for i=1:length(x) % [q(j,i),k(j,i),c(j,i)]=sedprop(u(j,i),p); % end

Once I had k, I can calculate f using pdeval. What do you think?

[Torsten]

Yes, that's correct. I just missed the call to "pdeval" in your code from above to get DuDx(j,i)