From 926aecdbc06f1d5eb8c7efdbe1f1b7aed690cc6d Mon Sep 17 00:00:00 2001 From: Franz Georg Fuchs Date: Tue, 8 Dec 2020 13:54:00 +0100 Subject: [PATCH 1/4] calculate forces --- .gitignore | 3 + apps/NavierStokes_2D/src/App.cpp | 153 ++++++------ libs/NavierStokes_VPM_2D/include/Point2d.hpp | 4 +- .../NavierStokes_VPM_2D/include/Structure.hpp | 2 + .../include/Structure_Circle.hpp | 4 + .../include/Structure_Ellipse.hpp | 2 + .../include/Structure_Flat.hpp | 2 + .../include/Structure_HalfCircle.hpp | 2 + .../include/Structure_Hilly.hpp | 2 + .../include/Structure_InverseCircle.hpp | 2 + libs/NavierStokes_VPM_2D/include/VPM2d.hpp | 11 + .../src/Split_Diffusion.cpp | 2 + .../src/Structure_Circle.cpp | 17 +- .../src/Structure_Ellipse.cpp | 10 + .../src/Structure_Flat.cpp | 10 + .../src/Structure_HalfCircle.cpp | 9 + .../src/Structure_Hilly.cpp | 9 + .../src/Structure_InverseCircle.cpp | 9 + libs/NavierStokes_VPM_2D/src/VPM2d.cpp | 221 ++++++++++++++++-- 19 files changed, 375 insertions(+), 99 deletions(-) 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..5884883 100644 --- a/apps/NavierStokes_2D/src/App.cpp +++ b/apps/NavierStokes_2D/src/App.cpp @@ -1,6 +1,7 @@ #include "VPM2d.hpp" #include "Point2d.hpp" #include "InitialConditions.hpp" +#include "Split_Advection.hpp" #include "Parameters.hpp" #include "Structure_Flat.hpp" #include "Structure_Hilly.hpp" @@ -10,16 +11,20 @@ #include "Structure_Ellipse.hpp" #include +#include +#define OMPI_SKIP_MPICXX 1 +#include + int m_numberofParticlesx = 50; 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 ax = 1.0; @@ -35,11 +40,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; @@ -80,10 +81,6 @@ int main(int argc, char** argv) time = std::atof(argv[i+1]); eat = 2; } - if( arg == "--dt" && (i+1) < argc ) { - dt = std::atof(argv[i+1]); - eat = 2; - } if( arg == "--domain" && (i+4) < argc ) { domain_ll = VPM::Point2d(std::atof(argv[i+1]),std::atof(argv[i+2])); domain_ur = VPM::Point2d(std::atof(argv[i+3]),std::atof(argv[i+4])); @@ -105,26 +102,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,48 +137,6 @@ int main(int argc, char** argv) } } - VPM::VPM2d* vpm = new VPM::VPM2d(argc, argv); - - switch (structure_num) - { - case (0): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (1): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (2): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (3): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (4): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (5): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - } - std::shared_ptr remeshParams = std::make_shared( remesh_isOn, remesh_steps, @@ -215,13 +154,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 +170,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()); @@ -253,12 +204,54 @@ int main(int argc, char** argv) save_init = false; } + VPM::VPM2d* vpm = new VPM::VPM2d(argc, argv); + + switch (structure_num) + { + case (0): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (1): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (2): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (3): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (4): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (5): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + } - vpm->run(pf, omega, time, dt, params, - outputFile, fn_count, save_init - ); + std::cerr<<" here1 "<run(pf, time, outputFile, fn_count, save_init, false); + std::cerr<<" here2 "< 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/VPM2d.cpp b/libs/NavierStokes_VPM_2D/src/VPM2d.cpp index f1afd5c..695e981 100644 --- a/libs/NavierStokes_VPM_2D/src/VPM2d.cpp +++ b/libs/NavierStokes_VPM_2D/src/VPM2d.cpp @@ -21,6 +21,8 @@ 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,144 @@ namespace VPM return ret; } + 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 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 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; + F += sign*pf.params.m_dx*( + pf.velocity[i_j].x*pf.velocity[i_j].y + + pf.velocity[i_j].y*pf.omega[i_j]*x + - x*(pf.velocity[i_j]-pf_old.velocity[i_j])/delta_t + + Reynolds_inv*(2*uxx + uyy + vxy)*x + - Reynolds_inv*(uy + vx)); + } + } + 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 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 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; + F += 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*(pf.velocity[i_j]-pf_old.velocity[i_j])/delta_t + + Reynolds_inv*(2*vyy + vxx + uxy)*y + + 2*Reynolds_inv*ux); + } + } +#else + for (int j=0; j structure) { m_advection->setStructure(structure); + m_structure = structure; + m_structure_set = true; } void VPM2d::run( @@ -82,6 +219,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 +230,23 @@ 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 = "<Euler_step( m_pf, delta_t); @@ -120,9 +261,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; } From e3cca12de2ea0d91724e2f790c7c9fdde15af2cf Mon Sep 17 00:00:00 2001 From: Franz Georg Fuchs Date: Tue, 8 Dec 2020 13:56:49 +0100 Subject: [PATCH 2/4] calculate forces --- libs/NavierStokes_VPM_2D/src/VPM2d.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/libs/NavierStokes_VPM_2D/src/VPM2d.cpp b/libs/NavierStokes_VPM_2D/src/VPM2d.cpp index 695e981..b20de11 100644 --- a/libs/NavierStokes_VPM_2D/src/VPM2d.cpp +++ b/libs/NavierStokes_VPM_2D/src/VPM2d.cpp @@ -58,6 +58,7 @@ 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, From bd4e0287436e4775b449dc50ee5c6159513d12e5 Mon Sep 17 00:00:00 2001 From: Franz Georg Fuchs Date: Fri, 11 Dec 2020 15:18:22 +0100 Subject: [PATCH 3/4] corrected bugs --- apps/NavierStokes_2D/src/App.cpp | 2 -- libs/NavierStokes_VPM_2D/src/VPM2d.cpp | 45 ++++++++++++++++++-------- 2 files changed, 32 insertions(+), 15 deletions(-) diff --git a/apps/NavierStokes_2D/src/App.cpp b/apps/NavierStokes_2D/src/App.cpp index 5884883..ca9b632 100644 --- a/apps/NavierStokes_2D/src/App.cpp +++ b/apps/NavierStokes_2D/src/App.cpp @@ -247,9 +247,7 @@ int main(int argc, char** argv) } - std::cerr<<" here1 "<run(pf, time, outputFile, fn_count, save_init, false); - std::cerr<<" here2 "< -#define C_VPM 0.125 +#define C_VPM 0.75 namespace VPM { @@ -112,17 +112,26 @@ namespace VPM 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; - F += sign*pf.params.m_dx*( - pf.velocity[i_j].x*pf.velocity[i_j].y - + pf.velocity[i_j].y*pf.omega[i_j]*x - - x*(pf.velocity[i_j]-pf_old.velocity[i_j])/delta_t - + Reynolds_inv*(2*uxx + uyy + vxy)*x - - Reynolds_inv*(uy + vx)); + 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; j Date: Fri, 18 Dec 2020 10:21:34 +0100 Subject: [PATCH 4/4] add Wing structure, and pos. to set Re --- apps/NavierStokes_2D/src/App.cpp | 110 ++++--- .../include/Structure_Wing.hpp | 22 ++ libs/NavierStokes_VPM_2D/include/VPM2d.hpp | 2 + .../src/Structure_Wing.cpp | 271 ++++++++++++++++++ libs/NavierStokes_VPM_2D/src/VPM2d.cpp | 50 +++- 5 files changed, 410 insertions(+), 45 deletions(-) create mode 100644 libs/NavierStokes_VPM_2D/include/Structure_Wing.hpp create mode 100644 libs/NavierStokes_VPM_2D/src/Structure_Wing.cpp diff --git a/apps/NavierStokes_2D/src/App.cpp b/apps/NavierStokes_2D/src/App.cpp index ca9b632..25d53c3 100644 --- a/apps/NavierStokes_2D/src/App.cpp +++ b/apps/NavierStokes_2D/src/App.cpp @@ -6,10 +6,10 @@ #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 #include #define OMPI_SKIP_MPICXX 1 @@ -26,6 +26,7 @@ int main(int argc, char** argv) std::string outputFile; double time = 0.1; double nu = 0.; + double Re = 0.; double ax = 1.0; double ay = 1.0; @@ -77,6 +78,10 @@ int main(int argc, char** argv) nu = std::atof(argv[i+1]); eat = 2; } + if( arg == "--Re" && (i+1) < argc ) { + Re = std::atof(argv[i+1]); + eat = 2; + } if( arg == "--T" && (i+1) < argc ) { time = std::atof(argv[i+1]); eat = 2; @@ -137,6 +142,67 @@ int main(int argc, char** argv) } } + + VPM::VPM2d* vpm = new VPM::VPM2d(argc, argv); + + switch (structure_num) + { + case (0): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (1): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (2): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (3): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (4): + { + std::shared_ptr structure = std::make_shared(); + vpm->setStructure(structure); + } + break; + case (5): + { + std::shared_ptr structure = std::make_shared(); + 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( remesh_isOn, remesh_steps, @@ -204,48 +270,6 @@ int main(int argc, char** argv) save_init = false; } - VPM::VPM2d* vpm = new VPM::VPM2d(argc, argv); - - switch (structure_num) - { - case (0): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (1): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (2): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (3): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (4): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - case (5): - { - std::shared_ptr structure = std::make_shared(); - vpm->setStructure(structure); - } - break; - } - vpm->run(pf, time, outputFile, fn_count, save_init, false); 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 2a7b637..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: 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=("<getCharacteristicLength(); + } + structureset = m_structure_set; + } void VPM2d::run( ParticleField & pf, @@ -256,6 +264,44 @@ namespace VPM if ( stepcount==0&&m_save_init ) { writeParticlesToFile(m_outputFile, m_filename_count, m_pf); + double ForceToCoeff = 2./m_pf.params.m_Uinfty.x/m_pf.params.m_Uinfty.x/m_structure->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++; } @@ -292,7 +338,7 @@ namespace VPM m_filename_count++; stepcount++; - if(stepcount%100==0) + if(stepcount%10==0) { double ForceToCoeff = 2./m_pf.params.m_Uinfty.x/m_pf.params.m_Uinfty.x/m_structure->getCharacteristicLength(); writeParticlesToFile(m_outputFile, m_filename_count, m_pf); @@ -311,7 +357,7 @@ namespace VPM double tmp; for (unsigned int i=0; i(&ft[i]), sizeof(double)); tmp = fx[i]; myfile_fx.write(reinterpret_cast(&tmp), sizeof(double));