MATLAB PDE Toolbox state.time not producing sensible values

I am using the new PDE workflow in MATLAB R2016a (but have been experiencing the same errors in the legacy workflow on R2015b, hence the upgrade to 2016a). I am solving a parabolic PDE on a 3D geometry. The f coefficient in the PDE is required to vary as a function of time. In particular, it needs to pulse on an off between 0 and a positive set of values. The relevant code is below:

duty_ratio = 0.5;
duty_period = 10;
duty_cycle = @(time) (mod(time, duty_period) < duty_period * duty_ratio);

fcoeff = @(region, state) duty_cycle(state.time) * constant * (pressure_interpolant([region.x' region.y' (region.z+z_offset)']).^2)';

specifyCoefficients(model,'m',0,...
                          'd',@dcoeff_new_matlab,...
                          'c',@ccoeff_new_matlab,...
                          'a',0,...
                          'f',fcoeff);

setInitialConditions(model, u0);

fprintf('%s: solving the heat PDE for power=%f watts\n', datestr(now),  power);
results = solvepde(model, tlist);

The duty_cycle(state.time) term is causing the problems. When that term is replaced by a constant, e.g. 1, the fcoeff function behaves as expected and the rest of the function computes fine - no problems there. The issue is at the level of the state.time parameter. Instead of producing sensible values (in this particular case between 0 and 70), it is producing nonsense values such as NaN and 2.2204e-16 (eps). I checked this by writing a dummy fcoeff function that just prints state.time and returns a constant. As a result, the duty_cycle function not work as it is not being supplied the correct value of state.time.

I have been fiddling with this for 24 hours now beleiving it to be my fault but it looks like MATLAB's - can anyone please suggest a workaround?

Also for information: the model has no specified boundaries and a scalar initial condition u0. The geometry is properly defined as I have been running this model without the duty_cycle term for several weeks.

Thanks Jack

Answer 1

[Alan Weiss]

Sorry about that. As I recently answered (contritely), the reason you are having trouble is probably my fault.

When solvers have a potentially nonlinear problem, they pass NaN values to see if the coefficients or boundary terms are NaN also. If so, then they use (internally) a nonlinear solver. So these NaN values are not nonsense, they are the way the solver probes to figure out what is going on. See the linked post for details (toward the end).

Also, you should know that the parabolic and hyperbolic equation solver can take any time in the stated range of times, not just the ones where you request solutions.

Rest assured that the documentation will be better, and more explicit, fairly soon. Again, my apologies for making this information hard to find.

Alan Weiss

MATLAB mathematical toolbox documentation

[Jack F]

Thanks Alan. Unfortunately the state.time parameter never seems to have an effect, even when the solver is actually 'solving' rather than 'probing'. Irrespective of the values I use, the state.time is always smaller than the rhs of the inequality and so the fcoeff is always 'on'. Could this be due to the fact that when the solver is probing with different values, it never tests a value for state.time that causes fcoeff to change and thus assumes that fcoeff actually has no time dependency and so cuts out state.time altogether?

Jack

[Alan Weiss]

Perhaps you missed a detail in the linked posting. You need to have your fcoeff function return NaN of the correct size when the solver passes a NaN value for state.time. If you do that, then rest assured that the solver knows that you have a time-dependent problem.

If you have done that already, then I am not sure what is going on. You say that the time never gets large enough? Does the solver finish and give you values for your entire tlist? If so, how did it get to the end of tlist and not use high enough values? Beats me. Did you ask for values in tlist that definitely make fcoeff take the value you want?

Sorry, that's all I can think of. It seems you know how to use the debugger, that's what I would do at this point.

Alan Weiss

MATLAB mathematical toolbox documentation

[Miguel Angel Arnedo]

I got exactly the same problem. Thee values of state.time are just st =NaN, or the extrange value st = 2.2204e-16. In my code I have taken care of the NaN cases as well. My system is time dependent, with N = 3, Neumann conditions. Please help

 applyBoundaryCondition(model,'Neumann','u',@uBc_coeffunc_I, 'Edge',3);
 .....
 function uBc = uBc_coeffunc_I(region,state)

        nr = length(region.x);  
        uBc = zeros(3,nr);
      
        st = state.time
        
        if  (~ isnan(st)) 
           % overwrites uBc(1,:)
            uBc(1,:) = st * 3;   % 3 for example, is dummy by now                    
        end                                   
 end

[Giulia Ulpiani]

I am having the same issue! Any help? Anyone found a solution?

[Andreas Apostolatos]

Dear all,

Even though this post is relatively old, I would like to provide an answer to the questions above just in case somebody else also stumbles upon the same question.

As Alan mentioned above, the function handle that you use in order to define a state-, time-dependent coefficient or boundary condition should probe the variable 'state.time' and check whether it is NaN or not. Our solvers use NaN internally as an indication of a potentially nonlinear problem. In case variable 'state.time' has value NaN, then the state-, time-dependent coefficient or boundary condition should be defined as a matrix with the expected number of rows and columns that have NaN values. Otherwise, in case variable 'state.time' has a value different than NaN, then you can define the values of the custom state-, time-dependent coefficient or boundary condition as you wish. As Alan pointed out, this information is mentioned on our documentation page below,

https://www.mathworks.com/help/pde/ug/pde.thermalmodel.thermalproperties.html#mw_dd1c63c7-c951-4525-996b-adaae3a6304a

where it is mentioned the following:

"For thermal conductivity, your function must return a matrix thermalVal with number of rows equal to 1, Ndim, Ndim*(Ndim+1)/2, or Ndim*Ndim, where Ndim is 2 for 2-D problems and 3 for 3-D problems. The number of columns must equal the number of evaluation points, for example, M = length(location.y). For details about dimensions of the matrix, see c Coefficient for specifyCoefficients.

If properties depend on the time or temperature, ensure that your function returns a matrix of NaN of the correct size when state.u or state.time are NaN. Solvers check whether a problem is time dependent by passing NaN state values and looking for returned NaN values."

In summary, your code within the function handle should contain such lines as,

% Function handle for the state-, time-dependent coefficient or boundary condition 
function outp = myfun(location, state)
    if ~isnan(state.time)
        outp = ones(numCols, length(location.x)); % expected number of columns 'numCols'
    else
        outp = nan(numCols, length(location.x));
    end
end

I hope that this information is useful.

Kind regards,

Andreas