Index: mttroot/mtt/bin/mtt
==================================================================
--- mttroot/mtt/bin/mtt
+++ mttroot/mtt/bin/mtt
@@ -15,10 +15,13 @@
###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
+## Revision 1.321 2001/07/27 23:43:34 geraint
+## Added -cc to usage options (required for use with xmtt).
+##
## Revision 1.320 2001/07/27 23:38:38 geraint
## Changed comment below (# SUMMARY) - fixes xmtt.
##
## Revision 1.319 2001/07/27 23:29:10 geraint
## *** empty log message ***
@@ -1196,10 +1199,14 @@
rdae_is_dae=1 ;;
-i )
mtt_switches="$mtt_switches $1";
shift;
case $1 in
+ dassl)
+ integration_method=dassl;
+ mtt_switches="$mtt_switches dassl";
+ ;;
euler)
integration_method=euler;
mtt_switches="$mtt_switches euler";
;;
implicit)
@@ -1209,11 +1216,11 @@
rk4)
integration_method=rk4;
mtt_switches="$mtt_switches rk4";
;;
*)
- echo $1 is an unknown integration method - use euler, rk4 or implicit;
+ echo $1 is an unknown integration method - use dassl, euler, rk4 or implicit;
exit;;
esac;;
-ae )
mtt_switches="$mtt_switches $1";
case $2 in
@@ -1380,11 +1387,11 @@
echo ' -cc C++ code generation'
echo ' -cr Use cr before resolving equations'
echo ' -d <dir> use directory <dir>'
echo ' -dc Maximise derivative (not integral) causality'
echo ' -dc Maximise derivative (not integral) causality'
- echo ' -i <implicit|euler|rk4> Use implicit, euler or rk4 integration'
+ echo ' -i <implicit|euler|rk4|dassl> Use implicit, euler, rk4 or dassl integration'
echo ' -ae <reduce|hybrd|hooke> Solve algebraic equations with Reduce, hybrd (fsolve) or Hooke and Jeeves'
echo ' -o ode is same as dae'
echo ' -oct use oct files in place of m files where appropriate'
echo ' -opt optimise code generation'
echo ' -p print environment variables'
@@ -2335,10 +2342,20 @@
$1_ode2odes_implicit.o : $1_cseo.o $1_csex.o $1_smxa.o $1_smxax.o mtt_implicit.o
@echo "Creating $1_ode2odes_implicit.o"
ar -cr \$@ \$^
$1_ode2odes_implicit.oct: $1_cseo.oct $1_csex.oct $1_smxa.oct $1_smxax.oct mtt_implicit.oct
@echo > /dev/null
+
+$1_ode2odes_dassl.m : $1_ode.m $1_odeo.m
+ @echo > /dev/null
+$1_ode2odes_dassl.cc : $1_ode.cc $1_odeo.cc
+ @echo > /dev/null
+$1_ode2odes_dassl.o : $1_ode.o $1_odeo.o mtt_dassl.o
+ @echo "Creating \$@"
+ ar -cr \$@ \$^
+$1_ode2odes_dassl.oct : $1_ode.oct $1_odeo.oct mtt_dassl.oct
+ @echo > /dev/null
mtt_Solver.cc: mtt_Solver.hh
$1_ode2odes_${algebraic_solver}.cc: mtt_Solver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc
$1_ode2odes_${algebraic_solver}.o: mtt_Solver.o mtt_${algebraic_solver}.o
@echo "Creating $1_ode2odes_${algebraic_solver}.o"
Index: mttroot/mtt/bin/trans/make_ode2odes
==================================================================
--- mttroot/mtt/bin/trans/make_ode2odes
+++ mttroot/mtt/bin/trans/make_ode2odes
@@ -7,10 +7,13 @@
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
+## 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.
@@ -259,21 +262,29 @@
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` # Inputs
-
-if [ "$method" = "implicit" ]; then
- ode=csex
- odeo=cseo
- algorithm="mtt_implicit(x,dx,AA,AAx,ddt,$Nx,open_switches)"
-else
- ode=ode
- odeo=odeo
- algorithm="mtt_euler(x,dx,ddt,$Nx,open_switches)"
-fi
+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
@@ -434,22 +445,13 @@
ColumnVector &u,
const double &t,
ColumnVector &par);
EOF
-if [ "$method" != "implicit" ]; then
-cat <<EOF >> $filename
-extern ColumnVector Fmtt_euler (
- ColumnVector &x,
- const ColumnVector &dx,
- const double &ddt,
- const int &nx,
- const ColumnVector &open_switches);
-
-EOF
-else
-cat <<EOF >> $filename
+case "$method" in
+ "implicit")
+ cat <<EOF >> $filename
extern ColumnVector Fmtt_implicit (
ColumnVector &x,
ColumnVector &dx,
Matrix &AA,
ColumnVector &AAx,
@@ -468,11 +470,38 @@
ColumnVector &u,
const double &t,
ColumnVector &par);
EOF
-fi
+ ;;
+ "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
@@ -554,10 +583,84 @@
static octave_value_list args, f;
f = feval ("${sys}_numpar", args, 1);
return f(0).${vector_value} ();
#endif
}
+
+#ifdef STANDALONE
+ColumnVector
+Fmtt_residual (const ColumnVector &X,
+ const ColumnVector &DX,
+ double t)
+{
+#else // !STANDALONE
+DEFUN_DLD (mtt_residual, args, , "mtt_residual")
+{
+ static ColumnVector X (MTTNX+MTTNYZ);
+ static ColumnVector DX (MTTNX+MTTNYZ);
+ static double t;
+
+ X = args(0).${vector_value} ();
+ DX = args(1).${vector_value} ();
+ t = args(2).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,MTTNX);
+
+#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);
+
+#ifdef STANDALONE
+ return residual;
+#else // !STANDALONE
+ return octave_value (residual);
+#endif // STANDALONE
+}
+
inline Octave_map
mtt_simpar (void)
{
#ifdef STANDALONE
@@ -627,36 +730,13 @@
return f(0).${vector_value} ();
#endif
}
EOF
-if [ "$method" != "implicit" ];then
-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
-else
-cat <<EOF >> $filename
+case "$method" in
+ "implicit")
+ cat <<EOF >> $filename
inline ColumnVector
mtt_implicit (ColumnVector &x,
ColumnVector &dx,
Matrix &AA,
ColumnVector &AAx,
@@ -717,11 +797,70 @@
return f(0).${vector_value} ();
#endif
}
EOF
-fi
+ ;;
+ "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
+}
+
+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,
@@ -905,12 +1044,13 @@
if (0 == j)
{
mtt_write (t, x, y, nrows);
}
EOF
-if [ "$method" = "rk4" ]; then
-cat << EOF >> $filename
+case "$method" in
+ "rk4")
+ cat << EOF >> $filename
{
static ColumnVector
k1 (MTTNX,0.0),
k2 (MTTNX,0.0),
k3 (MTTNX,0.0),
@@ -930,23 +1070,26 @@
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
-else
-cat << EOF >> $filename
+ ;;
+ "dassl")
+ ;;
+ "implicit")
+ cat << EOF >> $filename
dx = mtt_rate (x, u, t, par);
-EOF
-fi
-
-if [ "$method" = "implicit" ]; then
-cat <<EOF >> $filename
-
AA = mtt_smxa (x, u, ddt, par);
AAx = mtt_smxax (x, u, ddt, par);
EOF
-fi
+ ;;
+ "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;
ADDED mttroot/mtt/lib/cc/mtt_dassl.cc
Index: mttroot/mtt/lib/cc/mtt_dassl.cc
==================================================================
--- /dev/null
+++ mttroot/mtt/lib/cc/mtt_dassl.cc
@@ -0,0 +1,86 @@
+
+#include <octave/oct.h>
+#include <octave/DASSL.h>
+
+
+#ifdef OCTAVE_DEV
+#include <octave/parse.h>
+#define VECTOR_VALUE column_vector_value
+#else // !OCTAVE_DEV
+#include <octave/toplev.h>
+#define VECTOR_VALUE vector_value
+#endif // OCTAVE_DEV
+
+
+#ifdef STANDALONE
+extern ColumnVector
+Fmtt_residual (const ColumnVector &X, const ColumnVector &DX, double t);
+#endif // STANDALONE
+
+
+ColumnVector
+mtt_residual (const ColumnVector &X, const ColumnVector &DX, double t)
+{
+#ifdef STANDALONE
+ return Fmtt_residual (X, DX, t);
+#else // !STANDALONE
+ static octave_value_list args, f;
+ args(0) = octave_value (X);
+ args(1) = octave_value (DX);
+ args(2) = octave_value (t);
+ f = feval ("mtt_residual", args, 1);
+ return f(0).VECTOR_VALUE ();
+#endif // STANDALONE
+}
+
+
+#ifdef STANDALONE
+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 &openx)
+{
+#else // !STANDALONE
+DEFUN_DLD (mtt_dassl, args, ,
+ "dassl integration method")
+{
+ ColumnVector x = args(0).VECTOR_VALUE();
+ const ColumnVector u = args(1).VECTOR_VALUE();
+ const double t = args(2).double_value();
+ const ColumnVector par = args(3).VECTOR_VALUE();
+ const ColumnVector dx = args(4).VECTOR_VALUE();
+ const double ddt = args(5).double_value();
+ const int Nx = static_cast<int> (args(6).double_value());
+ const int Nyz = static_cast<int> (args(7).double_value());
+ const ColumnVector openx = args(8).VECTOR_VALUE();
+#endif // STANDALONE
+
+ static DAEFunc fdae(mtt_residual);
+ static ColumnVector XX (Nx+Nyz);
+ XX.insert (x,0);
+
+ for (register int i = Nx; i < Nyz; i++)
+ XX(i) = 0.0;
+
+ double tout = t + ddt;
+
+ DASSL fdassl (XX, t, fdae);
+ x = fdassl.do_integrate (tout).extract (0,Nx-1);
+
+ for (register int i = 0; i < Nx; i++)
+ if (openx (i) > 0.5)
+ x (i) = 0.0;
+
+
+#ifdef STANDALONE
+ return x;
+#else
+ return octave_value (x);
+#endif
+}