MATLAB Answers

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

6 views (last 30 days)
M S
M S on 25 Feb 2016
Commented: Ravi Kumar on 22 Oct 2018
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

Accepted Answer

Alan Weiss
Alan Weiss on 25 Feb 2016
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
  5 Comments
M S
M S on 25 Mar 2016
Yes, this suggestion is the best solution to our problem. Thank you Alan Weiss.Thank you MATLAB.

Sign in to comment.

More Answers (2)

Bruno Sainte-Rose
Bruno Sainte-Rose on 10 Jan 2017
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
  2 Comments
Bruno Sainte-Rose
Bruno Sainte-Rose on 11 Jan 2017
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

Sign in to comment.


juan ugalde
juan ugalde on 19 Oct 2018
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.
  1 Comment
Ravi Kumar
Ravi Kumar on 22 Oct 2018
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

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!