[BUG/ISSUE] Matlab Fun() issues

[jimmielin]

I'm trying to get a small mechanism running with Matlab/Octave but I'm running into a weird bug in the latest dev commit.

I have a new Matlab mechanism using small_strato, specified by m_small_strato.kpp:

#MODEL      small_strato
#LANGUAGE   matlab
#DOUBLE     ON
#INTEGRATOR rosenbrock
#DRIVER     general
#JACOBIAN   SPARSE_LU_ROW

{This is the small_strato example from Chapter 2 of the KPP manual}

Running kpp this generates the code fine but it crashes on Fun() - looking at m_small_strato_Fun.m, the code is obviously wrong:

function  [ Aout ] =  m_small_strato_Fun ( V , F , RCT , Vdot )

% Local variables
% A - Rate for each equation

A=zeros(1,length(RCT));
Vdot=zeros(1,length(V));

% Computation of equation rates
   A(1) = RCT(1)*F(2) ;
   A(2) = RCT(2)*V(2)*F(2) ;
   A(3) = RCT(3)*V(3) ;
   A(4) = RCT(4)*V(2)*V(3) ;
   A(5) = RCT(5)*V(3) ;
   A(6) = RCT(6)*V(1)*F(1) ;
   A(7) = RCT(7)*V(1)*V(3) ;
   A(8) = RCT(8)*V(3)*V(4) ;
   A(9) = RCT(9)*V(2)*V(5) ;
   A(10) = RCT(10)*V(5) ;

% Aggregate function
   Vdot(1) = A(5)-A(6)-A(7) ;
   Vdot(2) = 2*A(1)-A(2)+A(3)-A(4)+A(6)-A(9)+A(10) ;
   Vdot(3) = A(2)-A(3)-A(4)-A(5)-A(7)-A(8) ;
   Vdot(4) = -A(8)+A(9)+A(10) ;
   Vdot(5) = A(8)-A(9)-A(10) ;

return

% End of Fun function

For some reason, Aout is set as the output argument. This is wrong. Looking at PACT-1D-HALOGENS' current code at https://github.com/PACT1D/PACT-1D-HALOGENS/blob/main/mechanism/mech_Fun.m

the declaration of Fun() should only have (V, F, RCT) and Vdot is the output argument:

function  [ Vdot ] =  mech_Fun ( V , F , RCT )

% Local variables                                                  
% A - Rate for each equation                                       

A=zeros(1,length(RCT)); 
Vdot=zeros(1,length(V)); 
...

I suspect that the function declaration at this line is not correct for MATLAB, but I'm not sure how to fix this. For some reason

  if( z_useAggregate ) {
    FunctionBegin( F_VAR, V, F, RCT, Vdot, Aout );
  }
  else {
    FunctionBegin( FSPLIT_VAR, V, F, RCT, Vdot, P_VAR, D_VAR, Aout );
  }

Results in matlab generated code to have Aout as the output arg, not Vdot. I think Vdot is expected since Fun_Chem() uses it:

%  This line calls the Matlab ODE function routine  
  P = m_small_strato_Fun( Y, FIX, RCONST );

[jimmielin]

It looks like a fix to the function declarations is in order. I think Aout is only used by Fortran 90 since gen.c only generates it for Fortran:

  // Add code to return equation rates via optional argument Aout
  //   -- Bob Yantosca (29 Apr 2022)
  NewLines(1);
  if ( useLang == F90_LANG ) {
    fprintf(functionFile, "  !### Use Aout to return equation rates\n");
    fprintf(functionFile, "  IF ( PRESENT( Aout ) ) Aout = A\n");
  }

I've patched the function declaration for Fun() as follows:

  if(useLang == F90_LANG) {
    F_VAR = DefFnc( "Fun", 5,
                    "time derivatives of variables - Aggregate form");

    FSPLIT_VAR = DefFnc( "Fun_SPLIT", 7,
                         "time derivatives of variables - Split form");
  }
  else {
    F_VAR = DefFnc( "Fun", 4,
                    "time derivatives of variables - Aggregate form");

    FSPLIT_VAR = DefFnc( "Fun_SPLIT", 6,
                         "time derivatives of variables - Split form");
  }

  // We have added the capability to return equation rates and the
  // time derivative of variable species from Fun via optional arguments
  // Aout and VdotOut (when z_useAggregate=1)
  //   -- Bob Yantosca (03 May 2022)
  if( z_useAggregate ) {
    if(useLang == F90_LANG) { // F90_LANG adds Aout, but not needed for other languages
      FunctionBegin( F_VAR, V, F, RCT, Vdot, Aout );
    }
    else {
      FunctionBegin( F_VAR, V, F, RCT, Vdot );
    }
  }
  else {
    if(useLang == F90_LANG) { // F90_LANG adds Aout, but not needed for other languages
      FunctionBegin( FSPLIT_VAR, V, F, RCT, Vdot, P_VAR, D_VAR, Aout );
    }
    else {
      FunctionBegin( FSPLIT_VAR, V, F, RCT, Vdot, P_VAR, D_VAR );
    }
  }

  if ( (useLang==MATLAB_LANG)&&(!z_useAggregate) )
     printf("\nWarning: in the function definition for Fun_SPLIT.m move P_VAR to output vars\n");

This builds the code correctly for Matlab.

[jimmielin]

For some reason, Fun_Chem() gives a row vector which throws an error in the ODE solver:

>> run("m_small_strato_Main.m")
Initial Mass = 1.6240e+17

P =

   1.0e+06 *

   -2.0294    6.2742    2.0851    1.6565   -1.6565

Error using odearguments (line 93)
M_SMALL_STRATO_FUN_CHEM must return a column vector.

This seemed to be fine before, maybe it's a recent Matlab update? Vdot's declaration has been Vdot=zeros(1,length(V)); for quite some time, so it looks more like a MATLAB change.

[RolfSander]

Sorry, @jimmielin, I hardly know anything about Matlab. Maybe we should ask our Matlab colleagues directly here?

[jimmielin]

Hi Rolf, thanks for the suggestion, I'm still looking through the changes we made for the Matlab fixes for PACT-1D as they seem to be related. I'll follow up in that thread shortly once I figure out what changes are needed and confirm they don't break PACT-1D.

[yantosca]

We can now close this issue since PR #99 has been merged into version 3.1.1.