From 683c0c2e2131c2668e231ee109027e13d11f87f9 Mon Sep 17 00:00:00 2001 From: "Jonathan Tu (work)" Date: Mon, 23 Jan 2017 11:56:51 -0800 Subject: [PATCH 1/6] Fixed compile error. Had to include to avoid compile error due to the use of strncopy. --- src/ScalarToTecplot.cc | 49 +++++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 24 deletions(-) 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 ); } From 3a3861db9fa4b5caa60250366b8a17f73dddda84 Mon Sep 17 00:00:00 2001 From: "Jonathan Tu (work)" Date: Mon, 23 Jan 2017 12:03:26 -0800 Subject: [PATCH 2/6] Added constructor for OutputTecplot that doesn't require TecplotAllGrids argument. --- src/OutputTecplot.cc | 30 ++++++++++++++++++------------ src/OutputTecplot.h | 20 +++++++++++++------- 2 files changed, 31 insertions(+), 19 deletions(-) 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; From 3a6d19ba0b9ffc36cd54631363e2412a34d0f435 Mon Sep 17 00:00:00 2001 From: "Jonathan Tu (work)" Date: Mon, 23 Jan 2017 12:07:57 -0800 Subject: [PATCH 3/6] Fixed compile bugs in examples caused by interface changes. NavierStokesModel now takes in BaseFlow objects rather than Flux objects. Also, the pitch and plunge examples now create their own output directories (as the Oseen example already did). --- examples/Oseen.cc | 27 +++++++++++++-------------- examples/pitching.cc | 22 +++++++++++++--------- examples/plunging.cc | 22 +++++++++++++--------- 3 files changed, 39 insertions(+), 32 deletions(-) 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..8f01da6 100644 --- a/examples/pitching.cc +++ b/examples/pitching.cc @@ -11,6 +11,7 @@ #include #include +#include #include "ibpm.h" using namespace std; @@ -22,7 +23,7 @@ int main(int argc, char* argv[]) { // lift and drag double lift = 0.; double drag = 0.; - + // Setup grid int nx = 100; int ny = 100; @@ -31,13 +32,13 @@ int main(int argc, char* argv[]) { 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; @@ -51,7 +52,7 @@ int main(int argc, char* argv[]) { 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(); @@ -66,16 +67,19 @@ int main(int argc, char* argv[]) { State x(grid, geom.getNumPoints()); x.omega = 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; diff --git a/examples/plunging.cc b/examples/plunging.cc index 88d32f2..395f564 100644 --- a/examples/plunging.cc +++ b/examples/plunging.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include #include "ibpm.h" @@ -24,7 +25,7 @@ int main(int argc, char* argv[]) { // lift and drag double lift = 0.; double drag = 0.; - + // Setup grid int nx = 200; int ny = 200; @@ -33,13 +34,13 @@ int main(int argc, char* argv[]) { 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; @@ -53,7 +54,7 @@ int main(int argc, char* argv[]) { 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(); @@ -68,16 +69,19 @@ int main(int argc, char* argv[]) { State x(grid, geom.getNumPoints()); x.omega = 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( &force, 1 ); logger.init(); logger.doOutput( x ); - + // Step const double PI = 4. * atan(1.); int numSteps = 250; From 644f356f5f15ddf8cfbc50432f95f3602802260f Mon Sep 17 00:00:00 2001 From: "Jonathan Tu (work)" Date: Mon, 23 Jan 2017 12:31:28 -0800 Subject: [PATCH 4/6] Ignore Windows executables. Also, ignore example output files by directory names rather than extensions. --- .gitignore | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) 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 From 72f390c9a56109b3828648d20a5074add7de2e22 Mon Sep 17 00:00:00 2001 From: "Jonathan Tu (work)" Date: Fri, 24 Feb 2017 08:33:33 -0800 Subject: [PATCH 5/6] Fixed pitching and plunging examples by correctly initializing state to zero. Previously, vorticity had been set to zero, but boundary forces and fluxes were not specified. This meant the associated arrays were full of random data, leading to nondeterministic and incorrect behavior. --- examples/pitching.cc | 11 ++++++----- examples/plunging.cc | 23 ++++++++++++----------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/examples/pitching.cc b/examples/pitching.cc index 8f01da6..25be0bd 100644 --- a/examples/pitching.cc +++ b/examples/pitching.cc @@ -28,12 +28,11 @@ int main(int argc, char* argv[]) { 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() ); @@ -46,10 +45,10 @@ 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 BaseFlow q_potential( grid, magnitude, alpha ); @@ -66,6 +65,8 @@ 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 ); @@ -89,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 395f564..885478f 100644 --- a/examples/plunging.cc +++ b/examples/plunging.cc @@ -27,31 +27,30 @@ int main(int argc, char* argv[]) { 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 BaseFlow q_potential( grid, magnitude, alpha ); @@ -61,13 +60,15 @@ int main(int argc, char* argv[]) { 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 ); @@ -77,7 +78,7 @@ int main(int argc, char* argv[]) { OutputForce force( "plunging_out/force.dat" ); Logger logger; // Output Tecplot file every few timesteps - logger.addOutput( &tecplot, 10 ); + logger.addOutput( &tecplot, 25 ); logger.addOutput( &force, 1 ); logger.init(); logger.doOutput( x ); From cd2f5565cef75cebc11098c36e80165ad43284da Mon Sep 17 00:00:00 2001 From: "Jonathan Tu (work)" Date: Fri, 24 Feb 2017 09:04:16 -0800 Subject: [PATCH 6/6] Fixed definition of pi. Was previously entered as a long decimal. Now uses arctan, as is done in other files. --- src/RigidBody.cc | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) 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);