[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Best ways to free memory within c++ DLD's
|
From: |
Jason Hoogland |
|
Subject: |
Best ways to free memory within c++ DLD's |
|
Date: |
Wed, 15 Feb 2006 17:01:16 +1000 |
|
User-agent: |
KMail/1.8 |
Q: What are the best ways to free up memory for Octave classes at the end of
DLD's?
I have an optimisation routine based on Michael Creel's samin.cc, that starts
with an octave script containing a call to samin which makes thousands of
calls to the objective function compress_obj.m, which in turn does the grunt
calculation of the physical process being modelled via two c++ DLD's.
compress_opt.m (script)
|
\|/
`
samin.oct -> compress_obj.m (fn)
|
\|/
`
compress_fn_oct.oct -> piston_fn_oct.oct
I made reference to Paul Thomas's excellent Dal Segno Al Coda and used
function pointers and reduced use of certain Octave classes to try and
streamline things. I've pasted the main .cc file below for reference
(contains no deliberate memory clean up). It all works well. However there
seems to be memory leakage which only became evident and debilitating when
the iterations became large. I'm relatively new to c++. My question is:
what are the best ways to free up memory for Octave classes? Any tips or
tricks?
Jason
--
==========================================================
Jason Hoogland hoogland at gmail dot com
Doctoral student ph(w) +61 7 3365 4457
Centre for Hypersonics ph(mob) +61 413 300 887
The University of Queensland UTC+10
Brisbane QLD 4072, Australia http://www.marsgravity.org
----------------------------------------------------------
compress_fn_oct.cc
#include <octave/oct.h>
#include <octave/ov-struct.h>
#include <octave/quit.h>
#include <octave/Cell.h>
#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdio>
#include <ext/stdio_filebuf.h>
const double pi = 3.1415926535897931;
const double eps = 2.2204460492503131e-16;
octave_function *lsode_ptr;
std::string strPistonFn( "piston_fn_oct" );
octave_value_list lsode_call( ColumnVector x0, ColumnVector t ) {
octave_value_list arg;
arg.append( octave_value( strPistonFn ) );
arg.append( octave_value( x0 ) );
arg.append( octave_value( t ) );
arg = lsode_ptr->do_multi_index_op ( 3, arg );
// int istate = arg( 1 ).int_value();
if ( error_state ) {
std::string msg = arg( 2 ).string_value();
octave_stdout << "error: lsode_call, msg= \"" << msg << "\"";
return octave_value( -1 );
} // if ( error_state )
return arg;
} // octave_value_list lsode_call
/* General Arrangement
|<--------------LCT1------------>|<-------------LCT2---------->|
Secondary x(t) ^
(inner) piston |-> v(t) | qdotdx Rubber buffer Primary
+------\-----+-+-----------------+---+ / diaphragm
| \ |/| |\\\|_/ |
| pr +---|/|---+ pCT +---+-------------------------+
| Tr | |/| | TCT |
| +---|/|---+ +---+-------------------------+
| \ |/| \ |\\\|
+----\-------+\+-----------\-----+---+
\ \ \
Reservoir \ Compression tube (CT) V2
(r) V1 Primary
(outer) piston
0 1 2 3 N
+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
Spatial discretisation N intervals
*/
DEFUN_DLD( compress_fn_oct , arg , ,
"Calculation for compression by two free pistons in a tube.\n\
argout = compress_fn_oct(argin)\n\
where argin and argout are parameter passing structures.") {
octave_value_list retval;// For this function
octave_value_list argout;// For piston ODE
// We don't care about the name of the structure being passed,
// only the top level fields.
Octave_map arg0 = arg(0).map_value();
std::string strLsode( "lsode" );
lsode_ptr = extract_function ( strLsode, "compress_fn_oct", "", "", "");
std::ofstream logfile;
FILE *pfile;
octave_stdout.setf( std::ios_base::scientific );
// *************** Output
int verbosity ( ( arg0.contents( arg0.seek( "verbosity" ) ) (0) ).\
int_value() );
int nlog = 100;
if ( verbosity > 0 ) {
nlog = ( arg0.contents( arg0.seek( "nlog" ) ) (0) ).\
int_value();
std::string pathLog ( ( arg0.contents( \
arg0.seek( "pathLog" ) ) (0) ).string_value() );
pfile = fopen( pathLog.c_str() , "w" );
/*
Ok, it took me a while to find the relevant info on the web
about this. I wanted to use a traditional C FILE for the
logfile for speed, but also use the octave_value.print
method to the same logfile for convenience when the verbosity
level = 3. Here we have created a C file pointer, "pfile",
and a C++ ofstream object "logfile". Windows C++
implementations tend to make available an attach method, e.g.
logfile.attach( pfile, std::ios::out ). For GCC, the library
has changed from version to version. To make logfile point to
the same file stream as pfile, I found this advice best:
http://list.wylug.org.uk/pipermail/wylug-help/2004-August/001572.html
http://gcc.gnu.org/onlinedocs/libstdc++/27_io/howto.html#11
Note #include <ext/stdio_filebuf.h> is needed. The following
two lines do the trick. pfile and logfile can later be closed
independently, without effecting the other.
*/
__gnu_cxx::stdio_filebuf<char> fileBuf( pfile , std::ios::out );
logfile.std::ios::rdbuf( &fileBuf );
logfile << "# Logfile created by compress_fn_oct.oct\n";
logfile << "Verbosity= " << verbosity << "\n";
logfile.flush();
if ( verbosity == 3 ) {
// Here the pointer to the original octave_stdout stream
// is obtained, and octave_stdout is redirected temporarily
// to the logfile.
std::streambuf *octout_orig = octave_stdout.std::ios::rdbuf();
octave_stdout.std::ios::rdbuf( &fileBuf );
// Iterate through the incoming structure and print to logfile
Octave_map::iterator p0;
for ( p0 = arg0.begin() ; p0 != arg0.end(); p0++ ) {
logfile << arg0.key( p0 ) << " =\n";
octave_value( arg0.contents( p0 ) ). \
print( octave_stdout , false );
}
// In case we have a shutdown later, ensure input values
// are flushed immediately to logfile.
logfile.flush();
octave_stdout.flush();
// Redirect octave_stdout back to original buffer
octave_stdout.std::ios::rdbuf( octout_orig );
} // if ( verbosity == 3 )
} // if ( verbosity > 0 )
// *************** Geometry
double dCT1 ( ( arg0.contents( arg0.seek( "dCT1" ) )
(0) ).scalar_value() );
double dCT2 ( ( arg0.contents( arg0.seek( "dCT2" ) )
(0) ).scalar_value() );
double A1 = pi * pow( dCT1 / 2.0, 2.0);// [m]^2 primary CT area
double A2 = pi * pow( dCT2 / 2.0, 2.0);// [m]^2 secondary CT area
double LCT1 ( ( arg0.contents( arg0.seek( "LCT1" ) )
(0) ).scalar_value() );
double LCT2 ( ( arg0.contents( arg0.seek( "LCT2" ) )
(0) ).scalar_value() );
double LCage ( ( arg0.contents( arg0.seek( "LCage" ) ) (0) ).\
scalar_value() );
double dCage ( ( arg0.contents( arg0.seek( "dCage" ) ) (0) ).\
scalar_value() );
// [m]^2 min area in cage/diaphragm
double ACage = pi * pow( dCage / 2.0, 2.0);
double Vr0 ( ( arg0.contents( arg0.seek( "Vr0" ) ) (0) ).scalar_value() );
double VCT0 = LCT1*A1 + LCT2*A2;// [m]^3 CT initial volume
Matrix D0 ( ( arg0.contents( arg0.seek( "D0" ) ) (0) ).matrix_value() );
// *************** Pistons
double m1 ( ( arg0.contents( arg0.seek( "m1" ) ) (0) ).scalar_value() );
double m2 ( ( arg0.contents( arg0.seek( "m2" ) ) (0) ).scalar_value() );
double grav = 9.8;// (9.8) [m]/[s]^2 acceleration due to gravity
double mufFudge ( ( arg0.contents( arg0.seek( "mufFudge" ) ) (0) ).\
scalar_value() );
double muf01 ( ( arg0.contents( arg0.seek( "muf01" ) ) (0) ).\
scalar_value() );
double muf02 ( ( arg0.contents( arg0.seek( "muf02" ) ) (0) ).\
scalar_value() );
double mufc ( ( arg0.contents( arg0.seek( "mufc" ) )
(0) ).scalar_value() );
// [m]^2 circumferential contact area of chevron seal
double ASeal = 0.02*pi*dCT2;
double dLeak ( ( arg0.contents( arg0.seek( "dLeak" ) ) (0) ).\
scalar_value() );
// [m]^2 circumferential projected area of chevron leak
double ALeak = dLeak*pi*dCT2;
// *************** Diaphragm
double DtOpening ( ( arg0.contents( arg0.seek( "DtOpening" ) ) (0) ).\
scalar_value() ) ;
double pBurst ( ( arg0.contents( arg0.seek( "pBurst" ) ) (0) ).\
scalar_value() ) ;
double DDiaphragmOpen ( ( arg0.contents( arg0.seek( "DDiaphragmOpen" ) ) \
(0) ).scalar_value() );
// [m]^2 avg area of diaphragm when open
double ADiaphragmOpen = pi*pow( DDiaphragmOpen/2.0, 2.0 );
// *************** Pressures
double pr0 ( ( arg0.contents( arg0.seek( "pr0" ) ) (0) ).scalar_value() );
double pCT0 ( ( arg0.contents( arg0.seek( "pCT0" ) )
(0) ).scalar_value() );
double pLossFudge ( ( arg0.contents( arg0.seek( "pLossFudge" ) ) (0) ).\
scalar_value() ) ;
// *************** Temperatures
double Tw ( ( arg0.contents( arg0.seek( "Tw" ) ) (0) ).scalar_value() );
double T0 = Tw;// (296) [K] common initial temperature for reservoir and
CT
// *************** Gas properties
double gr ( ( arg0.contents( arg0.seek( "gr" ) ) (0) ).scalar_value() );
double gCT ( ( arg0.contents( arg0.seek( "gCT" ) ) (0) ).scalar_value() );
double Rr ( ( arg0.contents( arg0.seek( "Rr" ) ) (0) ).scalar_value() );
double RCT ( ( arg0.contents( arg0.seek( "RCT" ) ) (0) ).scalar_value() );
double Pr ( ( arg0.contents( arg0.seek( "Pr" ) ) (0) ).scalar_value() );
double cpCT ( ( arg0.contents( arg0.seek( "cpCT" ) )
(0) ).scalar_value() );
double cvCT ( ( arg0.contents( arg0.seek( "cvCT" ) )
(0) ).scalar_value() );
// sonic throat constant for leak
double CeCT = pow( (gCT+1)/2.0, -(gCT+1)/(2.0*(gCT-1)) );
// *************** Sutherland viscosity parameters for CT
double TS0 ( ( arg0.contents( arg0.seek( "TS0" ) ) (0) ).scalar_value() );
double S1 ( ( arg0.contents( arg0.seek( "S1" ) ) (0) ).scalar_value() );
double mu0 ( ( arg0.contents( arg0.seek( "mu0" ) ) (0) ).scalar_value() );
double qdotFudge ( ( arg0.contents( arg0.seek( "qdotFudge" ) ) (0) ).\
scalar_value() ) ;
bool boolHeat ( ( bool ) ( arg0.contents( arg0.seek( "boolHeat" ) ) \
(0) ).int_value() );
bool boolLeak ( ( bool ) ( arg0.contents( arg0.seek( "boolLeak" ) ) \
(0) ).int_value() );
bool boolFriction ( ( bool )
( arg0.contents( arg0.seek( "boolFriction" ) )\
(0) ).int_value() );
bool boolPLoss ( ( bool ) ( arg0.contents( arg0.seek( "boolPLoss" ) ) \
(0) ).int_value() );
int flagPlanStop ( ( arg0.contents( \
arg0.seek( "flagPlanStop" ) ) (0) ).int_value() );
// flagPlanStop:
// = 0, no planned stop
// = 1, stop at peak
// = 2, stop at burst
// *************** Initialisation
// Numerical process
bool boolStop = 0;// Flag to stop simulation loop
bool boolBurst = 0;// Flag to indicate that primary diaphragm has burst
int valStage = 1;// used to properly create u(xdash), particle velocity
// = 1, primary + secondary pistons travelling together in CT1
// = 2, secondary piston alone travelling forward in CT2
// = 3, secondary piston rebounding from maximum compression
bool boolReachedPeak = 0;// Flag indicating whether pressure is falling
// beyond peak. It seems the primary diaphragm
// always bursts on the way down.
double pCTmax = 0.0;// Peak pCT
int indStage2 = 0;// time index at which stage 2 commences
int tNalloc ( ( arg0.contents( arg0.seek( "tNalloc" ) )
(0) ).int_value() );
ColumnVector dt( tNalloc );
double dt0 ( ( arg0.contents( arg0.seek( "dt0" ) ) (0) ).scalar_value() );
dt( 0 ) = dt0;
ColumnVector t( tNalloc );
t( 0 ) = 0.0;
int i = 0;
int N ( ( arg0.contents( arg0.seek( "N" ) ) (0) ).int_value() ) ;
// N = Number of points for heat transfer discretisation
ColumnVector dtmax( tNalloc );
dtmax( 0 ) = dt0;
// Piston variables
ColumnVector x( tNalloc );
ColumnVector v( tNalloc );
x( 0 ) = 0.0;
v( 0 ) = 0.0;
double dx = 0.0;
double dv = 0.0;
// Piston constants
double A = A1;
double mp = m1 + m2;
// Diaphragm
double ADiaphragm = 0.0;
// Volumes
double V1 = Vr0;
double V2 = VCT0;
double V120 = pi *( pow( ( 0.270 / 2.0 ), 2.0 ) - pow( 0.100, 2.0) ) *
0.5;
// Gas conditions
// Reservoir
ColumnVector pr( tNalloc );
ColumnVector Tr( tNalloc );
ColumnVector mr( tNalloc );
ColumnVector dmr( tNalloc );
pr( 0 ) = pr0;
Tr( 0 ) = T0;
mr( 0 ) = pr( 0 ) * Vr0 / ( Rr * Tr( 0 ) );
dmr( 0 ) = 0.0;
double ubar = 0.0;
// Compression tube
ColumnVector xdash( N );
ColumnVector u( N, 0.0 );// u, uA, M, Re, use vector of length <= N
ColumnVector M( N, 0.0 );
ColumnVector Re( N, 0.0 );
ColumnVector uA( N, 0.0 );
ColumnVector pCT( tNalloc );
ColumnVector dpCT( tNalloc );
ColumnVector TCT( tNalloc );
ColumnVector DTCT( tNalloc );
ColumnVector rhoCT( tNalloc );
ColumnVector mCT( tNalloc );
ColumnVector dmCT( tNalloc );
ColumnVector aCT( tNalloc );
ColumnVector dpTransition( tNalloc );
ColumnVector dpCage( tNalloc );
ColumnVector Ff( tNalloc );
ColumnVector qdot( tNalloc );
double uAp = 0.0;
int indxdashPiston = 0;
int Ndash;
// Some utility variables
// Spatial discretisation
double duAp;
double ubarSum = 0.0;
double D1, D2, x1, x2;
int jdash = 0;
// Heating
double qSum = 0.0;
double ubarIntegral;
double rhoCTEck = 0.0;
double muCTEck = 0.0;
double TCTEck = 0.0;
// Pressure loss
double Adp1, Adp2;
// Friction
// none
// Switches
bool indBurst = 0;
pCT( 0 ) = pCT0;
TCT( 0 ) = T0;
rhoCT( 0 ) = pCT( 0 ) / ( RCT *TCT( 0 ) );
mCT( 0 ) = rhoCT( 0 ) * VCT0;
aCT( 0 ) = sqrt( gCT * RCT * TCT( 0 ) );
dmCT( 0 ) = 0.0;
dpCT( 0 ) = 0.0;
dpTransition( 0 ) = 0.0;
dpCage( 0 ) = 0.0;
// Heat transfer
ColumnVector f( N );
ColumnVector qdotdx( N, 0.0 );
ColumnVector DTCT_W( tNalloc, 0.0 );
ColumnVector DTCT_Q( tNalloc, 0.0 );
ColumnVector DTCT_m( tNalloc, 0.0 );
// Pressure loss(es)
double Kt = 0.0;
// Piston ODE
ColumnVector tDE( 2 );
ColumnVector xDE0( 2 );
Matrix xDE( 2, 2);
int indMinReachedPeak = 10;// Min i before allow boolReachedPeak to be set
std::string strStatus( "OK" );
// CREATE SPATIAL DISCRETISATION AND DIAMETER PROFILE
double dxdash = ( LCT1 + LCT2 - x( 0 ) ) / ( N - 1.0 );
xdash( 0 ) = 0.0;
double min1, min2, diff1, diff2;
int indCT12 = 1;
int indCage = 1;
min1 = LCT1;
min2 = LCT1 + LCT2;
for ( int k=1 ; k < N ; k++ ) {
xdash( k ) = xdash( k - 1 ) + dxdash;
// We need the rest to find the closest discretised points to
// x = LCT1 and x = LCT1 + LCT2 - LCage
diff1 = fabs( xdash( k ) - LCT1 );
diff2 = fabs( xdash( k ) - ( LCT1 + LCT2 - LCage ) );
if ( k == 1 ) {
min1 = diff1;
min2 = diff2;
} // if ( k == 1 )
if ( diff1 < min1 ) {
min1 = diff1;
indCT12 = k;
} // if ( diff1 < min1 )
if ( diff2 < min2 ) {
min2 = diff2;
indCage = k;
} // if ( diff2 < min2 )
} // for ( int k=1 ; k < N ; k++ )
// Account for eps error in the addition process
xdash( N - 1 ) = LCT1 + LCT2;
indxdashPiston = 0;
Ndash = N;
ColumnVector D( N, -1.0 );
D( 0 ) = dCT1;
for ( int indD=1; indD < D0.rows() ; indD++ ) {
// Fill up diameter vector
D1 = D0( indD - 1, 1 );
D2 = D0( indD, 1 );
x1 = D0( indD - 1, 0 );
x2 = D0( indD, 0 );
for ( int indxdash=0; indxdash < N ; indxdash++ ) {
if ( x1 < xdash( indxdash ) && xdash( indxdash ) <= x2 ) {
if ( D0( indD, 2 ) == 0 ) // segment of constant diameter
D( indxdash ) = D2;
else if ( D0( indD, 2 ) == 1 ) // segment of linear diameter
D( indxdash ) = D1 + ( D2 - D1 ) * \
( xdash( indxdash ) - x1 ) / ( x2 - x1 );
} // if ( x1 < xdash( indxdash ) && xdash( indxdash ) <= x2 )
} // for ( int indxdash=0; indxdash < N ; indxdash++ )
} // for ( int indD=1; indD <= D0.rows ; indD++ )
/*
// Just keeping this as a code example
logfile.open ("/home/hoogland/phd/sof/octave/lib/fn/compress_test.log");
Matrix log( N, 3, 0.0 );
for ( int k=0; k < N; k++ )
log( k, 0 ) = k;
log.insert( D, 0, 1);
log.insert( xdash, 0, 2);
octave_value( log ).print( logfile , false );
logfile.close();
*/
while (!boolStop ) {
OCTAVE_QUIT; // Allow Ctrl-C halt
/* PARTICLE VELOCITY PROFILE ALONG TUBE (COMPRESSED SIDE ONLY)
First, find the 1D particle velocity along the tube
taking into account the area change. */
/* dtmax is the time taken for the return journey of a sound
wave between the piston and the primary diaphragm.
This only applies before the diaphragm has burst
(boolBurst = 0). dt0 is frozen to the value of dt at the
moment of burst. For the quasi-steady assumption to hold,
we should try and keep the time step dt below dtmax.
On the other hand, dtmax may be too large for a proper
solution of the piston ODE, so use the arbitrarily chosen
initial dt0 unless dtmax dictates a smaller value. */
dtmax( i ) = 2.0 * ( LCT1 + LCT2 - x( i ) ) / aCT( i );
dt( i ) = dt0;
if ( dtmax( i ) < dt( i ) && !boolBurst )
dt(i) = dtmax(i);
if ( v( i ) > 0.0 ) {
// rate of swept volume at piston face
uAp = v( i ) * A;
indxdashPiston = (int) ceil( x( i ) / dxdash ) - 1;
Ndash = N - indxdashPiston;
duAp = uAp / ( Ndash - 1.0 );
ubarSum = 0.0;
if ( boolHeat ) {
// Eckert reference temperature
TCTEck = 0.5 * ( TCT( i ) + Tw );
// viscosity in CT (same value throughout tube):
muCTEck = mu0 * pow( TCTEck / TS0, 1.5 ) * \
( ( TS0 + S1 ) / ( TCTEck + S1 ) );
rhoCTEck = rhoCT( i ) * TCT( i ) / TCTEck;
qSum = 0.0;
} // if ( boolHeat )
for (int j=0 ; j < Ndash ; j++ ) {
// rate of swept volume along gas in tube
jdash = j + indxdashPiston;
u( jdash ) = ( uAp - j * duAp ) \
/ ( pi * pow( D( jdash ) / 2.0, 2.0) );
// M = Mach number on compressed side
M( jdash ) = u( jdash ) / aCT( i );
// Re = Reynolds number
Re( jdash ) = rhoCTEck * D( jdash ) * u( jdash ) / muCTEck;
/*
octave_stdout << "i= " << i
<< " Ndash= " << Ndash
<< " j= " << j
<< " indxdashPiston= " << indxdashPiston
<< " jdash= " << jdash
<< "\n";
*/
if ( j > 0 ) {
ubarSum += ( xdash( jdash ) - \
xdash ( jdash - 1 ) ) * \
( u( jdash ) + u( jdash - 1 ) );
} // if j>0
// HEAT TRANSFER AT WALL FOR COMPRESSION SIDE ONLY
if ( boolHeat ) {
// Darcy-Weisbach friction factor
// No motion, Laminar, Transitional and Turbulent
if ( Re( jdash ) <= 0.001 )
f( jdash ) = 0.0; // No motion for practical purposes
else if ( Re( jdash ) > 1.0 && Re( jdash ) < 2000.0 )
f( jdash ) = 64.0 / Re( jdash );// Laminar
else if ( Re( jdash ) >= 2000.0 && Re( jdash ) <= 4000.0 )
f( jdash ) = 0.032 * pow( Re( jdash ) / 2000.0, \
0.3187 );
else if ( Re ( jdash ) > 4000.0 ) // Turbulent
f( jdash ) = pow( 1.8 * log10( Re( jdash ) ) \
- 1.5147, -2.0 );
qdotdx( jdash ) = rhoCT( i ) * cpCT * fabs( u( jdash ) ) \
* ( f( jdash ) / 8.0 ) * pow( Pr, -0.6667 ) * \
pi * D( jdash ) * ( Tw - TCT( i ) ) * \
qdotFudge;
if ( j > 0 ) {
qSum += ( xdash( jdash ) - \
xdash ( jdash - 1 ) ) \
* ( qdotdx( jdash ) + qdotdx( jdash - 1 ) );
} // if ( j > 0 )
} // if ( boolHeat )
} // for (int j=0 ; j < Ndash ; j++ )
ubarIntegral = (0.5 * ubarSum);
ubar = ubarIntegral / ( xdash( N - 1 ) - xdash
( indxdashPiston ) );
if ( boolHeat ) {
qdot( i ) = 0.5 * qSum;
DTCT_Q( i ) = dt( i ) * qdot( i ) / ( cvCT * mCT( i ) );
} // if ( boolHeat )
else
DTCT_Q( i ) = 0.0;
} // if ( v( i ) > 0.0 )
// TEMPERATURE DELTAS
// Temperature change due to mass transfer
DTCT_m( i ) = ( cvCT * TCT( i ) + 0.5 * pow( ubar, 2.0 ) ) * \
dmCT( i ) / ( cvCT * mCT( i ) );
// Temperature change due to compression work
DTCT_W( i ) = pCT( i ) * A * dx / ( cvCT * mCT( i ) );
// CT temperature change is sum of individual temperature changes
// due to (1) Work (2) Heat transfer and (3) Mass transfer.
DTCT( i ) = DTCT_W( i ) + boolHeat * DTCT_Q( i ) + boolLeak * \
DTCT_m( i );
/* PRESSURE LOSS OF COMPRESSED GAS DUE TO SUDDEN AREA CHANGE(S)
This is basically guesswork, based on available data, source:
"Pressure losses in flow through a sudden contraction of duct area",
ESDU Data Item 01016, March 2002 obtained as pdf electronically
via the UQ Cybrary. */
dpCT( i + 1 ) = 0.0;
if ( boolPLoss) {
// Deal with transition from CT1 to CT2
if ( valStage == 1 ) {
Adp1 = A1;
Adp2 = A2;
// Calculate incompressible Kt (Kti) first.
if ( Re( indCT12 ) < 0.001 ) // No motion
Kt = 0.0;
else if ( Re( indCT12 ) <= 10000.0 )
Kt = (0.32 + 159.0 / Re( indCT12 ) ) \
* ( 1.0 - pow( Adp2 / Adp1, 2.0) );
else if ( Re( indCT12 ) > 10000.0 )
Kt = 0.54 + 0.004 * ( Adp2 / Adp1 ) - \
1.285 * pow( Adp2 / Adp1, 2.0) + \
0.741 * pow( Adp2 / Adp1, 3.0);
// Transform Kti to Kt for compressible regime
if ( M( indCT12 ) > 0.3 )
Kt = Kt + 0.2 * pow( M( indCT12 ), 1.85)*\
sqrt( 1.23 - pow( M( indCT12 ), 6.3))*\
sqrt( 1.0 - pow( ( 0.6 - pow( Kt, 1.1) / 0.6 ),
2.0 ) );
dpTransition( i + 1 ) = -Kt * 0.5 * rhoCT( i ) * \
pow( ubar, 2.0);
dpCT( i + 1 ) += dpTransition( i + 1 ) * pLossFudge;
} // if ( valStage == 1 )
Adp1 = A2;
Adp2 = ACage;
// Account for area change at start of cage
// Calculate incompressible Kt (Kti) first.
if ( Re( indCage ) < 0.001 ) // No motion
Kt = 0;
else if ( Re( indCage ) <= 10000.0 )
Kt = ( 0.32 + 159.0 / Re( indCage ) ) * \
( 1.0 - pow( Adp2 / Adp1, 2.0) );
else if ( Re( indCage ) > 10000.0 )
Kt = 0.54 + 0.004 * ( Adp2 / Adp1 ) - \
1.285 * pow( Adp2 / Adp1, 2.0 ) + \
0.741 * pow( Adp2 / Adp1, 3.0 );
// Transform Kti to Kt for compressible regime
if ( M( indCage ) > 0.3 )
Kt = Kt + 0.2 * pow( M( indCage ), 1.85) * \
sqrt( 1.23 - pow( M( indCage ), 6.3) ) * \
sqrt( 1.0 - pow( ( 0.6 - pow( Kt, 1.1) / 0.6 ), 2.0 ) );
dpCage( i + 1 ) = -Kt * 0.5 * rhoCT( i )* pow( ubar, 2.0);
dpCT( i + 1 ) += dpCage( i + 1 ) * pLossFudge;
} // if ( boolPLoss)
// NEW GAS CONDITIONS
pr( i + 1 ) = mr( i ) * Rr * Tr( i ) / V1;
// Assume adiabatic reservoir expansion
Tr( i + 1 ) = Tr( i ) * pow( ( pr( i + 1 ) / pr( i ) ), \
( ( gr - 1.0 ) / gr ) );
pCT( i + 1 ) = mCT( i ) * RCT * TCT( i ) / V2;
pCT( i + 1 ) += dpCT( i + 1 );
TCT( i + 1 ) = TCT( i ) + DTCT( i );
aCT( i + 1 ) = sqrt( gCT * RCT * TCT( i + 1 ) );
// SLIDING FRICTION
// Sliding friction force Ff
if ( boolFriction ) {
if ( valStage == 1 ) // Friction proportional to the weight force
Ff( i ) = muf01 * mp * grav;
else if ( valStage == 2 ) {
// Sliding friction due to chevron seal at front of piston
// is proportional to downstream pressure.
Ff( i ) = mufc * ASeal * pCT( i ) + \
muf02 * mp * grav * mufFudge;
if ( x( i ) < ( LCT1 + 1.19 ) ) {
// Accounts for unusual friction within primary piston
// Ff( i ) = 139e3;// force due to sabot bending and
friction
} // if ( x( i ) < ( LCT1 + 1.19 ) )
} // if ( valStage == 2 )
} // if ( boolFriction )
else // boolFriction == 0
Ff( i ) = 0.0;
// PISTON MOTION ODE
// This ODE uses the adiabatic assumption which should
// hold close to the piston.
// These globals are defined for use by piston_fn_oct.oct
set_global_value( "piston_A", A );
set_global_value( "piston_p01", pr( i ) );
set_global_value( "piston_p02", pCT( i ) );
set_global_value( "piston_V01", V1 );
set_global_value( "piston_V02", V2 );
set_global_value( "piston_m", mp );
set_global_value( "piston_g1", gr );
set_global_value( "piston_g2", gCT );
set_global_value( "piston_x0", x( i ) );
set_global_value( "piston_Ff", Ff( i ) );
t( i + 1 ) = t( i ) + dt( i );
tDE( 0 ) = t( i );
tDE( 1 ) = t( i + 1 );
xDE0( 0 ) = x( i );
xDE0( 1 ) = v( i );
//[xDE,istate,msg] = lsode("piston_fn_oct",xDE0,tDE);
argout = lsode_call( xDE0, tDE );
xDE = Matrix( argout( 0 ).matrix_value() );
v( i + 1 ) = xDE( 1, 1);
x( i + 1 ) = xDE( 1, 0);
dx = x( i + 1 ) - x( i );
dv = v( i + 1 ) - v( i );
V1 += A * dx;
V2 += -A * dx;
// MASS TRANSFER
/* We only consider leakage around the secondary piston
(since the primary piston uses an o-ring seal, we assume
this is leak-proof for now). Furthermore, we assume
that the leakage acts like a simple sonic throat, which
is valid as long as the supply total pressure (pCT) is
only a small fraction greater than the external total
pressure (pr, behind the secondary piston in this case).
Finally, we assume that the V1 volume accepts the leaked
gas, instantly uniformly mixes and changes temperature
in proportion to the new mass. */
dmCT( i + 1 ) = 0.0;
dmr( i + 1 ) = 0.0;
// Chevron seal leakage
if ( boolLeak && valStage >= 2 ) {
dmCT( i + 1 ) = dmCT( i + 1 ) - \
ALeak * pCT( i + 1 ) * sqrt( gCT / ( RCT * TCT( i + 1 ) ) ) *
\
CeCT * dt( i );
dmr( i + 1 ) += -dmCT( i + 1 );
} // if ( boolLeak && valStage >= 2 )
// Mass flow through burst diaphragm
if ( boolBurst ) {
/* We progressively increase the size of ADiaphragm in a
linear fashion with time based on the opening duration
estimate DtOpening, until it reaches A2. This allows
the linear increase even if dt(i) is not constant. */
ADiaphragm += ADiaphragmOpen * dt( i ) / DtOpening;
if ( ADiaphragm > ADiaphragmOpen )
ADiaphragm = ADiaphragmOpen;
dmCT( i + 1 ) += -ADiaphragm * pCT( i + 1 ) * \
sqrt( gCT / ( RCT * TCT( i + 1 ) ) ) * CeCT * dt( i );
} // if ( boolBurst )
mCT( i + 1 ) = mCT( i ) + dmCT( i + 1 );
// Already accounting for temperature change in CT due to
// leakage in DTCT equation.
mr( i + 1 ) = mr( i ) + dmr( i + 1 );
// New temperature is weighted average of old and new parts
Tr( i + 1 ) = ( Tr( i + 1 ) * mr( i ) + TCT( i + 1 ) * dmr( i ) ) / \
mr( i + 1 );
rhoCT( i + 1 ) = mCT( i + 1 ) / V2;
// SWITCHES
/*
if ( i > 10 )
boolStop = 1;
*/
if ( i > tNalloc - 3 ) {
// We want the last element index to be iEnd = tNalloc - 1.
// When i = tNalloc - 2, pCT( i + 1 ) = pCT( iEnd ) etc have
// already been filled. So we want to stop i when it
// reaches tNalloc - 2, i.e. i > tNalloc - 3.
boolStop = 1;
strStatus = "Stopped due to i reaching tNalloc.";
} // if ( i > tNalloc - 1 )
if ( x( i + 1 ) >= LCT1 && valStage == 1 ) {
valStage = 2;
indStage2 = i;
A = A2;
mp = m2;
/* When the secondary piston is thrown from the primary
piston we have to decide what the conditions are in the
space behind the secondary. We know "it's not much"
compared with the conditions in the driver tube (CT2), so
choose the small but finite volume inside the rubber buffer,
and use the conditions in front of the primary piston. */
V1 = V120;
pr( i + 1 ) = pCT( i );
Tr( i + 1 ) = TCT( i );
mr( i + 1 ) = pr( i + 1 ) * V1 / ( RCT * Tr( i + 1 ) );
} // if ( x( i + 1 ) >= LCT1 && valStage == 1 )
if ( v( i + 1 ) <= 0.0 && valStage == 2 )
valStage = 3;// Go to stage 3 when secondary piston begins rebound
if ( valStage == 3 && ( v( i + 1 ) > 0 || x( i + 1 ) < LCT1 ) ) {
// Stop if piston begins moving forward or reaches the
// CT1-2 transition in stage 3.
boolStop = 1;
if ( v( i + 1 ) > 0 )
strStatus = "*** WARNING: Secondary piston moving" \
" forward after rebound ***";
if ( x( i + 1 ) < LCT1 )
strStatus = "Successful return of piston to CT12" \
" transition";
} // if ( valStage == 3 && ( v( i + 1 ) > 0 || x( i + 1 ) < LCT1 ) )
if ( pCT( i + 1 ) >= pCTmax )
pCTmax = pCT( i + 1 );
else if ( !boolReachedPeak && i > indMinReachedPeak \
&& valStage > 1 ) {
// It was found necessary to add the condition valStage > 1
// because occasionally little drops would occur in stage 1
// that would set boolReachedPeak. More sophisticated
// peak detection may be necessary.
boolReachedPeak = 1;
}
if ( boolReachedPeak && flagPlanStop == 1 ) {
boolStop = 1;
strStatus = "Planned stop as pCT began falling.";
} // if ( boolReachedPeak && flagPlanStop == 1 )
if ( boolReachedPeak && flagPlanStop == 2 \
&& pCTmax < pBurst ) {
// If pCT has reached peak, but flagPlanStop says to stop
// when diaphragm has burst, but diaphragm won't burst
// because peak is below pBurst. If left unchecked,
// simulation will continue until tNalloc is reached.
boolStop = 1;
strStatus = "Stop planned at burst but peak below pBurst.";
} // if ( boolReachedPeak && flagPlanStop == 2
// && pCTmax < pBurst )
if ( boolReachedPeak && !boolBurst && \
pCTmax >= pBurst && pCT( i + 1 ) <= pBurst ) {
/* Originally, the condition pCT(i+1) >= pBurst meant that
burst occurred as soon as the critical pressure was
reached. The deformation of the diaphragm, however,
usually results in burst occurring as the pressure is
falling. By creating a flag boolReachedPeak, switched to 1
when the pressure starts falling, and catching the
pressure as soon as it falls _below_ the critical level,
we can simulate this behaviour. */
boolBurst = 1;
ADiaphragm = 0.0;
indBurst = i;
// Now freeze dt0 at the current value
dt0 = dt( i );
} // if ( boolReachedPeak && !boolBurst &&
// pCTmax >= pBurst && pCT( i + 1 ) <= pBurst )
if ( boolBurst && flagPlanStop == 2 ) {
boolStop = 1;
strStatus = "Planned stop at diaphragm burst.";
} // if ( boolBurst && flagPlanStop == 2 )
if ( x( i ) >= ( LCT1 + LCT2 - LCage ) ) {
boolStop = 1;
strStatus = "*** WARNING: Piston contacted cage! ***";
} // if ( x( i ) >= ( LCT1 + LCT2 - LCage ) )
if ( pCT( i + 1 ) <= eps ) {
boolStop = 1;
strStatus = "*** WARNING: pCT dropped below eps! ***";
} // if ( pCT( i + 1 ) <= eps )
// OUTPUT
if ( verbosity >= 1 ) {
if ( i == 0 || i == 1 || i % nlog == 0 || boolStop ) {
if ( boolStop )
++i;// we want to output the very last elements
fprintf( pfile, \
">i= %6d dt= %12.4e x= %12.4e v= %12.4e"
" dx= %12.4e dv= %12.4e" \
" pCT= %12.4e TCT= %12.4e DTCT_Q= %12.4e mCT= %12.4e" \
" dpCT= %12.4e Kt= %12.4e rhoCT= %12.4e ubar= %12.4e" \
" indxdashPiston= %6d" \
" *********\n", \
i, dt( i ), x( i ), v( i ), dx, dv, pCT( i ), TCT( i ), \
DTCT_Q( i ), mCT( i ), dpCT( i ), Kt, rhoCT( i ), ubar, \
indxdashPiston );
fflush( pfile );
if ( verbosity == 2 ) {
for ( int k=0; k < N; k++ ) {
if ( k >= indxdashPiston )
fprintf( pfile, ">> " );
else
fprintf( pfile, ">>*" );
fprintf( pfile, "k= %6d xdash= %12.4e u= %12.4e "
"qdotdx= %12.4e Re= %12.4e f= %12.4e\n", \
k, xdash( k ), u( k ), qdotdx( k ), \
Re( k ), f( k ) );
fflush( pfile );
} // for ( int k=0; k < N; k++ )
} // if ( verbosity == 2 )
if ( boolStop )
--i;// but we need to shift the counter back too
} // ( i == 0 || i == 1 || i % nlog == 0 || boolStop )
} // if ( verbosity >= 1 )
++i;
} // while ( !boolStop )
t = t.extract( 0, i );
pCT = pCT.extract( 0, i );
TCT = TCT.extract( 0, i );
aCT = aCT.extract( 0, i );
x = x.extract( 0, i );
v = v.extract( 0, i );
dpCT = dpCT.extract( 0, i );
Ff = Ff.extract( 0, i );
qdot = qdot.extract( 0, i );
// additional
dpCT = dpCT.extract( 0, i );
++i; // Transform to Octave index range
arg0.assign( "indStage2", ( Cell ) octave_value( indStage2 ) );
arg0.assign( "strStatus", ( Cell ) octave_value( strStatus ) );
arg0.assign( "valStage", ( Cell ) octave_value( valStage ) );
arg0.assign( "boolBurst", ( Cell ) octave_value( boolBurst ) );
arg0.assign( "t", ( Cell ) octave_value( t ) );
arg0.assign( "pCT", ( Cell ) octave_value( pCT ) );
arg0.assign( "TCT", ( Cell ) octave_value( TCT ) );
arg0.assign( "aCT", ( Cell ) octave_value( aCT ) );
arg0.assign( "x", ( Cell ) octave_value( x ) );
arg0.assign( "v", ( Cell ) octave_value( v ) );
arg0.assign( "dpCT", ( Cell ) octave_value( dpCT ) );
arg0.assign( "Ff", ( Cell ) octave_value( Ff ) );
arg0.assign( "qdot", ( Cell ) octave_value( qdot ) );
arg0.assign( "D", ( Cell ) octave_value( D ) );
// additional
arg0.assign( "dpTransition", ( Cell ) octave_value( dpTransition ) );
retval.append( arg0 );
octave_stdout.unsetf( std::ios_base::scientific );
if ( verbosity >= 1 ) {
fclose( pfile );
logfile.close();
} // if ( verbosity >= 1 )
return retval;
} // DEFUN_DLD( compress_fn_oct , arg , ,
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------
- Best ways to free memory within c++ DLD's,
Jason Hoogland <=