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 use directory ' echo ' -dc Maximise derivative (not integral) causality' echo ' -dc Maximise derivative (not integral) causality' - echo ' -i Use implicit, euler or rk4 integration' + echo ' -i Use implicit, euler, rk4 or dassl integration' echo ' -ae 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 <> $filename -extern ColumnVector Fmtt_euler ( - ColumnVector &x, - const ColumnVector &dx, - const double &ddt, - const int &nx, - const ColumnVector &open_switches); - -EOF -else -cat <> $filename +case "$method" in + "implicit") + cat <> $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 <> $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 <> $filename +extern ColumnVector Fmtt_euler ( + ColumnVector &x, + const ColumnVector &dx, + const double &ddt, + const int &nx, + const ColumnVector &open_switches); + +EOF + ;; +esac cat <> $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 <> $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 <> $filename +case "$method" in + "implicit") + cat <> $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 <> $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 (nx)); + args (7) = octave_value (static_cast (nyz)); + args (8) = octave_value (open_switches); + f = feval ("mtt_dassl", args, 1); + return f(0).${vector_value} (); +#endif +} + +EOF + ;; + "euler" | "rk4" | *) +cat <> $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 <> $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 <> $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 <> $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 +#include + + +#ifdef OCTAVE_DEV +#include +#define VECTOR_VALUE column_vector_value +#else // !OCTAVE_DEV +#include +#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 (args(6).double_value()); + const int Nyz = static_cast (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 +}