Ok, I learned that I can’t access coordinates from boundary condition dynamics, so I switched to data processors.
I wrote a simple OneCellIndexed functional that does nothing and am trying to run it with each iteration, but I keep getting segmentation fault. Is there any method that allows me to define an indexed functional that’d be run with each iteration? (something like integrateProcessingFunctional?) Or should I create a new functional object for each iteration?
Another thing I’d like to do is to disable streaming between 2 columns (in this case nx/2 and nx/2+1). Can it be done? Or should I define 2 separate lattices and apply functional to the end of one and beginning of another?
#include "palabos2D.h"
#include "palabos2D.hh"
#include <cstdlib>
#include <iostream>
using namespace plb;
using namespace std;
typedef double T;
#define DESCRIPTOR descriptors::D2Q9Descriptor
template<typename T, template<typename U> class Descriptor>
class Gate : public OneCellIndexedFunctional2D<T,Descriptor> {
public:
Gate(plint ny_, bool left_)
: ny(ny_),
left(left_)
{ }
Gate<T,Descriptor>* clone() const {
return new Gate<T,Descriptor>(*this);
}
virtual void execute(plint iX, plint iY, Cell<T,Descriptor>& cell) const {
// T h = 0.;
// for (plint i=0; i<9; i++)
// h += cell[i];
// cout << "Height in (" << iX << "," << iY << ") = " << h << endl;
}
private:
plint ny;
bool left;
};
Gate<T, DESCRIPTOR> *gleft, *gright;
void canalSetup( MultiBlockLattice2D<T, DESCRIPTOR>& water,
T rho0,
T force )
{
int nx = water.getNx();
int ny = water.getNy();
// Bounce-back on bottom wall
defineDynamics(water, Box2D(0,nx-1, 0,0), new BounceBack<T, DESCRIPTOR>(rho0) );
// Bounce-back on top wall
defineDynamics(water, Box2D(0,nx-1, ny-1,ny-1), new BounceBack<T, DESCRIPTOR>(rho0) );
// External force
// setExternalVector(water, water.getBoundingBox(),DESCRIPTOR<T>::ExternalField::forceBeginsAt, Array<T,2>(0.,-force));
water.initialize();
}
int main(int argc, char *argv[])
{
plbInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");
const T omega = 1.0;
const int nx = 1200;
const int ny = 400;
T force = 0.15/(T)ny;
T rho0 = 0.;
const int maxIter = 16000;
const int saveIter = 100;
MultiBlockLattice2D<T, DESCRIPTOR> water (nx,ny, new BGKdynamics<T, DESCRIPTOR>(omega) );
water.periodicity().toggleAll(false);
canalSetup(water, rho0, force);
pcout << "Starting simulation" << endl;
for (int iT=0; iT<maxIter; ++iT) {
if (iT%saveIter==0) {
ImageWriter<T> imageWriter("leeloo");
imageWriter.writeScaledGif(createFileName("u", iT, 6), *computeVelocityNorm(water));
}
water.collide();
gleft = new Gate<T, DESCRIPTOR>(ny,true);
gright = new Gate<T, DESCRIPTOR>(ny,false);
applyIndexed(water, Box2D(nx/2, nx/2, 0, ny-1), gleft);
applyIndexed(water, Box2D(nx/2+1, nx/2+1, 0, ny-1), gright);
water.stream();
}
}