Forces

.gitignore

@@ -40,3 +40,6 @@
 *.aux
 CMakeInit.txt
 Session.vim
+*.dat
+*.png
+build/

apps/NavierStokes_2D/src/App.cpp

@@ -1,26 +1,32 @@
 #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 <string>
+#define OMPI_SKIP_MPICXX 1
+#include <mpi.h>
 
 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 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<VPM::Structure_Wing> structure = std::make_shared<VPM::Structure_Wing>();
+                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="<<nu<<std::endl;
+        }
     }
 
     std::shared_ptr<VPM::RemeshParams> remeshParams = std::make_shared<VPM::RemeshParams>(
@@ -215,13 +220,12 @@ int main(int argc, char** argv)
             bc_to_yr
             );
 
-    std::shared_ptr<VPM::Parameters> params = std::make_shared<VPM::Parameters>(
+    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<VPM::Point2d> positions;
     std::vector<double> 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<VPM::Point2d> positions;
+        std::vector<double> 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;
 
 }

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_ )
             :

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:
 

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;
 
     };
 

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:
 

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:
 

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:
 

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:
 

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:
 

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<Point2d> m_p;
+
+    };
+
+}

libs/NavierStokes_VPM_2D/include/VPM2d.hpp

@@ -34,6 +34,8 @@ namespace VPM
             const bool onestep = false
             );
         void setStructure(std::shared_ptr<VPM::Structure> 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<Split_Diffusion> m_diffusion;
         std::shared_ptr<Split_Advection> m_advection;
         std::shared_ptr<Split_Source> m_source;
 
+        bool m_structure_set;
+        std::shared_ptr<VPM::Structure> m_structure;
+
         ParticleField m_pf;
+        ParticleField m_pf_old;
 
         std::string m_outputFile;
         int m_filename_count;

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;
     }
 

libs/NavierStokes_VPM_2D/src/Structure_Circle.cpp

@@ -5,22 +5,33 @@ namespace VPM
 
     Structure_Circle::Structure_Circle()
     {
+       m_origo = Point2d(-.5, .111);
+       m_radius = .2;
     }
     Structure_Circle::~Structure_Circle()
     {
     }
 
     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;
+    }
+
+
 }