Reply To: WARNING: no discreteNormal is found

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