Coupling force in NS and AD lattice

[hhan]

Dear Commuity,

I created my own case based on thermalPorousPlate3d, which is a case of NS coupled AD lattice concentration diffusion. Initially, I used the coupling method Struct NavierStokesAdvectionDiffusionVelocityCoupling from navierStokesAdvectionDiffusionCoupling.h for coupling, as shown in the following code:

SuperLatticeCoupling coupling(
NavierStokesAdvectionDiffusionVelocityCoupling{},
names::NavierStokes{}, sLatticeNS,
names::Concentration0{}, sLatticeAD);

This code produced the desired results. However, when I used either the coupling force terms I created and added to navierStokesAdvectionDiffusionCoupling.h or the built-in NavierStokesAdvectionDiffusionCoupling method, the program failed to produce results. The generated images were blank, and the output showed the following warnings:

“./tmp/imageData/data/_EuklidNormphysVelocityiT0001200.p” line 21: warning: matrix contains missing or undefined values
^C”./tmp/imageData/data/_densityiT0001200.p” line 21: warning: matrix contains missing or undefined values
Warning: empty cb range [0:0], adjusting to [-1:1]
Warning: empty cb range [0:0], adjusting to [-1:1]

How can I solve this problem? Below are my complete code and the struct code I created.

The struct I created:

struct DensityDrivenBuoyancyCoupling {
static constexpr OperatorScope scope = OperatorScope::PerCellWithParameters;

struct GRAVITY : public descriptors::FIELD_BASE<0,1> { };
struct DENSITY_DIFFERENCE : public descriptors::FIELD_BASE<1> { };

using parameters = meta::list<GRAVITY, DENSITY_DIFFERENCE>;

template <typename CELLS, typename PARAMETERS>
void apply(CELLS& cells, PARAMETERS& parameters) any_platform
{
using V = typename CELLS::template value_t<names::NavierStokes>::value_t;
using DESCRIPTOR = typename CELLS::template value_t<names::NavierStokes>::descriptor_t;

auto& cellNSE = cells.template get<names::NavierStokes>();
auto& cellADE = cells.template get<names::Concentration0>();

auto force = cellNSE.template getFieldPointer<descriptors::FORCE>();
auto densityDifference = parameters.template get<DENSITY_DIFFERENCE>();
auto gravity = parameters.template get<GRAVITY>();
V rho = cellADE.computeRho(); // get concentration
V buoyancyForceFactor = densityDifference * rho;

for (unsigned iD = 0; iD < DESCRIPTOR::d; ++iD) {
force[iD] += buoyancyForceFactor * gravity[iD];
}
// compute force
V u[DESCRIPTOR::d] { };
cellNSE.computeU(u);
cellADE.template setField<descriptors::VELOCITY>(u);
}
};

.cpp file code:

#include “olb3D.h”
#include “olb3D.hh”
using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;

using T = FLOATING_POINT_TYPE;
typedef D3Q19<FORCE> NSDESCRIPTOR;
typedef D3Q7<VELOCITY> ADDESCRIPTOR;

const int iTperiod = 100; // amount of timesteps when new boundary conditions are reset and results are visualized
const T diffusion = 1.e-4; // diffusion coefficient for advection-diffusion equation
const T radius = 1.5e-04; // particles radius
const T partRho = 0.7; // particles density
const T maxPhysT = 10.; // max. simulation time in s, SI unit

const int N = 19;
const T latticeU = 0.01;
const T physLength = 1;
const T physVelocity = 0.01;
const T physViscosity = 1.5e-5;
const T physDensity = 1.26;
T inletRadius = 0.01;
const T Re = 50; // Reynolds number

void prepareGeometry( UnitConverter<T,NSDESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry )
{

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

clout << “Prepare Geometry …” << std::endl;
T nx = 6*inletRadius;
T deltaX = converter.getPhysDeltaX();

std::vector<T> origin { -nx/4, -nx/4, -nx/2 };
std::vector<T> extend { nx/4, nx/4, deltaX };

IndicatorCuboid3D<T> inlet(extend, origin);

superGeometry.rename( 0, 2 );
superGeometry.rename( 2, 1, {1,1,1});
superGeometry.rename( 2, 3, inlet);
superGeometry.clean();

superGeometry.checkForErrors();
superGeometry.print();

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

void prepareLattice( SuperLattice<T, NSDESCRIPTOR>& sLatticeNS,
SuperLattice<T, ADDESCRIPTOR>& sLatticeAD,
UnitConverter<T,NSDESCRIPTOR> const& converter,
SuperGeometry<T,3>& superGeometry,
T diffusion)
{

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

const T omega = converter.getLatticeRelaxationFrequency();
const T omegaAD = converter.getLatticeRelaxationFrequencyFromDiffusivity<ADDESCRIPTOR>( diffusion );
// Material=1 –> bulk dynamics
// Material=3 –> bulk dynamics (inflow)
auto inflowIndicator = superGeometry.getMaterialIndicator({1,3});
sLatticeNS.defineDynamics<ForcedBGKdynamics>(inflowIndicator);
sLatticeAD.defineDynamics<AdvectionDiffusionBGKdynamics>(inflowIndicator);

// Material=2 –> bounce-back / zero distribution dynamics
sLatticeNS.defineDynamics<ForcedBGKdynamics>(superGeometry, 2);
sLatticeAD.defineDynamics<ZeroDistributionDynamics>(superGeometry, 2);

// Setting of the boundary conditions
setInterpolatedVelocityBoundary<T,NSDESCRIPTOR>(sLatticeNS, omega, superGeometry, 3);
setZeroDistributionBoundary<T,ADDESCRIPTOR>(sLatticeAD, superGeometry, 2);
setAdvectionDiffusionTemperatureBoundary<T,ADDESCRIPTOR>(sLatticeAD, superGeometry, 3);

// Initial conditions
std::vector<T> zero(3,T());
AnalyticalConst3D<T,T> rho1( 1. );
AnalyticalConst3D<T,T> rho0( 1.e-8 );
AnalyticalConst3D<T,T> u1(0., 0., converter.getCharLatticeVelocity());
AnalyticalConst3D<T,T> u0(0., 0., 0.);
AnalyticalConst3D<T,T> u2(converter.getCharLatticeVelocity(), 0., 0.);
AnalyticalConst3D<T,T> force(zero);

// Initialize all values of distribution functions to their local equilibrium
// 1 bulk
sLatticeNS.defineRhoU( superGeometry, 1, rho1, u2 );
sLatticeNS.iniEquilibrium( superGeometry, 1, rho1, u2 );
sLatticeNS.defineField<descriptors::FORCE>(superGeometry, 1, force);
// 2 boundary
sLatticeNS.defineRhoU( superGeometry, 2, rho1, u0 );
sLatticeNS.iniEquilibrium( superGeometry, 2, rho1, u0 );
sLatticeNS.defineField<descriptors::FORCE>(superGeometry, 2, force);
// 3 inlet
sLatticeNS.defineRhoU( superGeometry.getMaterialIndicator({3}), rho1, u0 );
sLatticeNS.iniEquilibrium( superGeometry.getMaterialIndicator({3}), rho1, u0 );
sLatticeNS.defineField<descriptors::FORCE>(superGeometry, 3, force);

// 1 bulk
sLatticeAD.defineRho( superGeometry, 1, rho0 );
sLatticeAD.iniEquilibrium( superGeometry.getMaterialIndicator({1}), rho0, u2 );
// 2 boundary
sLatticeAD.defineRho( superGeometry, 2, rho0 );
sLatticeAD.iniEquilibrium( superGeometry.getMaterialIndicator({2}), rho0, u0 );
// 3 inlet
sLatticeAD.defineRho( superGeometry, 3, rho1 );
sLatticeAD.iniEquilibrium( superGeometry.getMaterialIndicator({3}), rho1, u0 );

sLatticeNS.setParameter<descriptors::OMEGA>(omega);
sLatticeAD.setParameter<descriptors::OMEGA>(omegaAD);

// Lattice initialize
sLatticeNS.initialize();
sLatticeAD.initialize();

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

void setBoundaryValues( SuperLattice<T, NSDESCRIPTOR>& sLatticeNS,
UnitConverter<T,NSDESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry )
{

OstreamManager clout( std::cout, “setBoundaryValues” );
}

void getResults( SuperLattice<T, NSDESCRIPTOR>& sLatticeNS,
SuperLattice<T, ADDESCRIPTOR>& sLatticeAD,
UnitConverter<T,NSDESCRIPTOR> const& converter, int iT,
SuperGeometry<T,3>& superGeometry,
Timer<double>& timer )
{

OstreamManager clout( std::cout, “getResults” );
SuperVTMwriter3D<T> vtmWriter( “bifurcation3d_fluid” );
SuperVTMwriter3D<T> vtmWriterAD( “bifurcation3d_particle” );
SuperLatticePhysVelocity3D<T, NSDESCRIPTOR> velocity( sLatticeNS, converter );
SuperLatticeVelocity3D<T, NSDESCRIPTOR> latticeVelocity( sLatticeNS);

SuperLatticePhysPressure3D<T, NSDESCRIPTOR> pressure( sLatticeNS, converter );
SuperLatticeDensity3D<T, ADDESCRIPTOR> particles( sLatticeAD );
SuperLatticePhysField3D<T, ADDESCRIPTOR, descriptors::VELOCITY> extField(
sLatticeAD, converter.getConversionFactorVelocity());

vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( pressure );
vtmWriterAD.addFunctor( particles );
vtmWriterAD.addFunctor( extField );

if ( iT == 0 ) {
SuperLatticeGeometry3D<T, NSDESCRIPTOR> geometry( sLatticeNS, superGeometry );
SuperLatticeCuboid3D<T, NSDESCRIPTOR> cuboid( sLatticeNS );
SuperLatticeRank3D<T, NSDESCRIPTOR> rank( sLatticeNS );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
vtmWriterAD.createMasterFile();

}
if ( iT % iTperiod == 0) {
vtmWriter.write( iT );
vtmWriterAD.write( iT );

// GIF Writer
SuperEuklidNorm3D<T> normVel( velocity );
HyperplaneLattice3D<T> gifLattice(
superGeometry.getCuboidGeometry(),
Hyperplane3D<T>()
.centeredIn(superGeometry.getCuboidGeometry().getMotherCuboid())
.normalTo({0, -1, 0}),
600);
BlockReduction3D2D<T> planeReductionVelocity( normVel, gifLattice, BlockDataSyncMode::ReduceOnly );
BlockReduction3D2D<T> planeReductionParticles( particles, gifLattice, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
BlockGifWriter<T> gifWriter;
gifWriter.write(planeReductionVelocity, iT, “vel” );
gifWriter.write(planeReductionParticles, iT, “density”);
heatmap::write(planeReductionVelocity, iT);
heatmap::write(planeReductionParticles, iT);
}
}

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

// === 1st Step: Initialization ===

olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout, “main” );
UnitConverterFromResolutionAndRelaxationTime<T,NSDESCRIPTOR> const converter(
int {N}, // resolution: number of voxels per charPhysL
(T) 0.557646, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) inletRadius*2., // charPhysLength: reference length of simulation geometry
(T) Re*1.5e-5/( inletRadius*2 ), // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 1.5e-5, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 1.225 // 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(“bifurcation3d”);

// === 2nd Step: Prepare Geometry ===
// Instantiation of an empty cuboidGeometry
T nx = 6*inletRadius;
T deltaX = converter.getPhysDeltaX();
std::vector<T> origin { -nx/2, -nx/2, -nx/2 };
std::vector<T> extend { nx, nx, nx };
IndicatorCuboid3D<T> cuboid(extend, origin);

const int noOfCuboids = singleton::mpi().getSize();
CuboidGeometry3D<T> cuboidGeometry(cuboid, deltaX, noOfCuboids);
cuboidGeometry.setPeriodicity(true, false, false);

// Instantiation of an empty loadBalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );
// Default instantiation of superGeometry
SuperGeometry<T,3> superGeometry( cuboidGeometry, loadBalancer);

prepareGeometry( converter, superGeometry );

// === 3rd Step: Prepare Lattice ===
SuperLattice<T, NSDESCRIPTOR> sLatticeNS( superGeometry );
SuperLattice<T, ADDESCRIPTOR> sLatticeAD( superGeometry );

//prepareLattice and setBoundaryConditions
prepareLattice( sLatticeNS, sLatticeAD, converter, superGeometry, diffusion );

SuperLatticeCoupling coupling(
DensityDrivenBuoyancyCoupling{},
names::NavierStokes{}, sLatticeNS,
names::Concentration0{}, sLatticeAD);
coupling.setParameter<DensityDrivenBuoyancyCoupling::GRAVITY>(
Vector<T,3>{0.,0.,9.8} );
coupling.setParameter<DensityDrivenBuoyancyCoupling::DENSITY_DIFFERENCE>(
converter.getPhysDensity() – partRho );
/*
SuperLatticeCoupling coupling(
NavierStokesAdvectionDiffusionCoupling{},
names::NavierStokes{}, sLatticeNS,
names::Temperature{}, sLatticeAD);
coupling.setParameter<NavierStokesAdvectionDiffusionCoupling::T0>(
0.);
coupling.setParameter<NavierStokesAdvectionDiffusionCoupling::FORCE_PREFACTOR>(
10 * Vector<T,3>{0.0,0.0,1.0});

SuperLatticeCoupling coupling(
NavierStokesAdvectionDiffusionVelocityCoupling{},
names::NavierStokes{}, sLatticeNS,
names::Concentration0{}, sLatticeAD);
*/
// === 4th Step: Main Loop with Timer ===
Timer<double> timer( converter.getLatticeTime( maxPhysT ),
superGeometry.getStatistics().getNvoxel() );
timer.start();

for (std::size_t iT = 0; iT <= converter.getLatticeTime(maxPhysT); ++iT) {
sLatticeAD.setParameter<descriptors::LATTICE_TIME>(iT);
getResults( sLatticeNS, sLatticeAD, converter, iT, superGeometry, timer );
setBoundaryValues( sLatticeNS, converter, iT, superGeometry );

sLatticeNS.collideAndStream();
coupling.execute();
sLatticeAD.collideAndStream();
}
timer.stop();
timer.printSummary();
}

Thank you for reply!

[Adrian]

At first glance, your new DensityDrivenBuoyancyCoupling coupling operator looks ok.

Are you sure that the problem is in the coupling? Your error / warning could just as well be caused by a problem in the image output only. How does the VTK look? It could also be the case that everything code-wise works but the simulation diverges (i.e. a model error occurs).

[hhan]

Dear Adrian,

Thank you for reply!

I cheak my vtk output and also can’t get results from vtk file. I don’t know what the problem is.

When i don’t use the force coupling and only couple the velocity between AD and NS lattice, I can get seemly right output gif. Can model errors still occur in this case?

I wonder if it is caused by too little coupled force. Hope you can suggest some solutions.

Hope your reply!

[mathias]

hard to tell what is going wrong, that would be a case for the advanced section at our next spring school or in a common project..

[Author]

Posts