Why do I have error during usage of state.time for non-constant boundary condition ?

Dear all,

I am trying to solve a parabolic PDE with imported 3D geometry. I create the model and everything correctly but I decided to implement on one of the faces a boundary condition that depends on the time. Basically, what it needs to do is for each integer time-step say state.time=k we need to access azimuth(state.time). For this purpose I synchronized the time steps as follows:

nframes = 5; tlist = linspace(1,nframes,nframes);

Previously in the code I have specified the function

myufun14= @(region,state)sin(azimuth(state.time))

where I recall that azimuth is a vector. It gave and error :

"Subscript indices must either be real positive integers or logicals"

I tried just to test whether state.time might not return something different from the time steps included in tlist=linspace(1,nframes,nframes). In doing so I just attempted to apply the ceiling function:

myufun14= @(region,state)sin(azimuth(ceil(state.time+1)))

I got the same error:"Subscript indices must either be real positive integers or logicals"

Please could you help me!

Best

Answer 1

[Alan Weiss]

You need to write your boundary conditions in the appropriate syntax. I suppose that you are setting a Dirichlet condition on your specified face.

As the documentation states, your boundary function must return an N1-by-M matrix of points, where N1 is the number of equations (probably 1 for your problem), and M is the number of evaluation points, meaning the length of region.x. So if sin(azimuth(state.time)) is a scalar, then your function should be

 myufun14 = @(region,state)sin(azimuth(state.time))*ones(1,length(region.x))

It is also possible that the solver would need times other than those you specify. So you might need to change azimuth(state.time) to azimuth(floor(state.time)).

Alan Weiss

MATLAB mathematical toolbox documentation

[M S]

We have tried your suggestion but precisely the same error "Subscript indices must either be real positive integers or logicals" persists. The most important parts of the code follow:

myufun14= @(region,state) sin(azimuth(floor(state.time)))*ones(1,length(region.x));

applyBoundaryCondition(model,'Face',1,'u',myufun14,'Vectorized','off'); ( we have attempted with 'on' as well)

nframes = 5;

tlist = linspace(1,nframes,nframes);

model.generateMesh('hmax',0.25/2);

pdemesh(model);

u0=0;

u1 = parabolic(u0,tlist,model,c,a,f,d);

We have specified c,a,f,d as simple scalars. Thanks in advance for your help.

[Alan Weiss]

Sorry, I do not know what is going on in your case. I suggest that you either give us more details so that we can reproduce the problem, or contact technical support.

Alan Weiss

MATLAB mathematical toolbox documentation

[M S]

We attach two files: workspace and pdeModel;

The file workspace contains the workspace with preloaded matrix azimuth and the executed code of the file pdeModel, which has imported a stl geometry, has initiated a PDE model with boundary conditions. At the line: u1 = parabolic(u0,tlist,model,c,a,f,d);

there appears a compilation error which obviously comes from the application of the boundary conditions via the line:

myufun15= @(region,state) sin(azimuth(floor(state.time)))*ones(1,length(region.x));

Thanks in advance!

[Alan Weiss]

I am happy to tell you that I finally figured out, with the help of a team member, what is going on. Even better, I have a solution for you.

To understand the problem, you need to know that the toolbox solvers such as parabolic test to see whether a problem is nonlinear by passing NaN values. They then look and see if there are any NaN in the output, and, if so, they treat the problem as nonlinear. What happened to your function is this. When the time passed is NaN, your function azimuth(time) errors. It does not return a NaN.

The workaround is clear. Rewrite your function as follows:

function u = mynewufun(region,state,azimuth)

if isnan(state.time)
    u = NaN(size(region.x));
else
    u = sin(azimuth(floor(state.time)))*ones(1,length(region.x));
end

Pass @(region,state)mynewufun(region,state,azimuth) as your boundary condition and all will be well.

It is entirely possible that this your troubles are my fault; I did not document everywhere that custom functions need to have this behavior. I did note it for coefficients (at the bottom of the page), and in one description of boundary conditions, but not in the more modern boundary condition section. I apologize.

I hope that you will now be able to continue your work.

Alan Weiss

MATLAB mathematical toolbox documentation

[M S]

Yes, this suggestion is the best solution to our problem. Thank you Alan Weiss.Thank you MATLAB.

Answer 2

[Bruno Sainte-Rose]

Dear all,

Thank you for your answer of this topic, they helped me a lot, but I still have a stupid error that I am not able to debug. My situation is quite similar to the one of M S.

I have the environmental data of the sun irradiance stored in the vector "Sun_HF_data" along one day in steps of quarter of an hour (1x96). I need to change the Neumann bc at each time step (tlist = 0:1:96)taking the value from Sun_HF_data.

I built my bc function in the following way as suggested above:

function bc = mygfun(region,state,Sun_HF_time)

if isnan(state.time) bc = NaN(size(region.x)); else bc = Sun_HF_time(floor(state.time))*ones(1,length(region.x)); end end

but I am getting the following error: Reference to non-existent field 'time'.

Help is much appreciated, need to fix it asap

[Alan Weiss]

Please tell us your MATLAB version, and show us your exact function calls.

Alan Weiss

MATLAB mathematical toolbox documentation

[Bruno Sainte-Rose]

My Matlab version is 2016b

Here you can find the mygfun.m file script

function bc = mygfun(region,state,Sun_HF_time)
 if isnan(state.time) 
     bc = NaN(size(region.x));       
 else
     bc = Sun_HF_time(floor(state.time))*ones(1,length(region.x));
 end 
end

And here below you can find the relevant code lines of the pde simulation:

%pde toolbox initialization
 numberOfPDE=1;                                                    
 model = createpde(numberOfPDE);

% Geometry definition
......

% Mesh Generation
......

% Steady State Solution
......%boundary condition and coefficient definition
results = solvepde(model);

% Transient Solution
Sun_HF_time = Sun_HF_data(:,str2double(INPUT{1,1})); %read sun irradiance data from input file
tlist = 0:1:96;                                      %time history every 15 min

% Boundary condition
.....
External_pipe_1200 = applyBoundaryCondition
(model,'neumann','Edge',5:6,'q',C_air,'g',@mygfun);

% Solver Transient Solution
c = k_hpde;
d = Cp_hpde*rho_p100;
specifyCoefficients(model,'m',0,'d',d,'c',c,'a',0,'f',0);
setInitialConditions(model,results);
Results_time = solvepde(model,tlist);

Where: C_air is a scalar Sun_HF_time is a vector Nx1.

Thank you for your help. Much appreciated

Answer 3

[juan ugalde]

Hello,

I'am having a similar problem with my heatSourceFun. When I use a time dependent variable inside the function, such as heat=QQ(floor(state.time)) / or/ heat=QQ(state.time), there is no output. I tried using the proposed solution using if isnan, but it still doesn't work. The code works when I specify a discrete time (x=1,100,...)heat=QQ(x). Any clues?

Here is the code:

function heatSourceValue = heatSourceFun(location,state)

load('total.mat')
Volume_cell= (0.07)*pi*(0.0105)^2;

if isnan(state.time)

      heatSourceValue= NaN(size(location.x));

elseif    state.time<=3500
         heat=QQ(floor(state.time));% QQ loaded from total.mat

           heatSourceValue= (heat/Volume_cell)*ones(size(location.x));

elseif   state.time>3500
         heatSourceValue= state.time*zeros(size(location.x));

end

end

Thanks in advance.

[Ravi Kumar]

Hi Juan,

Do you see the correct values of heatSourceValue for some intermediate value of time, say state.time near 1000?

If you insert a breakpoint in the second branch of the if block, you should be able set set a condition such as state.time>1000, then code will run till 1000 and stop at the first instant state.time exceeds 1000. Once you stop in the debug mode, make sure the values heatSourceValue is correct for that time range.

Regards, Ravi