[Anand]

Dear Adrian,
Thank you for your reply.
I noticed this mistake.
Now, everything is working properly on a cluster.

Regards
Ananda

[Anand]

Hi Adrian,
I am running the simulations on clusters.
Our system has prepared a module (installed the OpenLB).
After that, I unzip the software on my network drive and made changes in config.mk to run the simulation in parallel mode (activated mpi).
Then, I created a folder in the “examples/turbulent/Pipe”, where i put the my .cpp, .stl, Makefile, module.mk and .sh file.

My .sh file is as follows:

#!/bin/bash
#SBATCH –qos bbdefault
#SBATCH –ntasks 72
#SBATCH –nodes 1
#SBATCH –time 72:00:00
#SBATCH –mail-type ALL
#SBATCH –account=accountname
module purge; module load bluebear
module load OpenLB/1.4-0-foss-2021a
module load gnuplot/5.4.1-GCCcore-10.2.0
export OMP_NUM_THREADS=1
make clean
make cleanbuild
make
mpiexec -np 72 ./Pipe

I just submit it and it works.

I have noticed some issues with the cluster for example:
If I use the following line in the code, my simulation diverges but the same line works fine on my computer.
“CirclePowerLawTurbulent3D<T> uF( superGeometry,3,maxVelocity[0], 7, 0.05, T(0));”

So I just replace this line with the following lines and it worked very well with Clusters.

std::vector<T> origin = { pipelength, piperadius, piperadius};
std::vector<T> axis = { 1, 0, 0 };
CirclePowerLawTurbulent3D<T> POu(origin, axis, piperadius, maxVelocity[0], 7, 0.05, T(1));

Please let me know what you think.
Thank you
Regards
Ananda

[Anand]

Hi Adrian,
I use following line and it is working (Pipe.cpp is a code)…

“mpiexec -np 72 ./Pipe”

[Anand]

Issue resolved…….

[Anand]

hi,
what is the meaning of “mpiexec.hydra”?
thank you
regards
Anand

[Anand]

Dear Fedor,
This is really useful information.
Thank you
Regards
Ananda

[Anand]

Dear Fedor,
Thank you very much for the reply.
I checked the bifurcation example.
Can you please confirm about rho1 or rho0?
How does it relate to the volume fraction, I mean mathematical relation?
What is the volume fraction in the bifurcation example?

Thank you
Regards
Ananda

[Anand]

Hi Adrian,
One more thing I need to mention.
If I keep everything the same and just change the descriptor to either D3Q19 or D3Q27, the results are different.
At the outlet
D3Q27 gives me the same velocity profile as inlet (Poiseuille)
D3Q19 give me fluctuation at the outlet no particular pattern

[Anand]

1) I calculated using: Re = (inletvelocity*pipediameter)/(kinematicviscosity); (inletvelocity=0.9423m/s, Pipediameter=0.04m; kinematicviscosity=0.0044 (kg/m-s)/1168 (kg/m3))
I also see the same Reynolds number in the UnitConversion Information Parameters output.

2) I use N=50 and M=50

3) Timespan = 10 second

4) I use .pvd file in paraview and put a vertical line at the outlet and plot velocity vs diameter of pipe

5) I just modify your the existing examples.

[Anand]

I check the velocity profile at the end of the pipe.
I see the profile is exactly the same at inlet (Poiseuille flow profile) even though I Reynolds number is 10005 and I use SmagorinskyBGKdynamics.

Please can someone explain to me what is wrong with my code, am I missing something?

Thank you in advance.
Regards
Anand

[Anand]

I managed to resolve this problem by changing the resolutions.

[Anand]

Hi Adrian,
Thank you for quick reply….
Here I have copy past the code, if you wish, please have alook….

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

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;
using namespace std;
typedef double T;
const T maxPhysT = 10.; // max. simulation time in s, SI unit
int N = 60; // resolution of the model
const T density = 1168; // Liquid Density (kg/m3)
const T viscosity = 0.0044/density; // Liquid kinematic viscosity (m2/s)
const T inletvelocity = 0.774; // Inlet Velocity (m/s)
const T pipediameter = 0.0405; // Pipe Diameter (m)
const T piperadius = pipediameter/2; // Pipe Radius (m)
const T Re = (inletvelocity*pipediameter)/(viscosity); // defined as 1/kinematic viscosity
const T pipelength = 1.5*(1.359*pipediameter*std::pow(Re, 0.25)); // Pipe Length (m) calculated using turbulence entrace length formula plus 1.5 time higher

// Stores geometry information in form of material numbers
void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry3D<T>& superGeometry )
{

OstreamManager clout(std::cout, “prepareGeometry”);

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

Vector<T, 3> inletcenter(-converter.getConversionFactorLength()* 0.2, piperadius, piperadius);
Vector<T, 3> outletcenter(pipelength, piperadius, piperadius);
IndicatorCylinder3D<T> pipe(inletcenter, outletcenter, piperadius);

superGeometry.rename(0, 2);

superGeometry.rename(2, 1, pipe);

Vector<T, 3> origin(0, piperadius, piperadius);
Vector<T, 3> extend = origin;

// Set material number for inflow
origin[0] = -converter.getConversionFactorLength() * 2;
extend[0] = converter.getConversionFactorLength() * 2;
IndicatorCylinder3D<T> inflow(origin, extend, piperadius);
superGeometry.rename(2, 3, 1, inflow);

// Set material number for outflow
origin[0] = pipelength – 2 * converter.getConversionFactorLength();
extend[0] = pipelength + 2 * converter.getConversionFactorLength();
IndicatorCylinder3D<T> outflow(extend, origin, piperadius);
superGeometry.rename(2, 4, 1, outflow);

// Removes all not needed boundary voxels outside the surface
superGeometry.clean();
// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();

superGeometry.print();

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

// Set up the geometry of the simulation
void prepareLattice(SuperLattice3D<T, DESCRIPTOR>& sLattice,
UnitConverter<T, DESCRIPTOR>const& converter,
Dynamics<T, DESCRIPTOR>& bulkDynamics,
SuperGeometry3D<T>& superGeometry)
{

OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

const T omega = converter.getLatticeRelaxationFrequency();

// Material=0 –>do nothing
sLattice.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );

// Material=1 –>bulk dynamics
sLattice.defineDynamics( superGeometry, 1, &bulkDynamics );

// material=2 –> depends on boundary type

sLattice.defineDynamics( superGeometry, 2, &instances::getBounceBack<T, DESCRIPTOR>() );

// Material=3 –>bulk dynamics
sLattice.defineDynamics( superGeometry, 3, &bulkDynamics );
setInterpolatedVelocityBoundary<T,DESCRIPTOR>(sLattice, omega, superGeometry, 3);

// Material=4 –>bulk dynamics
sLattice.defineDynamics( superGeometry, 4, &bulkDynamics );
setInterpolatedPressureBoundary<T,DESCRIPTOR>(sLattice, omega, superGeometry, 4);
// Initial conditions
AnalyticalConst3D<T,T> ux( 0. );
AnalyticalConst3D<T,T> uy( 0. );
AnalyticalConst3D<T,T> uz( 0. );
AnalyticalConst3D<T,T> rho( 1. );
AnalyticalComposed3D<T,T> u( ux,uy,uz );

//Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU( superGeometry.getMaterialIndicator({1, 3, 4}), rho, u );
sLattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 3, 4}), rho, u );

// Make the lattice ready for simulation
sLattice.initialize();

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

}

int main( int argc, char* argv[] )
{
/// === 1st Step: Initialization ===
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);
OstreamManager clout( std::cout,”main” );
clout.setMultiOutput(true);

UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR> converter(
int {N}, // resolution: number of voxels per charPhysL
(T) 0.5018, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) pipediameter, // charPhysLength: reference length of simulation geometry
(T) inletvelocity, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) viscosity, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) density // physDensity: physical density in __kg / m^3__
);

// Prints the converter log as console output
converter.print();
// Writes the converter log in a file
converter.write(“Pipe”);

// === 2nd Step: Prepare Geometry ===

Vector<T, 3> inletcenter(0, piperadius, piperadius);
Vector<T, 3> outletcenter(pipelength, piperadius, piperadius);
IndicatorCylinder3D<T> pipe(inletcenter, outletcenter, piperadius);
IndicatorLayer3D<T> extendedDomain(pipe, converter.getConversionFactorLength());

// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = std::min( 16*N,2*singleton::mpi().getSize() );
#else
const int noOfCuboids = 10;
#endif // ifdef PARALLEL_MODE_MPI
CuboidGeometry3D<T> cuboidGeometry(extendedDomain, converter.getConversionFactorLength(), noOfCuboids);

// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);

// Instantiation of a superGeometry
SuperGeometry3D<T> superGeometry(cuboidGeometry, loadBalancer, 2);

prepareGeometry(converter, superGeometry);

// === 3rd Step: Prepare Lattice ===
SuperLattice3D<T, DESCRIPTOR> sLattice( superGeometry );
const T omega = converter.getLatticeRelaxationFrequency();

SmagorinskyBGKdynamics<T, DESCRIPTOR>bulkDynamics(omega, instances::getBulkMomenta<T, DESCRIPTOR>(),0.1);
}

[Author]

Posts