diff --git a/.gitignore b/.gitignore index df53fde..e3bd430 100644 --- a/.gitignore +++ b/.gitignore @@ -2,12 +2,10 @@ build/*.o build/libibpm.a build/checkgeom build/ibpm -examples/*.cmd -examples/*.force -examples/*.bin -examples/*.plt -examples/*.cholesky -TAGS +examples/oseen_out +examples/pitching_out +examples/plunging_out +plunging_outTAGS test/*.o test/runner test/runner.err @@ -30,3 +28,4 @@ doc/examples/*.out doc/snippets/*.out callgrind.out.* .depend +*.exe diff --git a/examples/Oseen.cc b/examples/Oseen.cc index 859cd3d..e48d87c 100644 --- a/examples/Oseen.cc +++ b/examples/Oseen.cc @@ -50,7 +50,7 @@ const double alpha = 1.256431208626170; int main(int argc, char* argv[]) { cout << "Test of Oseen vortex\n"; - + // Setup grid int nx = 100; int ny = 100; @@ -59,7 +59,7 @@ int main(int argc, char* argv[]) { double xOffset = -5; double yOffset = -5; Grid grid( nx, ny, ngrid, length, xOffset, yOffset ); - + // Empty geometry -- no body Geometry geom; int numPoints = 0; @@ -67,12 +67,11 @@ int main(int argc, char* argv[]) { // Setup equations to solve double Reynolds=100; // No background flow - Flux q0( grid ); - q0 = 0; + BaseFlow q0( grid ); NavierStokesModel model( grid, geom, Reynolds, q0 ); model.init(); - + // Setup timestepper double dt = 0.05; NonlinearIBSolver solver( grid, model, dt, Scheme::EULER ); @@ -87,22 +86,22 @@ int main(int argc, char* argv[]) { // Create output directory, if does not already exist mkdir( "oseen_out", S_IRWXU | S_IRWXG | S_IRWXO ); - + // Setup output routines OutputTecplot outputComputed( "oseen_out/ibpm%03d.plt", - "Oseen vortex, numerical, step %03d" ); + "Oseen vortex, numerical, step %03d" ); OutputTecplot outputExact( "oseen_out/exact%03d.plt", - "Oseen vortex, exact, step %03d"); + "Oseen vortex, exact, step %03d" ); OutputTecplot outputError( "oseen_out/error%03d.plt", - "Oseen vortex, error, step %03d"); - + "Oseen vortex, error, step %03d" ); + // Output initial condition outputComputed.doOutput( x ); outputExact.doOutput( exact ); outputError.doOutput( error ); - // Step - int numSteps = 200; + // Step + int numSteps = 200; int iskip = 20; for(int i=1; i <= numSteps; ++i) { cout << "step " << i << endl; @@ -183,7 +182,7 @@ void initializeOseenVortex( double Reynolds, State& x ) { - + double t0 = Reynolds / (4 * alpha); - computeExactSolution( Reynolds, t0, 0, x ); + computeExactSolution( Reynolds, t0, 0, x ); } diff --git a/examples/pitching.cc b/examples/pitching.cc index aa092bb..25be0bd 100644 --- a/examples/pitching.cc +++ b/examples/pitching.cc @@ -11,6 +11,7 @@ #include #include +#include #include "ibpm.h" using namespace std; @@ -22,22 +23,21 @@ int main(int argc, char* argv[]) { // lift and drag double lift = 0.; double drag = 0.; - + // Setup grid int nx = 100; int ny = 100; int ngrid = 1; - double length = 4.0; + double length = 4; double xOffset = -1; double yOffset = -2; Grid grid( nx, ny, ngrid, length, xOffset, yOffset ); - - + // Make a flat plate, length 1, with center at 1/4 chord RigidBody plate; plate.addLine( 0, 0, 1, 0, grid.Dx() ); plate.setCenter( 0.25, 0 ); - + // Set the motion to pitching, amplitude = 0.25, period 10 time units double amplitude = 0.25; double freq = 0.1; @@ -45,13 +45,13 @@ int main(int argc, char* argv[]) { plate.setMotion( motion ); Geometry geom; geom.addBody( plate ); - geom.moveBodies(0); + geom.moveBodies( 0 ); // Setup equations to solve - double Reynolds=100; + double Reynolds = 100; double magnitude = 1; double alpha = 0; // angle of background flow - Flux q_potential = Flux::UniformFlow( grid, magnitude, alpha ); + BaseFlow q_potential( grid, magnitude, alpha ); cout << "Setting up Navier Stokes model..." << flush; NavierStokesModel model( grid, geom, Reynolds, q_potential ); model.init(); @@ -65,17 +65,22 @@ int main(int argc, char* argv[]) { // Build the state variable, zero initial conditions State x(grid, geom.getNumPoints()); x.omega = 0.; + x.f = 0.; + x.q = 0.; + + // Create output directory, if does not already exist + mkdir( "pitching_out", S_IRWXU | S_IRWXG | S_IRWXO ); // Setup output routines - OutputForce force( "tecplot/force.dat" ); - OutputTecplot tecplot( "tecplot/pitch%03d.plt", "Pitching plate, step %03d" ); + OutputForce force( "pitching_out/force.dat" ); + OutputTecplot tecplot( "pitching_out/pitch%03d.plt", "Pitching plate, step %03d" ); Logger logger; // Output Tecplot file every few timesteps logger.addOutput( &tecplot, 25 ); - logger.addOutput( &force, 1 ); + logger.addOutput( &force, 1 ); logger.init(); logger.doOutput( x ); - + // Step const double PI = 4. * atan(1.); int numSteps = 250; @@ -85,7 +90,7 @@ int main(int argc, char* argv[]) { << " time = " << setw(5) << x.time << " theta = " << theta << endl; solver.advance( x ); - x.computeNetForce( drag, lift); + x.computeNetForce( drag, lift ); cout << " x force : " << setw(16) << drag*2 << " , y force : " << setw(16) << lift*2 << "\n"; logger.doOutput( x ); diff --git a/examples/plunging.cc b/examples/plunging.cc index 88d32f2..885478f 100644 --- a/examples/plunging.cc +++ b/examples/plunging.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include #include "ibpm.h" @@ -24,60 +25,64 @@ int main(int argc, char* argv[]) { // lift and drag double lift = 0.; double drag = 0.; - + // Setup grid - int nx = 200; - int ny = 200; + int nx = 100; + int ny = 100; int ngrid = 1; - double length = 4.0; + double length = 4; double xOffset = -1; double yOffset = -2; Grid grid( nx, ny, ngrid, length, xOffset, yOffset ); - - + // Make a flat plate, length 1, with center at 1/4 chord RigidBody plate; plate.addLine( 0, 0, 1, 0, grid.Dx() ); plate.setCenter( 0.25, 0 ); - - // Set the motion to plunging: amplitude = 0.1, period 0.25 time unit - double amplitude = 0.1; - double freq = 0.25; + + // Set the motion to plunging: amplitude = 0.25, period 10 time units + double amplitude = 0.25; + double freq = 0.1; PitchPlunge motion( 0, 0, amplitude, freq ); plate.setMotion( motion ); Geometry geom; geom.addBody( plate ); - geom.moveBodies(0); + geom.moveBodies( 0 ); // Setup equations to solve - double Reynolds=100; + double Reynolds = 100; double magnitude = 1; double alpha = 0; // angle of background flow - Flux q_potential = Flux::UniformFlow( grid, magnitude, alpha ); + BaseFlow q_potential( grid, magnitude, alpha ); cout << "Setting up Navier Stokes model..." << flush; NavierStokesModel model( grid, geom, Reynolds, q_potential ); model.init(); cout << "done" << endl; // Setup timestepper - double dt = 0.001; + double dt = 0.005; NonlinearIBSolver solver( grid, model, dt, Scheme::AB2 ); solver.init(); // Build the state variable, zero initial conditions State x(grid, geom.getNumPoints()); x.omega = 0.; + x.f = 0.; + x.q = 0; + + // Create output directory, if does not already exist + mkdir( "plunging_out", S_IRWXU | S_IRWXG | S_IRWXO ); // Setup output routines - OutputTecplot tecplot( "tecplot/plunge%03d.plt", "Plunging plate, step %03d" ); - OutputForce force( "tecplot/force.dat" ); + OutputTecplot tecplot( "plunging_out/plunge%03d.plt", "Plunging plate, step %03d" ); + OutputForce force( "plunging_out/force.dat" ); Logger logger; // Output Tecplot file every few timesteps - logger.addOutput( &tecplot, 10 ); - logger.addOutput( &force, 1 ); + logger.addOutput( &tecplot, 25 ); + logger.addOutput( &force, 1 ); logger.init(); logger.doOutput( x ); - + // Step const double PI = 4. * atan(1.); int numSteps = 250; diff --git a/src/OutputTecplot.cc b/src/OutputTecplot.cc index bd77362..576f026 100644 --- a/src/OutputTecplot.cc +++ b/src/OutputTecplot.cc @@ -26,37 +26,43 @@ using namespace std; namespace ibpm { +OutputTecplot::OutputTecplot( string filename, string title ) { + _filename = filename; + _title = title; + _TecplotAllGrids = false; +} + OutputTecplot::OutputTecplot( string filename, string title, bool TecplotAllGrids ) { _filename = filename; _title = title; _TecplotAllGrids = TecplotAllGrids; } - + bool OutputTecplot::doOutput(const State& state) { // Add timestep to filename and title char filename[256]; sprintf( filename, _filename.c_str(), state.timestep ); char title[256]; sprintf( title, _title.c_str(), state.timestep ); - bool status = false; - const Grid& grid = state.omega.getGrid(); - + bool status = false; + const Grid& grid = state.omega.getGrid(); + // Calculate velocities Scalar u( state.omega.getGrid() ); Scalar v( state.omega.getGrid() ); FluxToVelocity( state.q, u, v ); - + // Create vector of Scalar fields vector varVec; varVec.push_back( &u ); varVec.push_back( &v); varVec.push_back( &state.omega ); - + vector varNameVec; varNameVec.push_back( "u" ); varNameVec.push_back( "v" ); varNameVec.push_back( "Vorticity" ); - + // Write the tecplot file if(_TecplotAllGrids) { status = true; @@ -69,19 +75,19 @@ bool OutputTecplot::doOutput(const State& state) { else status = ScalarToTecplot( varVec, varNameVec, filename, title ); return status; } - + bool OutputTecplot::doOutput(const BaseFlow& q, const State& x) { - // Currently no use for baseflow, but this method is defined for future + // Currently no use for baseflow, but this method is defined for future // flexibility return doOutput(x); } - + void OutputTecplot::setFilename( string filename ) { _filename = filename; } - + void OutputTecplot::setTitle( string title ) { _title = title; } - + } // namespace ibpm diff --git a/src/OutputTecplot.h b/src/OutputTecplot.h index 99bd6dd..d252d89 100644 --- a/src/OutputTecplot.h +++ b/src/OutputTecplot.h @@ -9,7 +9,7 @@ using std::string; namespace ibpm { - + /*! \file OutputTecplot.h \class OutputTecplot @@ -24,26 +24,32 @@ namespace ibpm { \version $Revision$ */ - + class OutputTecplot : public Output { public: /// \brief Constructor /// \param[in] filename Filename in the standard printf format (e.g. "file%06d.plt"), where timestep will be supplied /// \param[in] title Title in the standard printf format + OutputTecplot( string filename, string title ); + + /// \brief Constructor + /// \param[in] filename Filename in the standard printf format (e.g. "file%06d.plt"), where timestep will be supplied + /// \param[in] title Title in the standard printf format + /// \param[in] TecplotAllGrids bool describing whether or not to save Tecplot data for all grid levels OutputTecplot( string filename, string title, bool TecplotAllGrids ); - + /// \brief Write the Tecplot file bool doOutput(const State& x); - + /// \brief Write the Tecplot file bool doOutput(const BaseFlow& q, const State& x); - + /// \brief Change the filename for the output file void setFilename( string filename ); - + /// \brief Change the title for the output file void setTitle( string title ); - + private: string _filename; string _title; diff --git a/src/RigidBody.cc b/src/RigidBody.cc index d1274ea..4145a43 100644 --- a/src/RigidBody.cc +++ b/src/RigidBody.cc @@ -5,9 +5,9 @@ // // The points on the body are stored with respect to a reference configuration // in _refPoints -// -// The current locations of the points on the body, defined by the associated -// Motion, are contained in _currentPoints, and are updated whenever +// +// The current locations of the points on the body, defined by the associated +// Motion, are contained in _currentPoints, and are updated whenever // moveBodies() is called. // // Author(s): @@ -111,7 +111,7 @@ void RigidBody::addCircle_n( double x = xc + radius * cos( i * dTheta ); double y = yc + radius * sin( i * dTheta ); addPoint( x, y ); - } + } } void RigidBody::addLine( @@ -141,7 +141,7 @@ void RigidBody::addLine_n( double y = y1 + i * deltaY; addPoint(x,y); } -} +} void RigidBody::addLine_aoa( double l, @@ -151,8 +151,8 @@ void RigidBody::addLine_aoa( int numPoints ) { double x0, y0, x0r, y0r, x1r, y1r, x1, y1; - double pi = 3.141592653589793238462643383279502884197169399375; - double alpha = aoa / 180 * pi; + double pi = 4. * atan(1.); + double alpha = aoa / 180 * pi; double cosa = cos(alpha); double sina = sin(alpha); double delta = l / (numPoints - 1); @@ -167,7 +167,7 @@ void RigidBody::addLine_aoa( y1 = y1r+yC; addPoint(x1,y1); } -} +} int RigidBody::getNumPoints() const { return _refPoints.size(); @@ -178,7 +178,7 @@ void RigidBody::saveRaw(ostream& out) { out << n; for(int i=0; i> x >> y; if ( in.fail() ) { return false; } - addPoint(x,y); + addPoint(x,y); } return true; } @@ -339,7 +339,7 @@ bool RigidBody::load(istream& in) { addLine_aoa( l, xC, yC, aoa, numPoints ); setCenter( xC, yC ); #ifdef DEBUG - cerr << "Add a line: length" << l << ", AoA" << aoa + cerr << "Add a line: length" << l << ", AoA" << aoa ", n = " << numPoints << endl; #endif } @@ -386,7 +386,7 @@ bool RigidBody::load(istream& in) { one_line >> AMP >> DUR >> startTime; RB_CHECK_FOR_ERRORS; Motion* m = new SigmoidalStep( AMP, DUR, startTime); - setMotion( *m ); + setMotion( *m ); } else if ( motionType == "lagstep1" ) { // First-order lag filtered square pulse @@ -494,7 +494,7 @@ bool RigidBody::load(istream& in) { RB_CHECK_FOR_ERRORS; addPoint( x, y ); #ifdef DEBUG - cerr << "Add a point: (" << x << ", " << y << ")" << endl; + cerr << "Add a point: (" << x << ", " << y << ")" << endl; #endif } else if ( cmd == "raw" ) { @@ -503,7 +503,7 @@ bool RigidBody::load(istream& in) { RB_CHECK_FOR_ERRORS; loadRaw( filename ); #ifdef DEBUG - cerr << "Read a raw file: " << filename << endl; + cerr << "Read a raw file: " << filename << endl; #endif } else { @@ -519,7 +519,7 @@ bool RigidBody::load(istream& in) { string RigidBody::getName() { return _name; } - + BoundaryVector RigidBody::toBoundaryVector(const vector list) { int n = list.size(); BoundaryVector BVList(n); diff --git a/src/ScalarToTecplot.cc b/src/ScalarToTecplot.cc index 0bbbb9b..774a53f 100644 --- a/src/ScalarToTecplot.cc +++ b/src/ScalarToTecplot.cc @@ -2,12 +2,13 @@ #include "ScalarToTecplot.h" #include #include +#include #include using namespace std; namespace ibpm { - + class VarList { public: void addVariable( const Scalar* var, string varName ) { @@ -17,34 +18,34 @@ class VarList { _vars.push_back( var ); _names.push_back( varName ); } - + int getNumVars() const { return _vars.size(); } - + string getName(int i) const { return _names[i]; } - + const Scalar* getVariable(int i) const { return _vars[i]; } - + private: vector _vars; vector _names; }; - + bool writeTecplotFileASCII( const char* filename, const char* title, const VarList& list, int lev ) { int numVars = list.getNumVars(); assert( numVars > 0 ); - + // Get grid information const Grid& grid = list.getVariable(0)->getGrid(); assert( lev < grid.Ngrid() ); int nx = grid.Nx(); int ny = grid.Ny(); - + // Write the header for the Tecplot file cerr << "Writing Tecplot file " << filename << endl; FILE *fp = fopen( filename, "w" ); @@ -63,7 +64,7 @@ bool writeTecplotFileASCII( const char* filename, const char* title, const VarLi fprintf( fp, "SINGLE "); } fprintf( fp, ")\n" ); - + // Write the data for (int j=1; j varVec, vector varNameVec, string filename, string title, int lev ) { assert( varVec.size() > 0 ); assert( varVec.size() == varNameVec.size() ); - + // Get grid dimensions const Grid& grid = varVec[0]->getGrid(); int nx = grid.Nx(); int ny = grid.Ny(); int ngrid = grid.Ngrid(); assert( lev < ngrid ); - + // Calculate the variables for output // Calculate the grid Scalar x(grid); @@ -102,7 +103,7 @@ bool ScalarToTecplot( vector varVec, vector varNameVec, s } } } - + // Store pointers to variables and corresponding names in vectors VarList list; list.addVariable( &x, "x" ); @@ -110,33 +111,33 @@ bool ScalarToTecplot( vector varVec, vector varNameVec, s for( unsigned int i = 0; i < varVec.size(); i++ ) { list.addVariable( varVec[i], varNameVec[i] ); } - + // Add timestep to filename and title char _filename[BUFSIZ]; strncpy( _filename, filename.c_str(), BUFSIZ-1 ); char _title[BUFSIZ]; strncpy( _title, title.c_str(), BUFSIZ-1 ); - + // Write the Tecplot file bool status = writeTecplotFileASCII( _filename, _title, list, lev ); return status; } - + bool ScalarToTecplot( const Scalar* var, string varName, string filename, string title, int lev ) { vector varVec; varVec.push_back( var ); - + vector varNameVec; varNameVec.push_back( varName ); - + bool status = ScalarToTecplot( varVec, varNameVec, filename, title, lev); return status; } - -bool ScalarToTecplot( vector varVec, vector varNameVec, string filename, string title ) { + +bool ScalarToTecplot( vector varVec, vector varNameVec, string filename, string title ) { return ScalarToTecplot( varVec, varNameVec, filename, title, 0 ); } - + bool ScalarToTecplot( const Scalar* var, string varName, string filename, string title ) { return ScalarToTecplot( var, varName, filename, title, 0 ); }