LeesEdwards3D

apps/showcases/LeesEdwards/CMakeLists.txt

@@ -8,6 +8,7 @@ waLBerla_generate_target_from_python(NAME LeesEdwardsGen
         LeesEdwards_Stream.cpp LeesEdwards_Stream.h
         LeesEdwards_Setter.cpp LeesEdwards_Setter.h
         LeesEdwards_PackInfo.cpp LeesEdwards_PackInfo.h
+        LeesEdwards_PackInfo_Velocity.cpp LeesEdwards_PackInfo_Velocity.h
         LeesEdwards_InfoHeader.h)
 
 waLBerla_add_executable(NAME LeesEdwards

apps/showcases/LeesEdwards/LeesEdwards.cpp

@@ -34,6 +34,7 @@
 using namespace walberla;
 
 using PackInfo_T  = lbm::LeesEdwards_PackInfo;
+using PackInfoVelocity_T  = pystencils::LeesEdwards_PackInfo_Velocity;
 using flag_t      = walberla::uint8_t;
 using FlagField_T = FlagField< flag_t >;
 
@@ -122,7 +123,7 @@ int main(int argc, char** argv)
    auto parameters = walberlaEnv.config()->getOneBlock("Parameters");
 
    const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
-   const real_t offset    = parameters.getParameter< real_t >("offset", real_c(10));
+   // const real_t offset    = parameters.getParameter< real_t >("offset", real_c(10));
 
    const double remainingTimeLoggerFrequency =
       parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
@@ -133,23 +134,39 @@ int main(int argc, char** argv)
    BlockDataID forceFieldID   = field::addToStorage< VectorField_T >(blocks, "force", real_t(0), field::fzyx);
    BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(1.0), field::fzyx);
 
-   pystencils::LeesEdwards_Setter setterSweep(densityFieldID, forceFieldID, pdfFieldID);
+   pystencils::LeesEdwards_Setter setterSweep(densityFieldID, forceFieldID, pdfFieldID, velFieldID);
    for (auto& block : *blocks)
       setterSweep(&block);
 
-   // create time loop
-   SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
-
    // create communication for PdfField
    blockforest::communication::UniformBufferedScheme< Stencil_T > communication(blocks);
    communication.addPackInfo(make_shared< PackInfo_T >(pdfFieldID));
+   communication.addPackInfo(make_shared< PackInfoVelocity_T >(velFieldID));
 
    pystencils::LeesEdwards_Collision CollisionSweep(densityFieldID, forceFieldID, pdfFieldID, velFieldID);
    pystencils::LeesEdwards_Stream StreamSweep(densityFieldID, forceFieldID, pdfFieldID, velFieldID);
 
+   auto normalTimeStep = [&]() {
+      communication.communicate();
+      for (auto& block : *blocks)
+      {
+         StreamSweep(&block);
+         CollisionSweep(&block);
+      }
+
+   };
+
+   // create time loop
+   SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
+
+   std::function< void() > timeStep;
+   timeStep = std::function< void() >(normalTimeStep);
+   timeloop.add() << BeforeFunction(timeStep) << Sweep([](IBlock*) {}, "time step");
+
+
    // add LBM sweep and communication to time loop
-   timeloop.add() << Sweep(CollisionSweep, "collision");
-   timeloop.add() << Sweep(StreamSweep, "stream") << AfterFunction(communication, "communication");
+//   timeloop.add() << Sweep(CollisionSweep, "collision");
+//   timeloop.add() << Sweep(StreamSweep, "stream") << AfterFunction(communication, "communication");
 
    // log remaining time
    timeloop.addFuncAfterTimeStep(timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),

apps/showcases/LeesEdwards/LeesEdwards.prm

 
 Parameters 
 {
-	timesteps       10001;
-	vtkWriteFrequency 500;
+	timesteps       101;
+	vtkWriteFrequency 100;
 	remainingTimeLoggerFrequency 3; // in seconds
 	offset 10;
 }
@@ -10,7 +10,7 @@ Parameters
 DomainSetup
 {
    blocks        <  1,    1, 1 >;
-   cellsPerBlock <  64, 64, 1 >;
-   periodic      <  1,    1, 0 >;
+   cellsPerBlock <  64, 64, 64 >;
+   periodic      <  1,    1, 1 >;
 }
 

apps/showcases/LeesEdwards/LeesEdwards.py

@@ -8,7 +8,7 @@ from lbmpy.advanced_streaming.utility import get_accessor, Timestep
 from lbmpy.creationfunctions import create_lb_update_rule
 from lbmpy.updatekernels import create_stream_pull_with_output_kernel
 from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
-from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_info_header, generate_pack_info_for_field
 from lbmpy_walberla import RefinementScaling, generate_boundary, generate_lb_pack_info
 
 import sympy as sp
@@ -17,7 +17,7 @@ import sympy as sp
 # For advanced streaming patterns (AA, EsoTwist) the timestep is seperated into Odd and Even steps. In each of these
 # steps a different streaming is executed. For more common two population methods this is not the case.
 # In lbmpy we indicate this with Timestep.BOTH
-streaming_pattern = "push"
+streaming_pattern = "pull"
 accessor = get_accessor(streaming_pattern, Timestep.BOTH)
 
 omega = 1.0  # relaxation rate of first component
@@ -25,7 +25,7 @@ shear_velocity = 0.5  # shear velocity
 shear_dir = 0  # direction of shear flow
 shear_dir_normal = 1  # direction normal to shear plane, for interpolation
 
-stencil = LBStencil(Stencil.D2Q9)
+stencil = LBStencil(Stencil.D3Q19)
 q = stencil.Q
 dim = stencil.D
 
@@ -42,7 +42,7 @@ u_p = sp.Piecewise((1, sp.And(type_all_numbers(counters[1] <= 0, 'int'), points_
                    (0, True)) * shear_velocity
 
 lb_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=omega, compressible=True,
-                      velocity_input=velocity_field.center_vector + sp.Matrix([u_p, 0]),
+                      velocity_input=velocity_field.center_vector + sp.Matrix([u_p, 0, 0]),
                       density_input=density_field, force_model=ForceModel.GUO,
                       force=force_field.center_vector, kernel_type='collide_only')
 lbm_opt = LBMOptimisation(symbolic_field=pdfs)
@@ -56,10 +56,10 @@ for s in to_insert:
     collision = collision.new_with_inserted_subexpression(s)
 ma = []
 for a, c in zip(collision.main_assignments, collision.method.stencil):
-    if c[shear_dir_normal] == -1:
-        b = (True, False)
-    elif c[shear_dir_normal] == 1:
+    if c[2] == -1:
         b = (False, True)
+    elif c[2] == 1:
+        b = (True, False)
     else:
         b = (False, False)
     a = Assignment(a.lhs, a.rhs.replace(points_down, b[0]))
@@ -71,7 +71,7 @@ stream = create_stream_pull_with_output_kernel(collision.method, pdfs, pdfs_tmp,
                                                {'density': density_field, 'velocity': velocity_field},
                                                accessor=accessor)
 
-init = macroscopic_values_setter(collision.method, velocity=(0, 0),
+init = macroscopic_values_setter(collision.method, velocity=velocity_field.center_vector,
                                  pdfs=pdfs, density=density_field.center,
                                  streaming_pattern=streaming_pattern)
 
@@ -88,6 +88,7 @@ with CodeGeneration() as ctx:
 
     # communication
     generate_lb_pack_info(ctx, 'LeesEdwards_PackInfo', stencil, pdfs, streaming_pattern=streaming_pattern)
+    generate_pack_info_for_field(ctx, 'LeesEdwards_PackInfo_Velocity', velocity_field)
 
     # Info header containing correct template definitions for stencil and field
     generate_info_header(ctx, 'LeesEdwards_InfoHeader',