Fixed CoVo

apps/showcases/COVO/COVO.cpp

@@ -49,8 +49,8 @@ class CovoInit
 {
  public:
 
-   CovoInit( real_t circulation, real_t radius, Vector3<real_t> center, real_t lattice_velocity ) :
-      circulation_( circulation ), radius_(radius), center_(center), lattice_velocity_(lattice_velocity){}
+   CovoInit( real_t circulation, real_t radius, Vector3<real_t> center, real_t lattice_velocity, real_t ss_lattice ) :
+      circulation_( circulation ), radius_(radius), center_(center), lattice_velocity_(lattice_velocity), ss_lattice_(ss_lattice){}
 
    void operator()(const shared_ptr< StructuredBlockForest >& SbF, BlockDataID densityFieldID, BlockDataID velFieldID ) const;
 
@@ -60,6 +60,7 @@ class CovoInit
    const real_t radius_;
    const Vector3<real_t> center_;
    const real_t lattice_velocity_;
+   const real_t ss_lattice_;
 }; // class ShearProfile
 
 void CovoInit::operator()(const shared_ptr< StructuredBlockForest >& SbF, BlockDataID densityFieldID, BlockDataID velFieldID ) const
@@ -74,7 +75,6 @@ void CovoInit::operator()(const shared_ptr< StructuredBlockForest >& SbF, BlockD
 
          real_t beta = 2 * (radius_ * radius_);
          real_t circ2 = circulation_ * circulation_;
-         real_t ss_lattice = 1 / sqrt(3.0);
          real_t delta_x = (globalCell[0] - center_[0]) * (globalCell[0] - center_[0]);
          real_t delta_y = (globalCell[1] - center_[1]) * (globalCell[1] - center_[1]);
          real_t var = exp((- delta_x - delta_y) / beta);
@@ -82,12 +82,12 @@ void CovoInit::operator()(const shared_ptr< StructuredBlockForest >& SbF, BlockD
          real_t d_psi_x = circulation_ * (2 * radius_ * (globalCell[1] - center_[1]) * var) / beta;
          real_t d_psi_y = circulation_ * (2 * radius_ * (globalCell[0] - center_[0]) * var) / beta;
 
-         real_t rho_loc = exp( -(circ2 / 2 * ss_lattice) * var * beta);
+         real_t rho_loc = exp( -(circ2 / 2 * ss_lattice_) * var * beta);
 
          velField->get(x, y, z, 0) = d_psi_x + lattice_velocity_;
          velField->get(x, y, z, 1) = -d_psi_y;
 
-         densityField->get(x, y, z) = rho_loc;
+         densityField->get(x, y, z) = 1.0 * rho_loc;  // lattice density is 1.0
       )
       // clang-format on
    }
@@ -106,13 +106,38 @@ int main(int argc, char** argv)
    // read parameters
    auto parameters = walberlaEnv.config()->getOneBlock("Parameters");
 
-   // const real_t omega     = parameters.getParameter< real_t >("omega", real_c(1.4));
-   // const real_t u_max     = parameters.getParameter< real_t >("u_max", real_t(0.05));
-   const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
-   const real_t circulation = parameters.getParameter< real_t >("circulation", 1.0);
-   const real_t radius = parameters.getParameter< real_t >("radius", 1.0);
-   const real_t lattice_velocity = parameters.getParameter< real_t >("lattice_velocity", 1.0);
-   const Vector3<real_t> center = parameters.getParameter<Vector3<real_t>>("center");
+   const real_t ReferenceLength = parameters.getParameter< real_t >("ReferenceLength", real_t(0.1));
+   const real_t ConvectionVelocity = parameters.getParameter< real_t >("ConvectionVelocity", real_t(170.0));
+   const real_t Radius = parameters.getParameter< real_t >("Radius", real_t(0.005));
+   // const real_t Density = parameters.getParameter< real_t >("Density", real_t(1.16));
+   const real_t Circulation = parameters.getParameter< real_t >("Circulation", real_t(34.728));
+   const real_t S_speed = parameters.getParameter< real_t >("S_speed", real_t(347.1));
+   const real_t Time = parameters.getParameter< real_t >("Time", real_t(0.01));
+
+   const real_t ReynoldsNumber = parameters.getParameter< real_t >("ReynoldsNumber", real_t(1082000.0));
+   // const real_t LatticeDensity = parameters.getParameter< real_t >("LatticeDensity", real_t(1));
+   const real_t LatticeVelocity = parameters.getParameter< real_t >("LatticeVelocity", real_t(0.05));
+   const real_t LatticeLength = parameters.getParameter< real_t >("LatticeLength", real_t(400));
+
+
+   // real_t  Cd = Density / LatticeDensity;
+   real_t Cl = ReferenceLength / LatticeLength;
+   real_t Cu = ConvectionVelocity / LatticeVelocity;
+   real_t Ct = Cl / Cu;
+
+   real_t LatticeCirculation = Circulation / Cl * Ct;
+   real_t LatticeRadius = Radius / Cl;
+
+   real_t SSLattice = S_speed / Cu;
+
+   real_t KinematicVicosity = (LatticeLength * LatticeVelocity) / ReynoldsNumber;
+   real_t omega = 2.0 / (6.0 * KinematicVicosity + 1);
+
+   uint_t timesteps = uint_t(Time / Ct);
+   WALBERLA_LOG_INFO_ON_ROOT("Simulating with " << timesteps << " timesteps");
+
+   Vector3<real_t> center(LatticeLength / 2, LatticeLength / 2, 0);
+
 
    const double remainingTimeLoggerFrequency =
       parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
@@ -120,10 +145,10 @@ int main(int argc, char** argv)
    // create fields
    BlockDataID pdfFieldID     = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
    BlockDataID velFieldID     = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0), field::fzyx);
-   BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.16), field::fzyx);
+   BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1), field::fzyx);
 
 
-   CovoInit InitData(circulation, radius, center, lattice_velocity);
+   CovoInit InitData(LatticeCirculation, LatticeRadius, center, LatticeVelocity, SSLattice);
    InitData(blocks, densityFieldID, velFieldID);
 
    pystencils::COVO_MacroSetter setterSweep(densityFieldID, pdfFieldID, velFieldID);
@@ -137,7 +162,7 @@ int main(int argc, char** argv)
    blockforest::communication::UniformBufferedScheme< Stencil_T > communication(blocks);
    communication.addPackInfo(make_shared< PackInfo_T >(pdfFieldID));
 
-   pystencils::COVO_Sweep UpdateSweep(densityFieldID, pdfFieldID, velFieldID);
+   pystencils::COVO_Sweep UpdateSweep(densityFieldID, pdfFieldID, velFieldID, omega);
 
    // add LBM sweep and communication to time loop
    timeloop.add() << BeforeFunction(communication, "communication") << Sweep(UpdateSweep, "LB stream & collide");

apps/showcases/COVO/COVO.prm

@@ -8,7 +8,7 @@ Parameters
     Density 1.16;               // Kg / m^3
     Circulation 34.728;         // m^2 / s
     S_speed 347.1;              // m / s
-    Time 0.01;                  // s
+    Time 0.001;                  // s
 
     ReynoldsNumber 1082000.0;
 
@@ -23,6 +23,6 @@ Parameters
 DomainSetup
 {
    blocks        <  1,    1, 1 >;
-   cellsPerBlock <  LatticeLength, LatticeLength, 1 >;
+   cellsPerBlock <  400, 400, 1 >;
    periodic      <  1,    1, 0 >;
 }
\ No newline at end of file