Re: Zero force

[kolotinsky]

Dear Mathias,
Of course, velocity must go to infinity, I missed up it yesterday. In my code I used forcedMRT dynamics.
Here is my code

//////////////////////////////////////////////////////////////////////////
#include “olb3D.h”
#include “olb3D.hh” // use only generic version!
#include <cstdlib>
#include <iostream>

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

typedef double T;
#define DESCRIPTOR ForcedMRTD3Q19Descriptor

// Parameters for the simulation setup
const int maxIter = 2000;
const int nx = 30;
const int ny = 30;
const int nz = 30;

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

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

// Sets material number for fluid
superGeometry.rename( 0,1 );

// 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,
Dynamics<T, DESCRIPTOR>& bulkDynamics1,
SuperGeometry3D<T>& superGeometry ) {

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

// Initial conditions

std::vector<T> v( 3,T() );
AnalyticalConst3D<T,T> zeroVelocity( v );
AnalyticalConst3D<T,T> Rho( 1. );

// Initialize force

std::vector<T> force(3, T());
force[0] = 0;
force[1] = 0;
force[2] = 0.00075;

//olb::BlockVTIreader3D< T,BaseType > datareader(“forcefield.vti”,”ForceField”);
//BlockDataF3D<T, BaseType> blockfield(datareader.getBlockData());
//AnalyticalF3D<T,T>* field;
//field = new SpecialAnalyticalFfromBlockF3D<T,T>(blockfield, datareader.getCuboid(),spacing);
AnalyticalConst3D<T,T> constant(force);
sLattice.defineExternalField( superGeometry, 1,
DESCRIPTOR<T>::ExternalField::forceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfForce, constant );

// Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU( superGeometry, 1, Rho, zeroVelocity );
sLattice.iniEquilibrium( superGeometry, 1, Rho, zeroVelocity );

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

//delete field;
//field = 0;
}

// Output to console and files
void getResults( SuperLattice3D<T, DESCRIPTOR>& sLattice, int iT,
SuperGeometry3D<T>& superGeometry, Timer<T>& timer ) {

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

SuperVTMwriter3D<T> vtmWriter( “ionwake3d” );
SuperLatticeDensity3D<T, DESCRIPTOR> density( sLattice );
vtmWriter.addFunctor( density );

const int statIter = 1;
vtmWriter.createMasterFile();
vtmWriter.write();

// Writes output on the console
if ( iT%statIter==0 ) {
// Timer console output
timer.update( iT );
timer.printStep();

// Lattice statistics console output
sLattice.getStatistics().print( iT,iT );
}
}

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

// === 1st Step: Initialization ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout,”main” );

const T omega1 = 0.0041;

// === 2rd Step: Prepare Geometry ===

const int noOfCuboids = 1;

CuboidGeometry3D<T> cuboidGeometry( 0, 0, 0, 1, nx, ny, nz, noOfCuboids );

// Periodic boundaries in x- and y- and z-direction
cuboidGeometry.setPeriodicity( true, true, true );

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

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

prepareGeometry( superGeometry );

// === 3rd Step: Prepare Lattice ===
SuperLattice3D<T, DESCRIPTOR> sLattice( superGeometry );

ForcedMRTdynamics<T, DESCRIPTOR> bulkDynamics1 ( omega1, instances::getBulkMomenta<T, DESCRIPTOR>() );

prepareLattice( sLattice, bulkDynamics1, superGeometry );

// === 4th Step: Main Loop ===
int iT = 0;
clout << “starting simulation…” << endl;
Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();

//std::vector<T> spacing(3, T());
//spacing[0] = 1;
//spacing[1] = 1;
//spacing[2] = 1;
for ( iT = 0; iT < maxIter; ++iT ) {

// === 5th Step: Collide and Stream Execution ===
//olb::BlockVTIreader3D< T,BaseType > datareader(“forcefield.vti”,”ForceField”);
//BlockDataF3D<T, BaseType> blockfield(datareader.getBlockData());
//AnalyticalF3D<T,T>* field;
//field = new SpecialAnalyticalFfromBlockF3D<T,T>(blockfield, datareader.getCuboid(),spacing);
//sLattice.defineExternalField( superGeometry, 1,
// DESCRIPTOR<T>::ExternalField::forceBeginsAt,
// DESCRIPTOR<T>::ExternalField::sizeOfForce, *field );
//sLattice.collideAndStream();
//delete field;
//field = 0;

// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, iT, superGeometry, timer );
}

timer.stop();
timer.printSummary();
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

Best regards,
Kolotinskiy