Skip to content

Incomplete indicator shape

Viewing 8 posts - 1 through 8 (of 8 total)
  • Author
  • #6219

    Dear everyone,
    Recently I have tried to simulate flow in finned(sin-shaped fin) channel,but I have difficulty in creating the sin-shaped fin.The post-process image both in openlb and paraview seem to be weird.Here is my code:

    template <typename S>
    IndicatorSin2D<S>::IndicatorSin2D(Vector<S,2> _center, S _radius, S _theta, S _w, S _A)
    : center(_center),
    this->_myMin = _center – radius;
    this->_myMax = _center + radius;

    template <typename S>
    bool IndicatorSin2D<S>::operator()(bool output[], const S input[])
    output[0] = (input[0] <= this->_myMax[0] && input[0] >= this->_myMin[0]) &&
    ( input[1] <= this->A * sin(this->w * (input[0] – this->theta))
    && input[1] >= 0);
    return output[0];

    #include “olb2D.h”
    #ifndef OLB_PRECOMPILED // Unless precompiled version is used,
    #include “olb2D.hh” // include full template code
    #include <vector>
    #include <cmath>
    #include <iostream>
    #include <fstream>

    using namespace olb;
    using namespace olb::descriptors;
    using namespace olb::graphics;
    using namespace olb::util;
    using namespace std;

    typedef double T;


    //#define M_PI 3.14159265358979323846;

    // Parameters for the simulation setup
    const T lx = 0.005; // length of the channel
    const T ly = 0.001; // height of the channel
    const int N = 60; // resolution
    const T maxPhysT = 0.5; // max. simulation time in s, SI unit
    const T epsilon = 1.e-4; // precision of the convergence (residuum)

    const T Thot = 289; // temperature of the lower wall in Kelvin
    const T Tcold = 288; // temperature of the fluid and top wall in Kelvin
    //parameters of sinfin
    const T amplitude = 0.0002;
    const T w = 5e3 * M_PI;
    const T theta = 0.0024;

    /// Stores geometry information in form of material numbers
    void prepareGeometry(SuperGeometry2D<T>& superGeometry,
    ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter)

    OstreamManager clout(std::cout,”prepareGeometry”);
    clout << “Prepare Geometry …” << std::endl;

    superGeometry.rename(2,1,1,1 );
    Vector<T, 2>center(0.5*lx, 0.5*amplitude);
    IndicatorSin2D<T>sinfin(center, 0.5*amplitude, theta, w, amplitude);

    std::vector<T> extend( 2, T(0) );
    extend[0] = lx;
    extend[1] = converter.getPhysLength(1);
    std::vector<T> origin( 2, T(0) );
    IndicatorCuboid2D<T> bottom(extend, origin);

    origin[1] = ly-converter.getPhysLength(1);
    IndicatorCuboid2D<T> top(extend, origin);

    std::vector<T> extendwallleft(2,T(0));
    extendwallleft[0] = converter.getPhysLength(1);
    extendwallleft[1] = ly;
    std::vector<T> originwallleft(2,T(0));
    originwallleft[0] = 0.0;
    originwallleft[1] = 0.0;
    IndicatorCuboid2D<T> inflow(extendwallleft, originwallleft);

    std::vector<T> extendwallright(2,T(0));
    extendwallright[0] = converter.getPhysLength(1);
    extendwallright[1] = ly;
    std::vector<T> originwallright(2,T(0));
    originwallright[0] = lx-converter.getPhysLength(1);
    originwallright[1] = 0.0;
    IndicatorCuboid2D<T> outflow(extendwallright, originwallright);

    /// Set material numbers for bottom, top and pertubation
    // superGeometry.rename(1,5,test);
    superGeometry.rename(2,3,1,inflow );
    superGeometry.rename(2,4,1,outflow );

    /// Removes all not needed boundary voxels outside the surface
    /// Removes all not needed boundary voxels inside the surface


    clout << “Prepare Geometry … OK” << std::endl;

    void prepareLattice( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
    SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
    SuperLattice2D<T, TDESCRIPTOR>& ADlattice,
    Dynamics<T, NSDESCRIPTOR> &bulkDynamics,
    Dynamics<T, TDESCRIPTOR>& advectionDiffusionBulkDynamics,
    Dynamics<T, TDESCRIPTOR>& sourcedAdvectionDiffusionBGKdynamics,
    SuperGeometry2D<T>& superGeometry )

    OstreamManager clout(std::cout,”prepareLattice”);

    T Tomega = converter.getLatticeThermalRelaxationFrequency();
    const T omega = converter.getLatticeRelaxationFrequency();
    /// define lattice Dynamics
    clout << “defining dynamics” << endl;
    ADlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, TDESCRIPTOR>());
    NSlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, NSDESCRIPTOR>());
    ADlattice.defineDynamics(superGeometry, 1, &advectionDiffusionBulkDynamics);
    ADlattice.defineDynamics(superGeometry, 2, &advectionDiffusionBulkDynamics);
    ADlattice.defineDynamics(superGeometry, 3, &advectionDiffusionBulkDynamics);
    ADlattice.defineDynamics(superGeometry, 4, &advectionDiffusionBulkDynamics);
    ADlattice.defineDynamics(superGeometry, 5, &advectionDiffusionBulkDynamics);
    ADlattice.defineDynamics(superGeometry, 6, &advectionDiffusionBulkDynamics);
    //ADlattice.defineDynamics(superGeometry, 5, &sourcedAdvectionDiffusionBGKdynamics);

    NSlattice.defineDynamics(superGeometry, 1, &bulkDynamics);
    NSlattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, NSDESCRIPTOR>());
    NSlattice.defineDynamics(superGeometry, 3, &bulkDynamics);//xinjia
    NSlattice.defineDynamics(superGeometry, 4, &bulkDynamics);
    NSlattice.defineDynamics(superGeometry, 5, &instances::getBounceBack<T, NSDESCRIPTOR>());
    NSlattice.defineDynamics(superGeometry, 6, &instances::getBounceBack<T, NSDESCRIPTOR>());
    //NSlattice.defineDynamics(superGeometry, 7, &instances::getNoDynamics<T, NSDESCRIPTOR>());
    /// sets boundary
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 3);
    //setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 4);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 5);
    //setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 5);
    setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 6);
    //setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 7);
    setLocalPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 4);
    //setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({2, 3}));
    setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);

    // Vector<T, 2>center(0.5*lx, 0.5*amplitude);
    // IndicatorSin2D<T>sinfin(center, 0.5*amplitude, theta, w, amplitude);
    // setBouzidiZeroVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 5, sinfin);
    // setBouzidiVelocityBoundary<T,TDESCRIPTOR>(ADlattice, superGeometry, 5, sinfin);

    AnalyticalConst2D<T,T> rho(1.225);
    AnalyticalConst2D<T,T> u0(0.0, 0.0);
    AnalyticalConst2D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
    AnalyticalConst2D<T,T> T_hot(converter.getLatticeTemperature(Thot));

    T maxVelocity = converter.getCharLatticeVelocity();
    T distance2Wall = converter.getConversionFactorLength()/2;
    Poiseuille2D<T> u( superGeometry, 3, maxVelocity, distance2Wall );

    /// for each material set Rho, U and the Equilibrium
    NSlattice.defineRhoU(superGeometry, 1, rho, u0);
    NSlattice.iniEquilibrium(superGeometry, 1, rho, u0);
    NSlattice.defineRhoU(superGeometry, 2, rho, u0);
    NSlattice.iniEquilibrium(superGeometry, 2, rho, u0);
    NSlattice.defineRhoU(superGeometry, 3, rho, u);
    NSlattice.iniEquilibrium(superGeometry, 3, rho, u);
    NSlattice.defineRhoU( superGeometry, 4, rho, u0 );
    NSlattice.iniEquilibrium( superGeometry, 4, rho, u0 );
    NSlattice.defineRhoU(superGeometry, 5, rho, u0);
    NSlattice.iniEquilibrium(superGeometry, 5, rho, u0);
    NSlattice.defineRhoU(superGeometry, 6, rho, u0);
    NSlattice.iniEquilibrium(superGeometry, 6, rho, u0);
    // NSlattice.defineRhoU(superGeometry, 7, rho, u0);
    // NSlattice.iniEquilibrium(superGeometry, 7, rho, u0);

    ADlattice.defineRho(superGeometry, 1, T_cold);
    ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
    ADlattice.defineRho(superGeometry, 2, T_cold);
    ADlattice.iniEquilibrium(superGeometry, 2, T_cold, u0);
    ADlattice.defineRho(superGeometry, 3, T_cold);
    ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u);
    ADlattice.defineRho(superGeometry, 4, T_cold);
    ADlattice.iniEquilibrium(superGeometry, 4, T_cold, u0);
    ADlattice.defineRho(superGeometry, 5, T_hot);
    ADlattice.iniEquilibrium(superGeometry, 5, T_hot, u0);
    ADlattice.defineRho(superGeometry, 6, T_hot);
    ADlattice.iniEquilibrium(superGeometry, 6, T_hot, u0);
    // ADlattice.defineRho(superGeometry, 7, T_hot);
    // ADlattice.iniEquilibrium(superGeometry, 7, T_hot, u0);

    clout << “Prepare Lattice … OK” << std::endl;

    void setBoundaryValues(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
    SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
    SuperLattice2D<T, TDESCRIPTOR>& ADlattice,
    int iT, SuperGeometry2D<T>& superGeometry)


    void getResults(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
    SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
    SuperLattice2D<T, TDESCRIPTOR>& ADlattice, int iT,
    SuperGeometry2D<T>& superGeometry,
    Timer<T>& timer,
    bool converged)

    OstreamManager clout(std::cout,”getResults”);

    SuperVTMwriter2D<T> vtkWriter(“rectangle2″);
    SuperLatticePhysVelocity2D<T, NSDESCRIPTOR> velocity(NSlattice, converter);
    SuperLatticePhysPressure2D<T, NSDESCRIPTOR> presure(NSlattice, converter);
    SuperLatticePhysTemperature2D<T, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);
    SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry1(NSlattice, superGeometry);
    vtkWriter.addFunctor( presure );
    vtkWriter.addFunctor( velocity );
    vtkWriter.addFunctor( temperature );
    vtkWriter.addFunctor( geometry1 );

    const int saveIter = converter.getLatticeTime(0.005);

    if (iT == 0) {
    SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
    SuperLatticeCuboid2D<T, NSDESCRIPTOR> cuboid(NSlattice);
    SuperLatticeRank2D<T, NSDESCRIPTOR> rank(NSlattice);
    /// Writes the VTK files and prints statistics
    if (iT%saveIter == 0 || converged) {
    /// Timer console output

    /// Lattice statistics console output
    BlockReduction2D2D<T> planeReduction( temperature, 600, BlockDataSyncMode::ReduceOnly );
    // // write output as JPEG
    heatmap::plotParam<T> jpeg_Param;
    jpeg_Param.maxValue = Thot;
    jpeg_Param.minValue = Tcold;
    heatmap::write(planeReduction, iT,jpeg_Param);

    SuperEuklidNorm2D<T, NSDESCRIPTOR> normVel( velocity );
    BlockReduction2D2D<T> planeReduction1( normVel, 600, BlockDataSyncMode::ReduceOnly );
    // write output as JPEG
    heatmap::plotParam<T> jpeg_Param2;
    jpeg_Param2.maxValue = converter.getCharPhysVelocity()*1.5;
    heatmap::write(planeReduction1, iT, jpeg_Param2);


    int main(int argc, char *argv[])

    /// === 1st Step: Initialization ===
    OstreamManager clout(std::cout,”main”);
    olbInit(&argc, &argv);

    ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> converter(
    (T) ly/N, // physDeltaX
    (T) 4.63e-7,
    (T) ly, // charPhysLength
    (T) 0.005, // charPhysVelocity
    (T) 1.47e-5, // physViscosity
    (T) 1.225, // physDensity
    (T) 0.025105, // physThermalConductivity
    (T) 1006.3, // physSpecificHeatCapacity
    (T) 0.00425,//Ra * 1e-5 * 1e-5 / Pr / 9.81 / (Thot – Tcold) / pow(ly, 3), // physThermalExpansionCoefficient
    (T) Tcold, // charPhysLowTemperature
    (T) Thot // charPhysHighTemperature


    /// === 2nd Step: Prepare Geometry ===
    std::vector<T> extend(2,T());
    extend[0] = lx;
    extend[1] = ly;
    std::vector<T> origin(2,T());
    IndicatorCuboid2D<T> cuboid(extend, origin);

    /// Instantiation of a cuboidGeometry with weights
    const int noOfCuboids = singleton::mpi().getSize();
    const int noOfCuboids = 1;//7–>1
    CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), noOfCuboids);

    //cuboidGeometry.setPeriodicity(true, false);

    HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);

    SuperGeometry2D<T> superGeometry(cuboidGeometry, loadBalancer, 2);

    prepareGeometry(superGeometry, converter);

    /// === 3rd Step: Prepare Lattice ===

    SuperLattice2D<T, TDESCRIPTOR> ADlattice(superGeometry);
    //maybe we can specific the geometry domain,eg inlet
    SuperLattice2D<T, NSDESCRIPTOR> NSlattice(superGeometry);

    ForcedBGKdynamics<T, NSDESCRIPTOR> NSbulkDynamics(

    AdvectionDiffusionBGKdynamics<T, TDESCRIPTOR> TbulkDynamics (
    SourcedAdvectionDiffusionBGKdynamics<T, TDESCRIPTOR> SinDynamics (

    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
    // This coupling must be necessarily be put on the Navier-Stokes lattice!!
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//

    std::vector<T> dir{0.0, 1.0};//1–>5

    T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
    converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();

    NavierStokesAdvectionDiffusionCouplingGenerator2D<T,NSDESCRIPTOR> coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(ly), boussinesqForcePrefactor, converter.getLatticeTemperature(Tcold), 1., dir);

    NSlattice.addLatticeCoupling(coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 1, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 2, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 3, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 4, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 5, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 6, coupling, ADlattice);
    NSlattice, ADlattice,
    NSbulkDynamics, TbulkDynamics,SinDynamics,
    superGeometry );
    //xinjia sboundarycondition

    /// === 4th Step: Main Loop with Timer ===
    Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() );

    util::ValueTracer<T> converge(converter.getLatticeTime(0.005),epsilon);//50–>10
    for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT); ++iT) {

    if (converge.hasConverged()) {
    clout << “Simulation converged.” << endl;
    getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());

    clout << “Time ” << iT << “.” << std::endl;


    /// === 5th Step: Definition of Initial and Boundary Conditions ===
    //setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry);

    /// === 6th Step: Collide and Stream Execution ===


    /// === 7th Step: Computation and Output of the Results ===
    getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());


    I have two questions.The first is how to make the color of fin the same as its temperature.The other one is why the sin-fin shape is incomplete.I would appreciate it if anyone could help me solve this!

    Yours sincerely,


    Dear Mike,

    I recomment using Paraview for more involved psotprocessing. What Do you mean by “incomplete”?



    Hello Mathias,
    Thanks for your apply! By “incomplete” I mean the fin shape showed in the post-process image is not what I imagine, it seems part of the fin is missing due to resolution setting or something else.
    The image is listed blow,I intend to design a sin-shaped fin, but it came out as a triangle like.
    temperature image

    • This reply was modified 2 years, 3 months ago by Mike.

    Does the shape move towards what you expect if you increase the resolution (Set as deltaX in the unit converter)? Note that the problem could also be the discretization used by the heatmap writer.

    How does the raw VTK output look?


    Dear Adrian,

    Thanks for your reply! Yes, the fin-shape get finer if I increase the resolution, which does not make sense and prolong the simulation time.Also, I have tried to find if I can solve this problem by change some parameters in heatmap namespace,but I can not.What’s more, I notice that the fin-shape in paraview is similar to what it looks like in openlb getresults,which persuade me it may has nothing to do with heatmap, I don’t know if I get this right……

    As for VTK output, I know little about it.And it is too long and unreadable that I hesitate to put it here.But if it helps, I would be glad to do this!


    • This reply was modified 2 years, 3 months ago by Mike.

    I just tried your code and increasing the resolution seems to work correctly. Due to the use of regular grids in LBM, increasing the resolution is the primary way of resolving geometry details. The corresponding increase in (sequential) simulation time is handled by the very good parallelizability. As your example currently only uses ~20k cells you can significantly increase the resolution even on a single compute node.

    w.r.t. VTK: There is no need to look at the raw files for this issue, you already checked the output in ParaView (which visualizes the VTK files)


    Dear Adrian,

    Thank you again for your reply! Earlier I was thinking something goes wrong within my code about IndicatorSin2D, or there exist a parameter I don’t find deep inside the Indicator like a geometry resolution or something, which can be adjusted to refine the Indicator shape.

    Also I’m curious about what lead to the problem that the fin-shape is not the same as I suppose.Is it because that I created in a wrong way or is it about the display problem or something else.What’s more, when I use ug or freecad to create 3D geometry model and import it to openlb, the fin-shape is pretty nice.What is the mechanism behind this.I would appreciate it if you can help!



    As you implemented a custom indicator for your fin shape (which is the correct way to go for non-primitive shapes that are not easily reproduced using indicator arithmetic) there is likely something wrong in the custom operator(). Can you tell us more about the specific “fin” geometry you want to model?

    The reason why this is not an issue for CAD-imported geometries is that the STL reader indicator is designed to correctly import arbitrary STLs.

Viewing 8 posts - 1 through 8 (of 8 total)
  • You must be logged in to reply to this topic.