PDEPE function - BCFUN error

Hello all,

Can someone help me with the BCFUN error?

Thank you!

Error using pdepe

Unexpected output of BCFUN. For this problem BCFUN must return four column vectors of length 1.

Error in run_model (line 53)

sol = pdepe(0, @pdefun, @icfun, @bcfun, z_mesh, t_mesh);

Error in driver (line 58)

[z_mesh, t_mesh, T] = run_model(alpha, z_sensor, time, temp, env, hc, k, eps, sigma, A, Hs);

function [z_mesh, t_mesh, T] = run_model(alpha, z_sensor, time, t_mesh, temp, env, hc, k, eps, sigma, A,Hs)
% Create the time mesh using the measurement times [days].
% t_mesh is shifted to start at 0.
t_mesh = convertTo(time, 'datenum') - convertTo(time(1), 'datenum');
% Create the space mesh using even spacing between the top and 
% bottom sensors, plus the locations of the intermediate sensors.
z_mesh = unique([linspace(z_sensor(1), z_sensor(end), 100), z_sensor]);
% Solve the heat equation.
sol = pdepe(0, @pdefun, @icfun, @bcfun, z_mesh, t_mesh);

% Extract the first solution component as T.
T = sol(:,:,1);
%------------------------------------------------
% pdefun
%------------------------------------------------
    function [c, f, s] = pdefun(z, t, T, dTdz)
        c = 1;
        f = alpha * dTdz;
        s = 0;
    end
%------------------------------------------------
% icfun
%
% Notes:
% o The inital temperature profile is given by
%   interpolating the initial field data.
%------------------------------------------------
    function T0 = icfun(z)
        T0 = interp1(z_sensor, temp(1, :), z, 'linear');
    end
%------------------------------------------------
% bcfun
%
% Notes:
% o The boundary conditions through time are 
%   given by the flux and temperature at the bottom sensors. 
%------------------------------------------------
    function [ptop, qtop, pbot, qbot] = bcfun(ztop, Ttop, zbot, Tbot,t)
        % hc and Hs are vectors
        ptop = hc.*(Ttop-env(:,1))+ eps*sigma*Ttop^4-A*Hs;
        qtop = k;
        pbot = Tbot - interp1(t_mesh, temp(:,6), t, 'linear');
        qbot = 0;
    end
end

Answer 1

[Torsten]

Use

function [ptop, qtop, pbot, qbot] = bcfun(ztop, Ttop, zbot, Tbot,t)
% hc and Hs are vectors
ptop = hc.*(Ttop-env(:,1))+ eps*sigma*Ttop^4-A*Hs;
qtop = k;
pbot = Tbot - interp1(t_mesh, temp(:,6), t, 'linear');
qbot = 0;
size(ptop)
end

and see what the size of ptop is. It should be 1x1, but I doubt it because I guess that env(:,1) is a vector.

[Sanley Guerrier]

Yes, hc, env(:,1), Hs, and temp(:,6) are vectors of same length.

[Sanley Guerrier]

I tried size(ptop), it didn't work.

[Torsten]

Yes, hc, env(:,1), Hs, and temp(:,6) are vectors of same length.

ptop = interp1(t_mesh,hc,t)*(Ttop-interp1(t_mesh,env(:,1),t))+ eps*sigma*Ttop^4-A*interp1(t_mesh,Hs,t);
qtop = k

if you want to set the condition

interp1(t_mesh,hc,t)*(Ttop-interp1(t_mesh,env(:,1),t))+ eps*sigma*Ttop^4-A*interp1(t_mesh,Hs,t) + k*alpha*dT/dz = 0

at the top.

[Sanley Guerrier]

Thank you,

This command is working now but the output is only one row. Something is still not working.

ptop = interp1(t_mesh,hc,t)*(Ttop-interp1(t_mesh,env(:,1),t))+ eps*sigma*Ttop^4-A*interp1(t_mesh,Hs,t);

qtop = k

[Torsten]

This command is working now but the output is only one row.

What do you mean ? ptop must be one single value, namely you boundary condition evaluated at time t.

[Sanley Guerrier]

I meant the solution T should be a matrix length (z_mesh X t_mesh) but I run the code it only returned a row vector T.

[Sanley Guerrier]

When I used this top boundary condition it works fine. But when I switched to other one it didn't work.

function [ptop, qtop, pbot, qbot] = bcfun(ztop, Ttop, zbot, Tbot, t)

ptop = Ttop - interp1(t_mesh, temp(:,1), t, 'linear');

qtop = 0;

pbot = Tbot - interp1(t_mesh, temp(:,6), t, 'linear');

qbot = 0;

end

[Torsten]

I meant the solution T should be a matrix length (z_mesh X t_mesh) but I run the code it only returned a row vector T.

Don't you get a message that the integration stopped because the stepsize became too small or something similar ? If integration didn't proceed till the second output time, the only thing that pdepe returns are the initial conditions for T.

[Sanley Guerrier]

Warning: Failure at t=1.958356e+00. Unable to meet integration tolerances without reducing the step size below

the smallest value allowed (3.552714e-15) at time t.

> In ode15s (line 725)

In pdepe (line 291)

In run_model (line 53)

In driver (line 58)

Warning: Time integration has failed. Solution is available at requested time points up to t=1.958333e+00.

> In pdepe (line 305)

In run_model (line 53)

In driver (line 58)

Elapsed time is 18.348141 seconds.

[Torsten]

Yes, that's what I meant - my guess is that t_mesh(2) > 1.958356, isn't it ? So you have an explanation why you only get one row for T.

Concerning the error I cannot help. It seems the reason is the boundary condition at the top you supply. Check whether it's correct. Note that it reads ptop + qtop*f = 0 where f in your case is equal to alpha*dT/dz. Thus you have ptop + k*alpha*dT/dz = 0 as boundary condition.

[Sanley Guerrier]

Thank you, I got it.

[Sanley Guerrier]

I forgot to convert hc to m/hr because every was in m/hr.

Do you know why when I used qtop = -k works but qtop = k doesn't?

[Torsten]

You must decide why the boundary condition is

ptop - k*alpha*dT/dz = 0

and not

ptop + k*alpha*dT/dz = 0

It's equivalent to an inversion of the direction of the heat flux at the top.

[Sanley Guerrier]

That makes sence.

Thanks a lot for helping out.