#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.68 2002/04/15 10:54:31 geraint
## Statically declare outputs and initialise to zero.
## This is necessary to prevent spurious values from being output when no assignments are made (i.e. when "y(i) := 0 for all u" (Reduce:see NERO)).
##
## Revision 1.67 2002/04/09 12:04:21 geraint
## Replaced ios:: with std::ios:: for g++-3.0 compatability.
##
## Revision 1.66 2002/03/26 11:58:58 geraint
## Added cputime monitoring.
##
## Revision 1.65 2001/11/15 06:24:11 geraint
## Updated (-i dassl) residual function to use new DAEFunc (octave-2.1.35).
## YZ residual dependency on Ui still requires some work.
##
## Revision 1.64 2001/08/08 02:15:00 geraint
## Rationalisation of solver code, beginning with algebraic solvers.
##
## Revision 1.63 2001/08/07 04:39:24 geraint
## Consolidated dassl and residual functions.
##
## Revision 1.62 2001/08/01 22:14:32 geraint
## Bug fix for dassl.
##
## Revision 1.61 2001/08/01 04:06:07 geraint
## Added -i dassl for -cc and -oct.
##
## Revision 1.60 2001/07/16 22:23:00 geraint
## Fixed misleading variable name in .cc rep.
##
## Revision 1.59 2001/07/13 04:54:04 geraint
## Branch merge: numerical-algebraic-solution back to main.
##
## Revision 1.58 2001/07/13 00:51:39 geraint
## Fixed generation of odes.sg from .m and .oct simulations.
## .cc, .m and .oct simulations now all write mtt_data (lower case).
##
## Revision 1.57.2.5 2001/07/13 04:02:31 geraint
## Implemented numerical algebraic solution for _ode2odes.oct.
##
## Revision 1.57.2.4 2001/07/02 00:34:56 geraint
## gcc-3.0 compatibility.
##
## Revision 1.57.2.3 2001/06/25 23:28:29 geraint
## Generic mtt_rate and mtt_output - allows method independent calls.
##
## Revision 1.57.2.2 2001/06/05 03:20:40 geraint
## added -ae option to select algebraic equation solution method.
##
## Revision 1.57.2.1 2001/05/04 04:07:24 geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.57 2001/04/01 03:38:54 geraint
## Reset row to zero after write to file, ready for subsequent runs.
## Eliminates SIGSEGV in Octave when _ode2odes called multiple times.
##
## Revision 1.56 2001/03/30 15:13:58 gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.55 2001/03/27 13:21:59 geraint
## Octave version compatibility for save_ascii_data(_for_plotting).
##
## Revision 1.54 2001/03/27 01:14:27 geraint
## Improved determination of Octave version.
##
## Revision 1.53 2001/03/21 03:24:59 geraint
## Calculate inputs before outputs (.cc).
##
## Revision 1.52 2001/03/19 02:28:52 geraint
## Branch merge: merging-ode2odes-exe back to MAIN.
##
## Revision 1.51.2.7 2001/03/17 09:51:07 geraint
## Implemented Runge-Kutta IV fixed-step method (-i rk4).
##
## Revision 1.51.2.6 2001/03/16 03:56:13 geraint
## Removed psignal/siginfo.h - problematic and unnecessary.
##
## Revision 1.51.2.5 2001/03/12 23:16:37 geraint
## Minor improvements to signal handling (.exe).
##
## Revision 1.51.2.4 2001/03/12 03:59:30 geraint
## SIGINT (C-c C-c) now causes simulation data to be dumped to MTT.core.
## SIGQUIT (C-c C-\) as for SIGINT, then raises default SIGQUIT.
## SIGFPE as for SIGINT, then raises default SIGABRT.
##
## Revision 1.51.2.3 2001/03/07 04:06:55 geraint
## Irix: catch SIGFPE and write data before aborting (.exe).
## GNU/Linux: nada.
##
## Revision 1.51.2.2 2001/03/02 00:45:21 geraint
## Separated Euler and Implicit methods in .cc code and dependencies.
##
## Revision 1.51.2.1 2001/03/01 05:05:53 geraint
## Minor revisions.
##
## Revision 1.51 2001/02/19 06:33:19 geraint
## Removed operation form loop.
##
## Revision 1.50 2001/02/18 09:18:49 geraint
## Removed temporary Matrices from mtt_implicit.cc
##
## Revision 1.49 2001/02/14 06:06:34 geraint
## Removed octave_value_list wrappers from standalone.exe - speed improvements
##
## Revision 1.48 2001/02/11 07:08:59 geraint
## Static declarations of octave_value_lists: small .exe speed improvement
##
## Revision 1.47 2001/02/11 05:25:52 geraint
## Reduced number of matrix operations during .oct simulation data write
##
## Revision 1.46 2001/02/05 08:32:31 geraint
## typo
##
## Revision 1.45 2001/02/05 04:32:35 geraint
## Octave version 2.1.x compatability and #ifdef statements for standalone rep
##
## Revision 1.46 2001/01/08 06:21:59 geraint
## #ifdef STANDALONE stuff
##
## Revision 1.45 2001/01/07 01:25:49 geraint
## Compatibility with Octave 2.1.33
##
## Revision 1.44 2000/12/05 12:11:45 peterg
## Changed function name to name()
##
## Revision 1.43 2000/12/04 10:59:40 peterg
## *** empty log message ***
##
## Revision 1.42 2000/11/10 14:19:50 peterg
## Corrected the csex and cseo functions
##
## Revision 1.41 2000/11/09 17:06:39 peterg
## Now does euler for cc
##
## Revision 1.40 2000/10/17 09:55:00 peterg
## Replaced switchopen by logic
##
## Revision 1.39 2000/10/14 08:04:40 peterg
## Changed arguments to _inout for consistency
##
## Revision 1.38 2000/10/11 09:08:08 peterg
## cse --> csex
##
## Revision 1.37 2000/08/01 12:25:06 peterg
## Now includes euler
##
## Revision 1.36 2000/05/19 17:48:16 peterg
## Argument to state
##
## Revision 1.35 2000/05/18 18:59:40 peterg
## Removed the First time stuff
##
## Revision 1.34 2000/05/16 18:56:14 peterg
## *** empty log message ***
##
## Revision 1.33 2000/05/11 19:33:18 peterg
## Uniform version for _sim.m
##
## Revision 1.32 2000/05/11 08:30:00 peterg
##
## Revision 1.31 2000/05/10 18:33:25 peterg
## Use smxa and smxax in place of smx
##
## Revision 1.30 2000/04/18 11:24:19 peterg
## Removed _numpar.
##
## Revision 1.29 2000/04/07 19:10:57 peterg
## *** empty log message ***
##
## Revision 1.28 1999/12/08 05:56:52 peterg
## Reordered the writing of the input and output.
## Note that last value now discarded.
##
## Revision 1.27 1999/11/15 22:47:53 peterg
## Generates method-specific code.
##
## Revision 1.26 1999/10/20 01:31:43 peterg
## *** empty log message ***
##
## Revision 1.25 1999/08/29 06:55:26 peterg
## Removed [MTTu] = zero_input($Nu); # Zero the input
## to avoide the p2c bug ????
##
## Revision 1.24 1999/08/27 06:02:16 peterg
## removed zero_input to avoid p2c bug
##
## Revision 1.23 1999/08/02 13:39:19 peterg
## Replaced zero_vector by zero_input
##
## Revision 1.22 1999/04/20 06:16:07 peterg
## Removed initialisation of AA and AAx
## Remove _switch calls -- uses _switchopen exclusively
##
## Revision 1.21 1999/04/02 06:29:25 peterg
## New implicit method - solves numerical prob with ISW
##
## Revision 1.20 1999/04/02 02:13:58 peterg
## Back to RCS
##
## Revision 1.19 1999/03/30 21:39:25 peterg
## In implicit approach, set derivatives to zero (when switch is off)
## before update. This seems to stop numerical leakage though non-return
## switches.
##
## Revision 1.18 1999/03/15 01:17:07 peterg
## Removed some spurious debugging code
##
## Revision 1.17 1999/03/15 01:09:15 peterg
## Fixed bugs when Nx=0 (no state)
##
## Revision 1.16 1999/03/06 02:28:38 peterg
## Rearranged evaluation to: state - input - output - write
##
## Revision 1.15 1999/03/06 02:19:43 peterg
## Changed args to _input
##
## Revision 1.14 1998/10/01 16:02:01 peterg
## Integration with switches handled separately fro Euler and Implicit.
##
## Revision 1.13 1998/09/30 17:41:24 peterg
## Implicit method now allows for switches via _switchA
##
## Revision 1.12 1998/08/27 08:55:18 peterg
## Mods to integration methods
##
## Revision 1.11 1998/08/25 12:28:31 peterg
## Move initila switch to after initial input
##
## Revision 1.10 1998/08/25 12:22:45 peterg
## Put _switch after update and also at initilisation
##
## Revision 1.9 1998/08/15 13:46:59 peterg
## New versions of integration routines
##
## Revision 1.8 1998/08/11 13:28:03 peterg
## Lowercase mttLAST etc
##
## Revision 1.7 1998/07/30 11:29:54 peterg
## Added implicit integration stuff
##
## Revision 1.6 1998/07/30 10:44:37 peterg
## INcluded othe integration methods.
##
## Revision 1.5 1998/07/26 11:02:20 peterg
## Put mtt or MTT in front of variable names to avoid clashes with
## globals
##
## Revision 1.4 1998/07/25 20:14:00 peterg
## update code added for flexibility and octave efficiency
##
###############################################################
# Bourne shell script: make_ode2odes
# Copyright (c) P.J.Gawthrop July 1998.
# Tell user
sys=$1
lang=$2
filename=${sys}_ode2odes.${lang}
if [ -n "$3" ]; then
method=$3
else
method=implicit
fi
if [ -n "$4" ]; then
algebraic_solver=$4
else
algebraic_solver="Reduce_Solver"
fi
echo Creating $filename with $method integration method
# Find system constants
Nx=`mtt_getsize $sys x` # States
Nu=`mtt_getsize $sys u` # Inputs
Ny=`mtt_getsize $sys y` # Outputs
case "$method" in
"implicit")
ode=csex
odeo=cseo
algorithm="mtt_implicit(x,dx,AA,AAx,ddt,$Nx,open_switches)"
;;
"dassl")
ode=ode
odeo=odeo
algorithm="mtt_dassl(x,u,t,par,dx,ddt,MTTNX,MTTNYZ,open_switches)"
;;
"euler" | "rk4" | *)
ode=ode
odeo=odeo
algorithm="mtt_euler(x,dx,ddt,$Nx,open_switches)"
;;
esac
make_m()
{
#lang_header $1 ode2odes m 'x,par,simpar' '[Y,X,t]' > $filename
mtt_header ${sys} ode2odes m > $filename
cat <<EOF >> $filename
global MTT_data;
if nargin<3
simpar = ${sys}_simpar();
[simpar.dt] = mtt_simpar_update;
endif
if nargin<2
par = ${sys}_numpar();
[par] = mtt_numpar_update(par);
endif
if nargin<1
[x0] = ${sys}_state(par);
[x0] = mtt_state_update(x);
endif
## Initialise
t = 0.0;
ddt = simpar.dt/simpar.stepfactor;
ilast = round(simpar.last/ddt)+1; # Total number of steps
x = x0;
## Following removed due to p2c bug
## [u] = zero_input($Nu); # Zero the input
for MTTi=1:$Ny
y(MTTi) = 0;
endfor;
mttj = 0;
for it = 1:ilast #Integration loop
[u] = ${sys}_input(x,y,t,par); # Input
[y] = ${sys}_$odeo(x,u,t,par); # Output
if mttj==0
mtt_write(t,x,y,$Nx,$Ny); # Write it out
endif
EOF
if [ "$method" = "rk4" ]; then
cat << EOF >> $filename
[k1] = ddt * ${sys}_${ode}(x,u,t,par);
[k2] = ddt * ${sys}_${ode}(x+k1/2,u,t+ddt/2,par);
[k3] = ddt * ${sys}_${ode}(x+k2/2,u,t+ddt/2,par);
[k4] = ddt * ${sys}_${ode}(x+k3,u,t+ddt,par);
[dx] = [k1 + 2.0 * [k2 + k3] + k4] / (6.0 * ddt);
EOF
else
cat << EOF >> $filename
[dx] = ${sys}_$ode(x,u,t,par); # State derivative
EOF
fi
if [ "$method" = "implicit" ]; then
cat<< EOF >> $filename
[AA] = ${sys}_smxa(x,u,ddt,par); # (I-Adt) and (I-Adt)x
[AAx] = ${sys}_smxax(x,u,ddt,par); # (I-Adt) and (I-Adt)x
EOF
fi
cat <<EOF >> $filename
[open_switches] = ${sys}_logic(x,u,t,par); # Switch logic
[x] = $algorithm; # Integration update
t = t + ddt; # Time update
mttj = mttj+1; # Increment counter
if mttj==simpar.stepfactor
mttj = 0; # Reset counter
endif
endfor; # Integration loop
## Create the output data
mtt_data = MTT_data;
endfunction
EOF
} # make_m
make_cc()
{
# get octave version
case `$MATRIX --version | awk -F\. '{print $2}'` in
0) # stable
vector_value=vector_value
feval_header=toplev.h
save_ascii_data_function=save_ascii_data
;;
1) # development
vector_value=column_vector_value
feval_header=parse.h
save_ascii_data_function=save_ascii_data_for_plotting
;;
*)
vector_value=column_vector_value
feval_header=parse.h
save_ascii_data_function=save_ascii_data_for_plotting
;;
esac
cat <<EOF > $filename
#include <octave/oct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/ov-struct.h>
#include <octave/variables.h>
#ifndef STANDALONE
#include <octave/${feval_header}>
#endif
#include "${sys}_def.h"
#include "${sys}_sympar.h"
#include "mtt_${algebraic_solver}.hh"
#include <fstream>
#ifdef STANDALONE
#include <csignal>
extern ColumnVector F${sys}_ae (
ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par);
extern ColumnVector F${sys}_input (
ColumnVector &x,
ColumnVector &y,
const double &t,
ColumnVector &par);
extern ColumnVector F${sys}_logic (
ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par);
extern ColumnVector F${sys}_numpar (
void);
extern Octave_map F${sys}_simpar (
void);
extern ColumnVector F${sys}_state (
ColumnVector &par);
extern ColumnVector F${sys}_${ode} (
ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par);
extern ColumnVector F${sys}_${odeo} (
ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par);
EOF
case "$method" in
"implicit")
cat <<EOF >> $filename
extern ColumnVector Fmtt_implicit (
ColumnVector &x,
ColumnVector &dx,
Matrix &AA,
ColumnVector &AAx,
const double &ddt,
const int &nx,
const ColumnVector &open_switches);
extern Matrix F${sys}_smxa (
ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par);
extern ColumnVector F${sys}_smxax (
ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par);
EOF
;;
"dassl")
cat <<EOF >> $filename
extern ColumnVector Fmtt_dassl (
ColumnVector &x,
const ColumnVector &u,
const double &t,
const ColumnVector &par,
const ColumnVector &dx,
const double &ddt,
const int nx,
const int nyz,
const ColumnVector &open_switches);
EOF
;;
"euler" | "rk4" | *)
cat <<EOF >> $filename
extern ColumnVector Fmtt_euler (
ColumnVector &x,
const ColumnVector &dx,
const double &ddt,
const int &nx,
const ColumnVector &open_switches);
EOF
;;
esac
cat <<EOF >> $filename
void set_signal_handlers (void);
#endif // STANDALONE
ColumnVector
mtt_ae (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_ae(x,u,t,par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_ae", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_input (ColumnVector &x,
ColumnVector &y,
const double &t,
ColumnVector &par)
{
static MTT::${algebraic_solver} ae(MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);
static ColumnVector u (MTTNU);
#ifdef STANDALONE
u = F${sys}_input (x, y, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (y);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_input", args, 1);
u = f(0).${vector_value} ();
#endif
if (MTTNYZ == 0)
{
return u;
}
else
{
return ae.solve(x,u,t,par);
}
}
inline ColumnVector
mtt_logic (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_logic (x, u, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_logic", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_numpar (void)
{
#ifdef STANDALONE
return F${sys}_numpar ();
#else
static octave_value_list args, f;
f = feval ("${sys}_numpar", args, 1);
return f(0).${vector_value} ();
#endif
}
inline Octave_map
mtt_simpar (void)
{
#ifdef STANDALONE
return F${sys}_simpar ();
#else
static octave_value_list args;
static Octave_map f;
f["first"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["first"];
f["dt"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["dt"];
f["last"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["last"];
f["stepfactor"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["stepfactor"];
f["wmin"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["wmin"];
f["wmax"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["wmax"];
f["wsteps"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["wsteps"];
f["input"] = feval ("${sys}_simpar", args, 1)(0).map_value ()["input"];
return (f);
#endif
}
inline ColumnVector
mtt_state (ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_state (par);
#else
static octave_value_list args, f;
args (0) = octave_value (par);
f = feval ("${sys}_state", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_rate (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_${ode} (x, u, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_${ode}", args, 1);
return f(0).${vector_value} ();
#endif
}
inline ColumnVector
mtt_output (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_${odeo} (x, u, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_${odeo}", args, 1);
return f(0).${vector_value} ();
#endif
}
EOF
case "$method" in
"implicit")
cat <<EOF >> $filename
inline ColumnVector
mtt_implicit (ColumnVector &x,
ColumnVector &dx,
Matrix &AA,
ColumnVector &AAx,
const double &ddt,
const int &nx,
const ColumnVector &open_switches)
{
#ifdef STANDALONE
return Fmtt_implicit (x, dx, AA, AAx, ddt, nx, open_switches);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (dx);
args (2) = octave_value (AA);
args (3) = octave_value (AAx);
args (4) = octave_value (ddt);
args (5) = octave_value ((double)nx);
args (6) = octave_value (open_switches);
f = feval ("mtt_implicit", args, 1);
return f(0).${vector_value} ();
#endif
}
inline Matrix
mtt_smxa (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_smxa (x, u, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_smxa", args, 1);
return f(0).matrix_value ();
#endif
}
inline ColumnVector
mtt_smxax (ColumnVector &x,
ColumnVector &u,
const double &t,
ColumnVector &par)
{
#ifdef STANDALONE
return F${sys}_smxax (x, u, t, par);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
f = feval ("${sys}_smxax", args, 1);
return f(0).${vector_value} ();
#endif
}
EOF
;;
"dassl")
cat <<EOF >> $filename
inline ColumnVector
mtt_dassl (ColumnVector &x,
const ColumnVector &u,
const double &t,
const ColumnVector &par,
const ColumnVector &dx,
const double &ddt,
const int &nx,
const int &nyz,
const ColumnVector &open_switches)
{
#ifdef STANDALONE
return Fmtt_dassl (x, u, t, par, dx, ddt, nx, nyz, open_switches);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (u);
args (2) = octave_value (t);
args (3) = octave_value (par);
args (4) = octave_value (dx);
args (5) = octave_value (ddt);
args (6) = octave_value (static_cast<double> (nx));
args (7) = octave_value (static_cast<double> (nyz));
args (8) = octave_value (open_switches);
f = feval ("mtt_dassl", args, 1);
return f(0).${vector_value} ();
#endif
}
#ifdef STANDALONE
ColumnVector
Fmtt_residual (const ColumnVector &X,
const ColumnVector &DX,
double t,
int &ires)
{
#else // !STANDALONE
DEFUN_DLD (mtt_residual, args, , "mtt_residual")
{
static ColumnVector X (MTTNX+MTTNYZ);
static ColumnVector DX (MTTNX+MTTNYZ);
static double t;
static int &ires;
X = args(0).${vector_value} ();
DX = args(1).${vector_value} ();
t = args(2).double_value ();
ires = static_cast<int>(args(3).double_value ());
#endif // STANDALONE
static ColumnVector residual (MTTNX+MTTNYZ);
static ColumnVector U (MTTNU+MTTNYZ);
static ColumnVector u (MTTNU);
static ColumnVector y (MTTNY,0.0);
static ColumnVector par (MTTNPAR);
static ColumnVector dx(MTTNX);
static ColumnVector yz(MTTNYZ);
static ColumnVector x (MTTNX);
static ColumnVector ui (MTTNYZ);
static octave_value_list new_args;
x = X.extract (0,MTTNX-1);
if (MTTNYZ > 0)
ui = X.extract (MTTNX,MTTNX+MTTNYZ-1);
#ifdef STANDALONE
par = F${sys}_numpar();
u = F${sys}_input(x,y,t,par);
#else
par = feval ("${sys}_numpar", new_args, 1)(0).${vector_value} ();
new_args(0) = octave_value (x);
new_args(1) = octave_value (u);
new_args(2) = octave_value (t);
new_args(3) = octave_value (par);
u = feval ("${sys}_input", new_args, 1)(0).${vector_value} ();
#endif
U.insert (u,0);
if (MTTNYZ > 0)
U.insert (ui,MTTNU);
#ifdef STANDALONE
dx = F${sys}_${ode} (x,U,t,par);
yz = F${sys}_ae (x,U,t,par);
#else
new_args(1) = octave_value (U);
dx = feval ("${sys}_${ode}", new_args, 1)(0).${vector_value} ();
yz = feval ("${sys}_ae", new_args, 1)(0).${vector_value} ();
#endif
for (register int i = 0; i < MTTNX; i++)
residual (i) = dx(i) - DX(i);
if (MTTNYZ > 0)
{
residual.insert (yz,MTTNX);
// XXX:
// DASSL needs residual to be dependent on Ui and Uidot
// mtt_dassl always sets the initial Ui to zero, so
// Ui - h*Uidot should equal zero BUT, we don't know h
static double t_old;
double step;
if (t != t_old)
{
step = t - t_old;
t = t_old;
}
else
step = t;
for (register int i = MTTNX; i < MTTNX+MTTNYZ; i++)
residual(i) += X(i) - DX(i)*step;
}
#ifdef STANDALONE
return residual;
#else // !STANDALONE
return octave_value (residual);
#endif // STANDALONE
}
EOF
;;
"euler" | "rk4" | *)
cat <<EOF >> $filename
inline ColumnVector
mtt_euler (ColumnVector &x,
const ColumnVector &dx,
const double &ddt,
const int &nx,
const ColumnVector &open_switches)
{
#ifdef STANDALONE
return Fmtt_euler (x, dx, ddt, nx, open_switches);
#else
static octave_value_list args, f;
args (0) = octave_value (x);
args (1) = octave_value (dx);
args (2) = octave_value (ddt);
args (3) = octave_value ((double)nx);
args (4) = octave_value (open_switches);
f = feval ("mtt_euler", args, 1);
return f(0).${vector_value} ();
#endif
}
EOF
;;
esac
cat <<EOF >> $filename
inline void
mtt_write (const double &t,
ColumnVector &x,
ColumnVector &y,
const int &nrows,
const bool dump_data = false,
std::ostream &file = std::cout)
{
static Matrix data;
static int row;
if (dump_data)
{
if (row > 0)
{
Matrix written_data = data.extract (0, 0, row-1, data.cols ()-1);
$save_ascii_data_function (file, written_data, "mtt_data");
}
return;
}
const int nx = x.length (), ny = y.length ();
register int col = 0;
if (0 == row)
data = Matrix (nrows, 1+ny+1+nx, 0.0);
data.elem (row, col) = t;
for (register int i = 0; i < ny; i++)
data.elem (row, ++col) = y.elem (i);
data.elem (row, ++col) = t;
for (register int i = 0; i < nx; i++)
data.elem (row, ++col) = x.elem (i);
row++;
static std::fstream fcputime ("MTT.cputime", std::ios::out | std::ios::trunc | std::ios::app);
static clock_t cputime0 = clock();
static clock_t cputime1 = cputime0;
clock_t cputime = clock();
fcputime << t << '\t'
<< static_cast <double> (cputime - cputime0) / CLOCKS_PER_SEC << '\t'
<< static_cast <double> (cputime - cputime1) / CLOCKS_PER_SEC << std::endl;
cputime1 = cputime;
if (nrows == row)
{
#ifdef STANDALONE
$save_ascii_data_function (file, data, "mtt_data");
// std::cout << data << std::endl;
#else // ! STANDALONE
set_global_value ("MTT_data", data);
#endif
row = 0;
fcputime.close();
}
}
#ifdef STANDALONE
void dump_data (std::ostream &file)
{
ColumnVector null (0);
mtt_write (0.0, null, null, 0, true, file);
}
void handle_signal (int signum)
{
// handle some signals to ensure data is written.
std::cerr << "# Writing data to MTT.core (signal " << signum << ")" << std::endl;
std::ofstream corefile ("MTT.core");
dump_data (corefile);
switch (signum)
{
case SIGFPE:
// Intel chips do not raise SIGFPE for DIVZERO :-(
// raise (SIGABRT);
break;
case SIGINT:
break;
case SIGQUIT:
signal (SIGQUIT, SIG_DFL);
raise (SIGQUIT);
break;
default:
std::cerr << "# Warning: make_ode2odes needs updating!" << std::endl;
signal (signum, SIG_DFL);
raise (signum);
break;
}
corefile.close ();
set_signal_handlers ();
}
void set_signal_handlers (void)
{
signal (SIGFPE, handle_signal);
signal (SIGINT, handle_signal);
signal (SIGQUIT, handle_signal);
}
int main (void) {
set_signal_handlers ();
#else
DEFUN_DLD (${sys}_ode2odes, args, ,
"Octave ode2odes representation of system with $method integration method\nUsage: mtt_data = ${sys}_ode2odes (x0, par, simpar)\n")
{
static octave_value_list retval;
#endif // STANDALONE
static ColumnVector x0;
static ColumnVector par;
static Octave_map simpar;
static double
first = 0.0,
dt = 0.0,
last = 0.0,
stepfactor = 0.0;
#ifndef STANDALONE
int nargin = args.length ();
switch (nargin)
{
case 3:
first = args (2).map_value ()["first"].double_value ();
dt = args (2).map_value ()["dt"].double_value ();
last = args (2).map_value ()["last"].double_value ();
stepfactor = args (2).map_value ()["stepfactor"].double_value ();
par = args (1).${vector_value} ();
x0 = args (0).${vector_value} ();
break;
case 2:
first = mtt_simpar ()["first"].double_value ();
dt = mtt_simpar ()["dt"].double_value ();
last = mtt_simpar ()["last"].double_value ();
stepfactor = mtt_simpar ()["stepfactor"].double_value ();
par = args (1).${vector_value} ();
x0 = args (0).${vector_value} ();
break;
case 1:
first = mtt_simpar ()["first"].double_value ();
dt = mtt_simpar ()["dt"].double_value ();
last = mtt_simpar ()["last"].double_value ();
stepfactor = mtt_simpar ()["stepfactor"].double_value ();
par = mtt_numpar ();
x0 = args (0).${vector_value} ();
break;
case 0:
#endif // ! STANDALONE
first = mtt_simpar ()["first"].double_value ();
dt = mtt_simpar ()["dt"].double_value ();
last = mtt_simpar ()["last"].double_value ();
stepfactor = mtt_simpar ()["stepfactor"].double_value ();
par = mtt_numpar ();
x0 = mtt_state (par);
#ifndef STANDALONE
break;
default:
usage("${sys}_ode2odes (x par simpar)", nargin);
error("aborting.");
}
#endif // STANDALONE
static ColumnVector dx (MTTNX, 0.0);
static ColumnVector x (MTTNX, 0.0);
static ColumnVector u (MTTNU, 0.0);
static ColumnVector y (MTTNY, 0.0);
static Matrix AA (MTTNX, MTTNX, 0.0);
static ColumnVector AAx (MTTNX, 0.0);
static ColumnVector open_switches (MTTNX, 0.0);
register double t = 0.0;
const double ddt = dt / stepfactor;
const int ilast = static_cast<int> (round ((last - first) / ddt)) + 1;
const int nrows = static_cast<int> (round ((last - first) / dt)) + 1;
for (register int i = 0; i < MTTNY; i++)
{
y (i) = 0.0;
}
for (register int i = 0; i < MTTNX; i++)
{
x (i) = x0 (i);
}
for (register int j = 0, i = 1; i <= ilast; i++)
{
u = mtt_input (x, y, t, par);
y = mtt_output (x, u, t, par);
if (0 == j)
{
mtt_write (t, x, y, nrows);
}
EOF
case "$method" in
"rk4")
cat << EOF >> $filename
{
static ColumnVector
k1 (MTTNX,0.0),
k2 (MTTNX,0.0),
k3 (MTTNX,0.0),
k4 (MTTNX,0.0);
const double
t1 = t + ddt/2.0,
t2 = t + ddt;
ColumnVector
x1 (x),
x2 (x),
x3 (x);
k1 = ddt * mtt_rate (x , u, t , par); x1 += k1 * 0.5;
k2 = ddt * mtt_rate (x1, u, t1, par); x2 += k2 * 0.5;
k3 = ddt * mtt_rate (x2, u, t1, par); x3 += k3;
k4 = ddt * mtt_rate (x3, u, t2, par);
dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
}
EOF
;;
"dassl")
;;
"implicit")
cat << EOF >> $filename
dx = mtt_rate (x, u, t, par);
AA = mtt_smxa (x, u, ddt, par);
AAx = mtt_smxax (x, u, ddt, par);
EOF
;;
"euler" | *)
cat << EOF >> $filename
dx = mtt_rate (x, u, t, par);
EOF
;;
esac
## Common stuff
cat <<EOF >> $filename
open_switches = mtt_logic (x, u, t, par);
x = $algorithm;
t += ddt;
j++;
j = (j == static_cast<int> (stepfactor)) ? 0 : j;
}
#ifdef STANDALONE
return 0;
#else
/*
retval (0) = octave_value (y);
retval (1) = octave_value (x);
retval (2) = octave_value (t);
*/
retval = octave_value (get_global_value ("MTT_data"));
return (retval);
#endif // STANDALONE
}
EOF
}
case ${lang} in
m)
make_m
;;
cc)
make_cc
;;
*)
echo Language ${lang} is not supported
esac