diff --git a/.gitignore b/.gitignore index a2a97e8..e199fab 100644 --- a/.gitignore +++ b/.gitignore @@ -40,3 +40,6 @@ *.aux CMakeInit.txt Session.vim +*.dat +*.png +build/ diff --git a/apps/NavierStokes_2D/src/App.cpp b/apps/NavierStokes_2D/src/App.cpp index a6246c7..25d53c3 100644 --- a/apps/NavierStokes_2D/src/App.cpp +++ b/apps/NavierStokes_2D/src/App.cpp @@ -1,14 +1,19 @@ #include "VPM2d.hpp" #include "Point2d.hpp" #include "InitialConditions.hpp" +#include "Split_Advection.hpp" #include "Parameters.hpp" #include "Structure_Flat.hpp" #include "Structure_Hilly.hpp" #include "Structure_Circle.hpp" +#include "Structure_Wing.hpp" #include "Structure_HalfCircle.hpp" #include "Structure_InverseCircle.hpp" #include "Structure_Ellipse.hpp" + #include +#define OMPI_SKIP_MPICXX 1 +#include int m_numberofParticlesx = 50; int m_numberofParticlesy = 50; @@ -16,11 +21,12 @@ int m_numberofParticlesy = 50; int main(int argc, char** argv) { + MPI_Init(&argc, &argv); std::string inputFile; std::string outputFile; double time = 0.1; - double dt = 0.01; double nu = 0.; + double Re = 0.; double ax = 1.0; double ay = 1.0; @@ -35,11 +41,7 @@ int main(int argc, char** argv) int order = 1; - int example_num = -1; int structure_num = -1; - double example_dist = 0.5; - double example_strength = 1.; - double example_core_radius = 0.15; double population_threshold = -1; @@ -76,12 +78,12 @@ int main(int argc, char** argv) nu = std::atof(argv[i+1]); eat = 2; } - if( arg == "--T" && (i+1) < argc ) { - time = std::atof(argv[i+1]); + if( arg == "--Re" && (i+1) < argc ) { + Re = std::atof(argv[i+1]); eat = 2; } - if( arg == "--dt" && (i+1) < argc ) { - dt = std::atof(argv[i+1]); + if( arg == "--T" && (i+1) < argc ) { + time = std::atof(argv[i+1]); eat = 2; } if( arg == "--domain" && (i+4) < argc ) { @@ -105,26 +107,10 @@ int main(int argc, char** argv) order = std::atoi(argv[i+1]); eat = 2; } - if( arg == "--example_num" && (i+1) < argc ) { - example_num = std::atoi(argv[i+1]); - eat = 2; - } if( arg == "--structure_num" && (i+1) < argc ) { structure_num = std::atoi(argv[i+1]); eat = 2; } - if( arg == "--example_dist" && (i+1) < argc ) { - example_dist = std::atof(argv[i+1]); - eat = 2; - } - if( arg == "--example_strength" && (i+1) < argc ) { - example_strength = std::atof(argv[i+1]); - eat = 2; - } - if( arg == "--example_core_radius" && (i+1) < argc ) { - example_core_radius = std::atof(argv[i+1]); - eat = 2; - } if( arg == "--bc_xl" && (i+2) < argc ) { bc_xl = argv[i+1]; bc_to_xl = std::atof(argv[i+2]); @@ -156,6 +142,7 @@ int main(int argc, char** argv) } } + VPM::VPM2d* vpm = new VPM::VPM2d(argc, argv); switch (structure_num) @@ -196,6 +183,24 @@ int main(int argc, char** argv) vpm->setStructure(structure); } break; + case (6): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + } + + if (Re>0) + { + double charlength; + bool structureset; + vpm->getCharacteristicLength(charlength, structureset); + if (structureset) + { + nu = charlength*Uinfty.x/Re; + std::cerr<<"Setting nu="< remeshParams = std::make_shared( @@ -215,13 +220,12 @@ int main(int argc, char** argv) bc_to_yr ); - std::shared_ptr params = std::make_shared( + VPM::Parameters params = VPM::Parameters( domain_ll, domain_ur, m_numberofParticlesx, m_numberofParticlesy, nu, - 0, population_threshold, remeshParams, bcParams, @@ -232,20 +236,33 @@ int main(int argc, char** argv) std::vector positions; std::vector omega; VPM::ParticleField pf; - int fn_count; - bool save_init; + int fn_count = 0; + bool save_init = true; if (inputFile.empty()) { - init(*params, example_num, example_dist, example_strength, example_core_radius, positions, omega); - pf.positions = positions; - fn_count = 0; - save_init = true; + std::vector positions; + std::vector omega; + init(params, -1, -1, -1, -1, positions, omega); + + pf.omega = omega; + pf.positions = positions; + pf.regular_positions = positions; + pf.params = params; + pf.cartesianGrid = true; + pf.velocity_correspondsTo_omega = false; + + VPM::Split_Advection advection; + advection.calculateVelocity(pf); + + // it is important to set this, because a structure can/will be set + pf.velocity_correspondsTo_omega = false; } else { - - readParticlesFromFile(inputFile, pf); + VPM::ParticleField pf; + bool random_velocity_dist=false; + readParticlesFromFile(inputFile, pf, random_velocity_dist);//params->m_N, positions, omega); std::size_t found = inputFile.find_last_of("_"); std::string num = inputFile.substr(found+2, inputFile.size()); @@ -254,11 +271,9 @@ int main(int argc, char** argv) } + vpm->run(pf, time, outputFile, fn_count, save_init, false); - vpm->run(pf, omega, time, dt, params, - outputFile, fn_count, save_init - ); - + MPI_Finalize(); return 0; } diff --git a/libs/NavierStokes_VPM_2D/include/Point2d.hpp b/libs/NavierStokes_VPM_2D/include/Point2d.hpp index 147e01e..d3b24f6 100644 --- a/libs/NavierStokes_VPM_2D/include/Point2d.hpp +++ b/libs/NavierStokes_VPM_2D/include/Point2d.hpp @@ -16,7 +16,7 @@ namespace VPM int x; int y; - IPoint2d( ) {} + IPoint2d( ): x(0), y(0) {} IPoint2d( const int xy_ ) : @@ -36,7 +36,7 @@ namespace VPM double x; double y; - Point2d( ) {} + Point2d( ): x(0.), y(0.) {} Point2d( const double xy_ ) : diff --git a/libs/NavierStokes_VPM_2D/include/Structure.hpp b/libs/NavierStokes_VPM_2D/include/Structure.hpp index 91c1147..c4c2da9 100644 --- a/libs/NavierStokes_VPM_2D/include/Structure.hpp +++ b/libs/NavierStokes_VPM_2D/include/Structure.hpp @@ -9,6 +9,8 @@ namespace VPM Structure() {}; ~Structure() {}; virtual bool isInside(const Point2d pos, const double pad) = 0; + virtual void getOrigo(Point2d & pos) = 0; + virtual double getCharacteristicLength() = 0; private: diff --git a/libs/NavierStokes_VPM_2D/include/Structure_Circle.hpp b/libs/NavierStokes_VPM_2D/include/Structure_Circle.hpp index 7c057d1..1073d50 100644 --- a/libs/NavierStokes_VPM_2D/include/Structure_Circle.hpp +++ b/libs/NavierStokes_VPM_2D/include/Structure_Circle.hpp @@ -9,8 +9,12 @@ namespace VPM Structure_Circle(); ~Structure_Circle(); bool isInside(const Point2d pos, const double pad); + void getOrigo(Point2d & pos); + double getCharacteristicLength(); private: + Point2d m_origo; + double m_radius; }; diff --git a/libs/NavierStokes_VPM_2D/include/Structure_Ellipse.hpp b/libs/NavierStokes_VPM_2D/include/Structure_Ellipse.hpp index 02e1308..c072d5d 100644 --- a/libs/NavierStokes_VPM_2D/include/Structure_Ellipse.hpp +++ b/libs/NavierStokes_VPM_2D/include/Structure_Ellipse.hpp @@ -11,6 +11,8 @@ namespace VPM Structure_Ellipse(Point2d origo, double semi_major_axis, double semi_minor_axis); ~Structure_Ellipse(); bool isInside(const Point2d pos, const double pad); + void getOrigo(Point2d & pos); + double getCharacteristicLength(); private: diff --git a/libs/NavierStokes_VPM_2D/include/Structure_Flat.hpp b/libs/NavierStokes_VPM_2D/include/Structure_Flat.hpp index 71dda6c..e39b6fd 100644 --- a/libs/NavierStokes_VPM_2D/include/Structure_Flat.hpp +++ b/libs/NavierStokes_VPM_2D/include/Structure_Flat.hpp @@ -9,6 +9,8 @@ namespace VPM Structure_Flat(); ~Structure_Flat(); bool isInside(const Point2d pos, const double pad); + void getOrigo(Point2d & pos); + double getCharacteristicLength(); private: diff --git a/libs/NavierStokes_VPM_2D/include/Structure_HalfCircle.hpp b/libs/NavierStokes_VPM_2D/include/Structure_HalfCircle.hpp index eb91300..432d2ac 100644 --- a/libs/NavierStokes_VPM_2D/include/Structure_HalfCircle.hpp +++ b/libs/NavierStokes_VPM_2D/include/Structure_HalfCircle.hpp @@ -9,6 +9,8 @@ namespace VPM Structure_HalfCircle(); ~Structure_HalfCircle(); bool isInside(const Point2d pos, const double pad); + void getOrigo(Point2d & pos); + double getCharacteristicLength(); private: diff --git a/libs/NavierStokes_VPM_2D/include/Structure_Hilly.hpp b/libs/NavierStokes_VPM_2D/include/Structure_Hilly.hpp index b237a06..b817a47 100644 --- a/libs/NavierStokes_VPM_2D/include/Structure_Hilly.hpp +++ b/libs/NavierStokes_VPM_2D/include/Structure_Hilly.hpp @@ -9,6 +9,8 @@ namespace VPM Structure_Hilly(); ~Structure_Hilly(); bool isInside(const Point2d pos, const double pad); + void getOrigo(Point2d & pos); + double getCharacteristicLength(); private: diff --git a/libs/NavierStokes_VPM_2D/include/Structure_InverseCircle.hpp b/libs/NavierStokes_VPM_2D/include/Structure_InverseCircle.hpp index 934a964..ecc6c86 100644 --- a/libs/NavierStokes_VPM_2D/include/Structure_InverseCircle.hpp +++ b/libs/NavierStokes_VPM_2D/include/Structure_InverseCircle.hpp @@ -9,6 +9,8 @@ namespace VPM Structure_InverseCircle(); ~Structure_InverseCircle(); bool isInside(const Point2d pos, const double pad); + void getOrigo(Point2d & pos); + double getCharacteristicLength(); private: diff --git a/libs/NavierStokes_VPM_2D/include/Structure_Wing.hpp b/libs/NavierStokes_VPM_2D/include/Structure_Wing.hpp new file mode 100644 index 0000000..c4731ca --- /dev/null +++ b/libs/NavierStokes_VPM_2D/include/Structure_Wing.hpp @@ -0,0 +1,22 @@ +#pragma once +#include "Structure.hpp" + +namespace VPM +{ + class Structure_Wing : public Structure + { + public: + Structure_Wing(); + ~Structure_Wing(); + bool isInside(const Point2d pos, const double pad); + void getOrigo(Point2d & pos); + double getCharacteristicLength(); + + private: + Point2d m_origo; + double m_charlength; + std::vector m_p; + + }; + +} diff --git a/libs/NavierStokes_VPM_2D/include/VPM2d.hpp b/libs/NavierStokes_VPM_2D/include/VPM2d.hpp index 0316336..a8e7955 100644 --- a/libs/NavierStokes_VPM_2D/include/VPM2d.hpp +++ b/libs/NavierStokes_VPM_2D/include/VPM2d.hpp @@ -34,6 +34,8 @@ namespace VPM const bool onestep = false ); void setStructure(std::shared_ptr structure); + void getCharacteristicLength(double & charlength, bool & structureset); + private: @@ -43,13 +45,24 @@ namespace VPM double & delta_t ); + Point2d getForces( + const ParticleField & pf, + const ParticleField & pf_old, + const double & delta_t + ); + + int m_dim; std::shared_ptr m_diffusion; std::shared_ptr m_advection; std::shared_ptr m_source; + bool m_structure_set; + std::shared_ptr m_structure; + ParticleField m_pf; + ParticleField m_pf_old; std::string m_outputFile; int m_filename_count; diff --git a/libs/NavierStokes_VPM_2D/src/Split_Diffusion.cpp b/libs/NavierStokes_VPM_2D/src/Split_Diffusion.cpp index b55d6a7..45d7080 100644 --- a/libs/NavierStokes_VPM_2D/src/Split_Diffusion.cpp +++ b/libs/NavierStokes_VPM_2D/src/Split_Diffusion.cpp @@ -53,6 +53,8 @@ namespace VPM const double dt ) { + std::cerr<<" ---> Split_Diffusion:: RK2_step not yet implemented\n"; + exit(0); return; } diff --git a/libs/NavierStokes_VPM_2D/src/Structure_Circle.cpp b/libs/NavierStokes_VPM_2D/src/Structure_Circle.cpp index 5f50650..a42b1b9 100644 --- a/libs/NavierStokes_VPM_2D/src/Structure_Circle.cpp +++ b/libs/NavierStokes_VPM_2D/src/Structure_Circle.cpp @@ -5,6 +5,8 @@ namespace VPM Structure_Circle::Structure_Circle() { + m_origo = Point2d(-.5, .111); + m_radius = .2; } Structure_Circle::~Structure_Circle() { @@ -12,15 +14,24 @@ namespace VPM bool Structure_Circle::isInside(const Point2d pos, const double pad) { - double circle_r = .2; - Point2d circle_m = Point2d(-.5, .111); bool ret = false; - if (L2Norm(Point2d(pos-circle_m)) <= circle_r+pad) + if (L2Norm(Point2d(pos-m_origo)) <= m_radius+pad) { ret = true; } return ret; } + void Structure_Circle::getOrigo(Point2d & pos) + { + pos = m_origo; + } + + double Structure_Circle::getCharacteristicLength() + { + return 2*m_radius; + } + + } diff --git a/libs/NavierStokes_VPM_2D/src/Structure_Ellipse.cpp b/libs/NavierStokes_VPM_2D/src/Structure_Ellipse.cpp index 6cfed67..140caf6 100644 --- a/libs/NavierStokes_VPM_2D/src/Structure_Ellipse.cpp +++ b/libs/NavierStokes_VPM_2D/src/Structure_Ellipse.cpp @@ -24,6 +24,16 @@ namespace VPM { } + void Structure_Ellipse::getOrigo(Point2d & pos) + { + pos = m_origo; + } + + double Structure_Ellipse::getCharacteristicLength() + { + return 0.2;// not correct + } + bool Structure_Ellipse::isInside(const Point2d pos, const double pad) { diff --git a/libs/NavierStokes_VPM_2D/src/Structure_Flat.cpp b/libs/NavierStokes_VPM_2D/src/Structure_Flat.cpp index 765c92d..dc5b085 100644 --- a/libs/NavierStokes_VPM_2D/src/Structure_Flat.cpp +++ b/libs/NavierStokes_VPM_2D/src/Structure_Flat.cpp @@ -19,6 +19,16 @@ namespace VPM ret = true; } return ret; + + } + + void Structure_Flat::getOrigo(Point2d & pos) + { + } + + double Structure_Flat::getCharacteristicLength() + { + return 0.;// not correct } } diff --git a/libs/NavierStokes_VPM_2D/src/Structure_HalfCircle.cpp b/libs/NavierStokes_VPM_2D/src/Structure_HalfCircle.cpp index e1a42be..fcefe55 100644 --- a/libs/NavierStokes_VPM_2D/src/Structure_HalfCircle.cpp +++ b/libs/NavierStokes_VPM_2D/src/Structure_HalfCircle.cpp @@ -25,4 +25,13 @@ namespace VPM return ret; } + void Structure_HalfCircle::getOrigo(Point2d & pos) + { + } + + double Structure_HalfCircle::getCharacteristicLength() + { + return 0.;// not correct + } + } diff --git a/libs/NavierStokes_VPM_2D/src/Structure_Hilly.cpp b/libs/NavierStokes_VPM_2D/src/Structure_Hilly.cpp index 2ee8d34..80b39fa 100644 --- a/libs/NavierStokes_VPM_2D/src/Structure_Hilly.cpp +++ b/libs/NavierStokes_VPM_2D/src/Structure_Hilly.cpp @@ -22,4 +22,13 @@ namespace VPM return ret; } + void Structure_Hilly::getOrigo(Point2d & pos) + { + } + + double Structure_Hilly::getCharacteristicLength() + { + return 0.;// not correct + } + } diff --git a/libs/NavierStokes_VPM_2D/src/Structure_InverseCircle.cpp b/libs/NavierStokes_VPM_2D/src/Structure_InverseCircle.cpp index 8d0be55..86e33e3 100644 --- a/libs/NavierStokes_VPM_2D/src/Structure_InverseCircle.cpp +++ b/libs/NavierStokes_VPM_2D/src/Structure_InverseCircle.cpp @@ -36,4 +36,13 @@ namespace VPM return ret; } + void Structure_InverseCircle::getOrigo(Point2d & pos) + { + } + + double Structure_InverseCircle::getCharacteristicLength() + { + return 0.;// not correct + } + } diff --git a/libs/NavierStokes_VPM_2D/src/Structure_Wing.cpp b/libs/NavierStokes_VPM_2D/src/Structure_Wing.cpp new file mode 100644 index 0000000..f1e464f --- /dev/null +++ b/libs/NavierStokes_VPM_2D/src/Structure_Wing.cpp @@ -0,0 +1,271 @@ +#include "Structure_Wing.hpp" +#include + +namespace VPM +{ + + Structure_Wing::Structure_Wing() + { + m_origo = Point2d(.45, .0); + m_charlength = .12; + m_p.clear(); + m_p.push_back(Point2d( 1.00000000, 0.00000000+0.0128347)); + m_p.push_back(Point2d( 0.99939800, 0.00003500+0.0128347)); + m_p.push_back(Point2d( 0.99759200, 0.00013700+0.0128347)); + m_p.push_back(Point2d( 0.99458800, 0.00029600+0.0128347)); + m_p.push_back(Point2d( 0.99039300, 0.00049700+0.0128347)); + m_p.push_back(Point2d( 0.98501600, 0.00071900+0.0128347)); + m_p.push_back(Point2d( 0.97847000, 0.00093500+0.0128347)); + m_p.push_back(Point2d( 0.97077200, 0.00111200+0.0128347)); + m_p.push_back(Point2d( 0.96194000, 0.00121200+0.0128347)); + m_p.push_back(Point2d( 0.95199500, 0.00119700+0.0128347)); + m_p.push_back(Point2d( 0.94096100, 0.00103300+0.0128347)); + m_p.push_back(Point2d( 0.92886400, 0.00069400+0.0128347)); + m_p.push_back(Point2d( 0.91573500, 0.00015700+0.0128347)); + m_p.push_back(Point2d( 0.90160400, -0.00060000+0.0128347)); + m_p.push_back(Point2d( 0.88650500, -0.00159200+0.0128347)); + m_p.push_back(Point2d( 0.87047600, -0.00282900+0.0128347)); + m_p.push_back(Point2d( 0.85355300, -0.00431400+0.0128347)); + m_p.push_back(Point2d( 0.83577900, -0.00604800+0.0128347)); + m_p.push_back(Point2d( 0.81719700, -0.00802700+0.0128347)); + m_p.push_back(Point2d( 0.79785000, -0.01024400+0.0128347)); + m_p.push_back(Point2d( 0.77778500, -0.01269000+0.0128347)); + m_p.push_back(Point2d( 0.75705100, -0.01535700+0.0128347)); + m_p.push_back(Point2d( 0.73569800, -0.01823200+0.0128347)); + m_p.push_back(Point2d( 0.71377800, -0.02128900+0.0128347)); + m_p.push_back(Point2d( 0.69134200, -0.02449500+0.0128347)); + m_p.push_back(Point2d( 0.66844500, -0.02781400+0.0128347)); + m_p.push_back(Point2d( 0.64514200, -0.03120700+0.0128347)); + m_p.push_back(Point2d( 0.62149000, -0.03463100+0.0128347)); + m_p.push_back(Point2d( 0.59754500, -0.03804300+0.0128347)); + m_p.push_back(Point2d( 0.57336500, -0.04139700+0.0128347)); + m_p.push_back(Point2d( 0.54900900, -0.04464200+0.0128347)); + m_p.push_back(Point2d( 0.52453400, -0.04771900+0.0128347)); + m_p.push_back(Point2d( 0.50000000, -0.05056300+0.0128347)); + m_p.push_back(Point2d( 0.47546600, -0.05309900+0.0128347)); + m_p.push_back(Point2d( 0.45099100, -0.05525700+0.0128347)); + m_p.push_back(Point2d( 0.42663500, -0.05697900+0.0128347)); + m_p.push_back(Point2d( 0.40245500, -0.05822400+0.0128347)); + m_p.push_back(Point2d( 0.37851000, -0.05897400+0.0128347)); + m_p.push_back(Point2d( 0.35485800, -0.05923600+0.0128347)); + m_p.push_back(Point2d( 0.33155500, -0.05904600+0.0128347)); + m_p.push_back(Point2d( 0.30865800, -0.05845900+0.0128347)); + m_p.push_back(Point2d( 0.28622200, -0.05754700+0.0128347)); + m_p.push_back(Point2d( 0.26430200, -0.05637600+0.0128347)); + m_p.push_back(Point2d( 0.24294900, -0.05499400+0.0128347)); + m_p.push_back(Point2d( 0.22221500, -0.05342700+0.0128347)); + m_p.push_back(Point2d( 0.20215000, -0.05169400+0.0128347)); + m_p.push_back(Point2d( 0.18280300, -0.04980500+0.0128347)); + m_p.push_back(Point2d( 0.16422100, -0.04777300+0.0128347)); + m_p.push_back(Point2d( 0.14644700, -0.04561000+0.0128347)); + m_p.push_back(Point2d( 0.12952400, -0.04332600+0.0128347)); + m_p.push_back(Point2d( 0.11349500, -0.04092900+0.0128347)); + m_p.push_back(Point2d( 0.09839600, -0.03843100+0.0128347)); + m_p.push_back(Point2d( 0.08426500, -0.03584300+0.0128347)); + m_p.push_back(Point2d( 0.07113600, -0.03317000+0.0128347)); + m_p.push_back(Point2d( 0.05903900, -0.03041600+0.0128347)); + m_p.push_back(Point2d( 0.04800500, -0.02758600+0.0128347)); + m_p.push_back(Point2d( 0.03806000, -0.02468500+0.0128347)); + m_p.push_back(Point2d( 0.02922800, -0.02172200+0.0128347)); + m_p.push_back(Point2d( 0.02153000, -0.01870700+0.0128347)); + m_p.push_back(Point2d( 0.01498400, -0.01564900+0.0128347)); + m_p.push_back(Point2d( 0.00960700, -0.01255900+0.0128347)); + m_p.push_back(Point2d( 0.00541200, -0.00944300+0.0128347)); + m_p.push_back(Point2d( 0.00240800, -0.00630800+0.0128347)); + m_p.push_back(Point2d( 0.00060200, -0.00316000+0.0128347)); + m_p.push_back(Point2d( 0.00000000, 0.00000000+0.0128347)); + m_p.push_back(Point2d( 0.00060200, 0.00316500+0.0128347)); + m_p.push_back(Point2d( 0.00240800, 0.00630600+0.0128347)); + m_p.push_back(Point2d( 0.00541200, 0.00941600+0.0128347)); + m_p.push_back(Point2d( 0.00960700, 0.01248000+0.0128347)); + m_p.push_back(Point2d( 0.01498400, 0.01548900+0.0128347)); + m_p.push_back(Point2d( 0.02153000, 0.01844100+0.0128347)); + m_p.push_back(Point2d( 0.02922800, 0.02134800+0.0128347)); + m_p.push_back(Point2d( 0.03806000, 0.02421900+0.0128347)); + m_p.push_back(Point2d( 0.04800500, 0.02706200+0.0128347)); + m_p.push_back(Point2d( 0.05903900, 0.02987400+0.0128347)); + m_p.push_back(Point2d( 0.07113600, 0.03264400+0.0128347)); + m_p.push_back(Point2d( 0.08426500, 0.03536000+0.0128347)); + m_p.push_back(Point2d( 0.09839600, 0.03801100+0.0128347)); + m_p.push_back(Point2d( 0.11349500, 0.04058500+0.0128347)); + m_p.push_back(Point2d( 0.12952400, 0.04307100+0.0128347)); + m_p.push_back(Point2d( 0.14644700, 0.04545700+0.0128347)); + m_p.push_back(Point2d( 0.16422100, 0.04772900+0.0128347)); + m_p.push_back(Point2d( 0.18280300, 0.04987400+0.0128347)); + m_p.push_back(Point2d( 0.20215000, 0.05188500+0.0128347)); + m_p.push_back(Point2d( 0.22221500, 0.05375300+0.0128347)); + m_p.push_back(Point2d( 0.24294900, 0.05547000+0.0128347)); + m_p.push_back(Point2d( 0.26430200, 0.05702600+0.0128347)); + m_p.push_back(Point2d( 0.28622200, 0.05841400+0.0128347)); + m_p.push_back(Point2d( 0.30865800, 0.05962900+0.0128347)); + m_p.push_back(Point2d( 0.33155500, 0.06066000+0.0128347)); + m_p.push_back(Point2d( 0.35485800, 0.06149700+0.0128347)); + m_p.push_back(Point2d( 0.37851000, 0.06213300+0.0128347)); + m_p.push_back(Point2d( 0.40245500, 0.06256200+0.0128347)); + m_p.push_back(Point2d( 0.42663500, 0.06277900+0.0128347)); + m_p.push_back(Point2d( 0.45099100, 0.06277400+0.0128347)); + m_p.push_back(Point2d( 0.47546600, 0.06253000+0.0128347)); + m_p.push_back(Point2d( 0.50000000, 0.06202900+0.0128347)); + m_p.push_back(Point2d( 0.52453400, 0.06125400+0.0128347)); + m_p.push_back(Point2d( 0.54900900, 0.06019400+0.0128347)); + m_p.push_back(Point2d( 0.57336500, 0.05884500+0.0128347)); + m_p.push_back(Point2d( 0.59754500, 0.05721800+0.0128347)); + m_p.push_back(Point2d( 0.62149000, 0.05534400+0.0128347)); + m_p.push_back(Point2d( 0.64514200, 0.05325800+0.0128347)); + m_p.push_back(Point2d( 0.66844500, 0.05099300+0.0128347)); + m_p.push_back(Point2d( 0.69134200, 0.04857500+0.0128347)); + m_p.push_back(Point2d( 0.71377800, 0.04602900+0.0128347)); + m_p.push_back(Point2d( 0.73569800, 0.04337700+0.0128347)); + m_p.push_back(Point2d( 0.75705100, 0.04064100+0.0128347)); + m_p.push_back(Point2d( 0.77778500, 0.03784700+0.0128347)); + m_p.push_back(Point2d( 0.79785000, 0.03501700+0.0128347)); + m_p.push_back(Point2d( 0.81719700, 0.03217600+0.0128347)); + m_p.push_back(Point2d( 0.83577900, 0.02934700+0.0128347)); + m_p.push_back(Point2d( 0.85355300, 0.02655400+0.0128347)); + m_p.push_back(Point2d( 0.87047600, 0.02381700+0.0128347)); + m_p.push_back(Point2d( 0.88650500, 0.02115300+0.0128347)); + m_p.push_back(Point2d( 0.90160400, 0.01858000+0.0128347)); + m_p.push_back(Point2d( 0.91573500, 0.01611300+0.0128347)); + m_p.push_back(Point2d( 0.92886400, 0.01376900+0.0128347)); + m_p.push_back(Point2d( 0.94096100, 0.01156200+0.0128347)); + m_p.push_back(Point2d( 0.95199500, 0.00950800+0.0128347)); + m_p.push_back(Point2d( 0.96194000, 0.00762200+0.0128347)); + m_p.push_back(Point2d( 0.97077200, 0.00591500+0.0128347)); + m_p.push_back(Point2d( 0.97847000, 0.00440100+0.0128347)); + m_p.push_back(Point2d( 0.98501600, 0.00309200+0.0128347)); + m_p.push_back(Point2d( 0.99039300, 0.00200100+0.0128347)); + m_p.push_back(Point2d( 0.99458800, 0.00113700+0.0128347)); + m_p.push_back(Point2d( 0.99759200, 0.00051000+0.0128347)); + m_p.push_back(Point2d( 0.99939800, 0.00012800+0.0128347)); + m_p.push_back(Point2d( 1.00000000, 0.00000000+0.0128347)); + + } + Structure_Wing::~Structure_Wing() + { + } + + namespace { + //struct Segment { + // Point2d p1; + // Point2d p2; + //}; + + //bool onSegment(Segment l1, Point2d p) + //{ //check whether p is on the Segment or not + // bool ret = false; + // if ( p.x <= std::max(l1.p1.x, l1.p2.x) && + // p.x <= std::min(l1.p1.x, l1.p2.x) && + // (p.y <= std::max(l1.p1.y, l1.p2.y) && p.y <= std::min(l1.p1.y, l1.p2.y)) + // ) + // { + // ret = true; + // } + // return ret; + //} + + //int direction(Point2d a, Point2d b, Point2d c) + //{ + // int ret; + // int val = (b.y-a.y)*(c.x-b.x)-(b.x-a.x)*(c.y-b.y); + // if (val == 0) + // { + // ret = 0; //colinear + // } + // else if(val < 0) + // { + // ret = 2; //anti-clockwise direction + // } + // else + // { + // ret = 1; //clockwise direction + // } + // return ret; + //} + + //bool isIntersect(Segment l1, Segment l2) + //{ + // bool ret = false; + // //four direction for two Segments and points of other Segment + // int dir1 = direction(l1.p1, l1.p2, l2.p1); + // int dir2 = direction(l1.p1, l1.p2, l2.p2); + // int dir3 = direction(l2.p1, l2.p2, l1.p1); + // int dir4 = direction(l2.p1, l2.p2, l1.p2); + + // if(dir1 != dir2 && dir3 != dir4) + // { + // ret = true; //they are intersecting + // } + // else if(dir1==0 && onSegment(l1, l2.p1)) + // { //when p2 of Segment2 is on Segment1 + // ret = true; + // } + // else if(dir2==0 && onSegment(l1, l2.p2)) //when p1 of Segment2 is on Segment1 + // { + // ret = true; + // } + // else if(dir3==0 && onSegment(l2, l1.p1)) //when p2 of Segment1 is on Segment2 + // { + // ret = true; + // } + // else if(dir4==0 && onSegment(l2, l1.p2)) //when p1 of Segment1 is on Segment2 + // { + // ret = true; + // } + // return ret; + //} + bool ccw(Point2d A, Point2d B, Point2d C) + { + return (C.y-A.y) * (B.x-A.x) > (B.y-A.y) * (C.x-A.x); + } + bool intersect(Point2d A, Point2d B, Point2d C, Point2d D) + { + return (ccw(A,C,D) != ccw(B,C,D)) && (ccw(A,B,C) != ccw(A,B,D)); + } + } + + bool Structure_Wing::isInside(const Point2d pos, const double pad) + { + + bool ret = false; + int count_intersections = 0; + if (L2Norm(pos-Point2d(0,0))<1e-12) + { + ret = true; + } + else + { + Point2d p1 = pos; + Point2d p2 = Point2d(pos.x, +0.1); + //Point2d p1 = Point2d(-0.1, pos.y); + //Point2d p2 = Point2d(+1.1, pos.y); + for ( int i = 0; i < m_p.size()-1; i++) + { + Point2d s1 = m_p[i]; + Point2d s2 = m_p[i+1]; + if (intersect(p1, p2, s1, s2)) + //if (isIntersect(s1, s2)) + { + count_intersections +=1; + } + } + if (count_intersections%2 != 0) + { + ret = true; + } + } + //std::cerr<<"P=("< -#define C_VPM 0.125 +#define C_VPM 0.75 namespace VPM { VPM2d::VPM2d(int argc, char** argv) + : + m_structure_set(false) { m_advection = std::make_shared(); m_diffusion = std::make_shared(); @@ -56,9 +58,172 @@ namespace VPM return ret; } + // implements equation (32) in https://www-ljk.imag.fr/membres/Chloe.Mimeau/paper_IJNMF.pdf + Point2d VPM2d::getForces( + const ParticleField & pf, + const ParticleField & pf_old, + const double & delta_t + ) + { + Point2d F = Point2d(0,0); + Point2d origo; + m_structure->getOrigo(origo); + double Reynolds_inv = pf.params.m_nu/(m_structure->getCharacteristicLength()*L2Norm(pf.params.m_Uinfty)); + if (m_structure_set) + { + if (!pf.cartesianGrid) + { + std::cerr<<" ---> calculating forces called on IRregular grid\n"; + exit(0); + return Point2d(); + } +#define forcesversion2 +#ifdef forcesversion2 + int c = 3; + for (int j=c; jisInside(pf.regular_positions[i_j], 0.) && m_structure->isInside(pf.regular_positions[i_j], .4)) + //{ + F -= (pf.velocity[i_j]-pf_old.velocity[i_j])/delta_t*pf.params.m_vol; + //} + } + } + for (int i=c-1; i0) + { + j=pf.params.m_regular_num_py-c; + sign=-1; + } + int i_j = i +int(pf.params.m_regular_num_px)*j; + int i_jm = i +int(pf.params.m_regular_num_px)*(j-1); + int i_jp = i +int(pf.params.m_regular_num_px)*(j+1); + int ip_jp = i+1+int(pf.params.m_regular_num_px)*(j+1); + int ip_jm = i+1+int(pf.params.m_regular_num_px)*(j-1); + int im_jp = i-1+int(pf.params.m_regular_num_px)*(j+1); + int im_jm = i-1+int(pf.params.m_regular_num_px)*(j-1); + int im_j = i-1+int(pf.params.m_regular_num_px)*j; + int ip_j = i+1+int(pf.params.m_regular_num_px)*j; + double x = (pf.regular_positions[i_j].x-origo.x); + double y = (pf.regular_positions[i_j].y-origo.y); + double uy = (pf.velocity[i_jp].x - pf.velocity[i_jm].x)/(2*pf.params.m_dy); + double vx = (pf.velocity[ip_j].y - pf.velocity[im_j].y)/(2*pf.params.m_dx); + double vy = (pf.velocity[i_jp].y - pf.velocity[i_jm].y)/(2*pf.params.m_dy); + double uxx = (pf.velocity[ip_j].x - 2*pf.velocity[i_j].x + pf.velocity[im_j].x)/pf.params.m_vol; + double uyy = (pf.velocity[i_jp].x - 2*pf.velocity[i_j].x + pf.velocity[i_jm].x)/pf.params.m_vol; + double vxy = (pf.velocity[ip_jp].y - pf.velocity[ip_jm].y - pf.velocity[im_jp].y + pf.velocity[im_jm].y)/pf.params.m_vol; + double ut = (pf.velocity[i_j].x-pf_old.velocity[i_j].x)/delta_t; + F.x += sign*pf.params.m_dx*( + pf.velocity[i_j].x*pf.velocity[i_j].y + + pf.velocity[i_j].y*pf.omega[i_j]*y + - y*ut + + Reynolds_inv*(2*uxx + uyy + vxy)*y + - Reynolds_inv*(uy + vx)); + F.y += sign*pf.params.m_dx*( + 0.5*(pf.velocity[i_j].y*pf.velocity[i_j].y - pf.velocity[i_j].x*pf.velocity[i_j].x) + - pf.velocity[i_j].y*pf.omega[i_j]*x + + x*ut + - Reynolds_inv*(2*uxx + uyy + vxy)*x + - 2*Reynolds_inv*vy); + } + } + for (int j=c-1; j0) + { + i=pf.params.m_regular_num_px-c; + sign=+1; + } + int i_j = i +int(pf.params.m_regular_num_px)*j; + int i_jm = i +int(pf.params.m_regular_num_px)*(j-1); + int i_jp = i +int(pf.params.m_regular_num_px)*(j+1); + int ip_jp = i+1+int(pf.params.m_regular_num_px)*(j+1); + int ip_jm = i+1+int(pf.params.m_regular_num_px)*(j-1); + int im_jp = i-1+int(pf.params.m_regular_num_px)*(j+1); + int im_jm = i-1+int(pf.params.m_regular_num_px)*(j-1); + int im_j = i-1+int(pf.params.m_regular_num_px)*j; + int ip_j = i+1+int(pf.params.m_regular_num_px)*j; + double x = (pf.regular_positions[i_j].x-origo.x); + double y = (pf.regular_positions[i_j].y-origo.y); + double ux = (pf.velocity[ip_j].x - pf.velocity[im_j].x)/(2*pf.params.m_dx); + double uy = (pf.velocity[i_jp].x - pf.velocity[i_jm].x)/(2*pf.params.m_dy); + double vx = (pf.velocity[ip_j].y - pf.velocity[im_j].y)/(2*pf.params.m_dx); + double vxx = (pf.velocity[ip_j].y - 2*pf.velocity[i_j].y + pf.velocity[im_j].y)/pf.params.m_vol; + double vyy = (pf.velocity[i_jp].y - 2*pf.velocity[i_j].y + pf.velocity[i_jm].y)/pf.params.m_vol; + double uxy = (pf.velocity[ip_jp].x - pf.velocity[ip_jm].x - pf.velocity[im_jp].x + pf.velocity[im_jm].x)/pf.params.m_vol; + double vt = (pf.velocity[i_j].y-pf_old.velocity[i_j].y)/delta_t; + F.x += sign*pf.params.m_dx*( + 0.5*(pf.velocity[i_j].y*pf.velocity[i_j].y-pf.velocity[i_j].x*pf.velocity[i_j].x) + - pf.velocity[i_j].x*pf.omega[i_j]*y + - y*vt + + Reynolds_inv*(2*vyy + vxx + uxy)*y + + 2*Reynolds_inv*ux); + F.y += sign*pf.params.m_dx*( + - pf.velocity[i_j].x*pf.velocity[i_j].y + + pf.velocity[i_j].x*pf.omega[i_j]*x + + x*vt + - Reynolds_inv*(2*vyy + vxx + uxy)*x + + Reynolds_inv*(vx + uy)); + } + } +#else + for (int j=0; j structure) { m_advection->setStructure(structure); + m_structure = structure; + m_structure_set = true; + } + void VPM2d::getCharacteristicLength (double & charlength, bool & structureset) + { + if (m_structure_set) + { + charlength = m_structure->getCharacteristicLength(); + } + structureset = m_structure_set; } void VPM2d::run( @@ -82,6 +247,8 @@ namespace VPM m_filename_count = filename_count; m_save_init = save_init; + std::vector ft, fx, fy; + int stepcount = 0; double delta_t = 0.01; bool brexit = false; @@ -91,21 +258,61 @@ namespace VPM // ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| // std::cerr<<"------------- Step " < calculate Linfty norm of the velocity field.\n"; - m_advection->getLinftyGradVelocity(m_pf); - std::cerr<<" ---> max Linfty of velocity = "< calculate Linfty norm of the velocity field.\n"; + m_advection->getLinftyGradVelocity(m_pf); + std::cerr<<" ---> max Linfty of velocity = "<getCharacteristicLength(); + std::string t = std::to_string(m_filename_count); + std::ofstream myfile_ft; + std::ofstream myfile_fx; + std::ofstream myfile_fy; + std::ofstream myfile_dx; + std::ofstream myfile_dy; + myfile_ft.open(m_outputFile+"_t"+t+"_ft.dat", std::ofstream::binary); + myfile_fx.open(m_outputFile+"_t"+t+"_fx.dat", std::ofstream::binary); + myfile_fy.open(m_outputFile+"_t"+t+"_fy.dat", std::ofstream::binary); + myfile_dx.open(m_outputFile+"_t"+t+"_dx.dat", std::ofstream::binary); + myfile_dy.open(m_outputFile+"_t"+t+"_dy.dat", std::ofstream::binary); + double tmp; + for (unsigned int i=0; i(&ft[i]), sizeof(double)); + tmp = fx[i]; + myfile_fx.write(reinterpret_cast(&tmp), sizeof(double)); + tmp *= ForceToCoeff; + myfile_dx.write(reinterpret_cast(&tmp), sizeof(double)); + tmp = fy[i]; + myfile_fy.write(reinterpret_cast(&tmp), sizeof(double)); + tmp *= ForceToCoeff; + myfile_dy.write(reinterpret_cast(&tmp), sizeof(double)); +#else + myfile_ft.write(reinterpret_cast(&ft[i]), sizeof(double)); + myfile_fx.write(reinterpret_cast(&fx[i]), sizeof(double)); + myfile_fy.write(reinterpret_cast(&fy[i]), sizeof(double)); + myfile_dx.write(reinterpret_cast(&dx[i]), sizeof(double)); + myfile_dy.write(reinterpret_cast(&dy[i]), sizeof(double)); +#endif } + myfile_ft.close(); + myfile_fx.close(); + myfile_fy.close(); + myfile_dx.close(); + myfile_dy.close(); + m_filename_count++; + } - //if (stepcount==0) - //{ - // m_pf.Linfty_gradVelocity = C_VPM / delta_t; - //} + //if (stepcount==0) + //{ + // m_pf.Linfty_gradVelocity = C_VPM / delta_t; + //} - brexit = new_stepsize(final_time, m_pf, delta_t); + brexit = new_stepsize(final_time, m_pf, delta_t); + + m_pf_old = pf; // 2) Diffusion m_diffusion->Euler_step( m_pf, delta_t); @@ -120,9 +327,61 @@ namespace VPM } m_advection->calculateVelocity(m_pf); + Point2d F = getForces(m_pf, m_pf_old, delta_t); + + std::cerr<<"Forces="<getCharacteristicLength(); + writeParticlesToFile(m_outputFile, m_filename_count, m_pf); + + std::string t = std::to_string(m_filename_count); + std::ofstream myfile_ft; + std::ofstream myfile_fx; + std::ofstream myfile_fy; + std::ofstream myfile_dx; + std::ofstream myfile_dy; + myfile_ft.open(m_outputFile+"_t"+t+"_ft.dat", std::ofstream::binary); + myfile_fx.open(m_outputFile+"_t"+t+"_fx.dat", std::ofstream::binary); + myfile_fy.open(m_outputFile+"_t"+t+"_fy.dat", std::ofstream::binary); + myfile_dx.open(m_outputFile+"_t"+t+"_dx.dat", std::ofstream::binary); + myfile_dy.open(m_outputFile+"_t"+t+"_dy.dat", std::ofstream::binary); + double tmp; + for (unsigned int i=0; i(&ft[i]), sizeof(double)); + tmp = fx[i]; + myfile_fx.write(reinterpret_cast(&tmp), sizeof(double)); + tmp *= ForceToCoeff; + myfile_dx.write(reinterpret_cast(&tmp), sizeof(double)); + tmp = fy[i]; + myfile_fy.write(reinterpret_cast(&tmp), sizeof(double)); + tmp *= ForceToCoeff; + myfile_dy.write(reinterpret_cast(&tmp), sizeof(double)); +#else + myfile_ft.write(reinterpret_cast(&ft[i]), sizeof(double)); + myfile_fx.write(reinterpret_cast(&fx[i]), sizeof(double)); + myfile_fy.write(reinterpret_cast(&fy[i]), sizeof(double)); + myfile_dx.write(reinterpret_cast(&dx[i]), sizeof(double)); + myfile_dy.write(reinterpret_cast(&dy[i]), sizeof(double)); +#endif + } + myfile_ft.close(); + myfile_fx.close(); + myfile_fy.close(); + myfile_dx.close(); + myfile_dy.close(); + + } if (onestep) { brexit = true; }