VTK Output - Average Density

Palabos General
#1

I am trying to get the “average density” to be written to the VTK file for the code listed below. I have made some attempts but am having trouble and get a ‘unary *’ error if I add:

vtkOut.writeData(*computeAverageDensity(lattice), “AVGdensity”, 1.);

or this error:

In function void writeVTK(plb::MultiBlockLattice2D<T, plb::descriptors::ForcedShanChenD2Q9Descriptor>&, plb::plint)': LiquidBlob.cpp:93: error: no matching function for call tocomputeAverage(plb::MultiScalarField2D&, const char[11], double)’
/home/DLRoseum/palabos/src/dataProcessors/dataAnalysisWrapper2D.hh:1171: note: candidates are: T plb::computeAverage(plb::MultiScalarField2D&, plb::MultiScalarField2D&, int) [with T = T]

if I add:

vtkOut.writeData(*computeAverage(*computeDensity(lattice)), “AVGdensity”, 1.);

to the void writeVTK(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter) section of code.

Please help!!!

Thanks,

Dan


#include "palabos2D.h"
#include "palabos2D.hh"
#include <cstdlib>
#include <iostream>
#include <iomanip>


using namespace plb;
using namespace std;

typedef double T;

#define DESCRIPTOR descriptors::ForcedShanChenD2Q9Descriptor

const plint maxIter = 40001; // Number of iteration steps.
const plint nx = 600;       // Choice of lattice dimensions.
const plint ny = 600;
const T omega = 1.;        // Choice of the relaxation parameter

// All cells have a uniform initially density to represent a vacuum or a liquid ...
const T RhoAll = 75;
// .. except for those in the center circle which have a different value
const T RhoLiquid = 400; //Density in the center circle
//T RhoLiquid2 = 400; //2nd density for liquid - not used right now.

//external body force (gravity) acting on the fluid
const T force = 0;//1.5/(T)ny;  //used to be 0.15

Array<T,2> u0(0,0);

void initializeConstRho(plint iX, plint iY, T& rho, Array<T,2>& u) {
    u = u0;
    rho = RhoAll;
}

void initializeRhoOnCircle(plint iX, plint iY, T& rho, Array<T,2>& u) {
    plint radius = nx/6;
    plint centerX = nx/2;
    plint centerY = ny/2;
    u = u0;
    if( (iX-centerX)*(iX-centerX) + (iY-centerY)*(iY-centerY) < radius*radius) {
        rho = RhoLiquid;
    }
    else {
        rho = RhoAll;
    }
}

// Initialize the lattice at zero velocity and constant density, except
//   for a slight density excess on a circular sub-domain.
void defineInitialDensityAtCenter(MultiBlockLattice2D<T,DESCRIPTOR>& lattice)
{
    // Initialize constant density everywhere.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), RhoAll, u0 );

    // And slightly higher density in the central box.
    initializeAtEquilibrium (
           lattice, lattice.getBoundingBox(), initializeRhoOnCircle );

    // Let's have gravity acting on the fluid.
    setExternalVector(lattice, lattice.getBoundingBox(),
                      DESCRIPTOR<T>::ExternalField::forceBeginsAt, Array<T,2>(0.,-force));

    lattice.initialize();
}

// Write a VTK file which can be post-processed for example in ParaView.
void writeVTK(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
{
    VtkImageOutput2D<T> vtkOut(createFileName("vtk", iter, 6));
    vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);
}

int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/75_400/");

    // For the choice of the parameters G, rho0, and psi0, we refer to the book
    //   Michael C. Sukop and Daniel T. Thorne (2006),
    //   Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
    //   Springer-Verlag Berlin/Heidelberg.

    const int saveIter = 100;
    const int statIter = 100;

    const T G      = -120.0;

    const T rho0 = 200.0;
    const T psi0 = 4.0;

    MultiBlockLattice2D<T, DESCRIPTOR> lattice (
            nx,ny, new ExternalMomentBGKdynamics<T, DESCRIPTOR>(omega) );

    //Periodic horizontal boundary condition
    //lattice.periodicity().toggle(0,true);

    //Periodic boundary conditions on all sides
    lattice.periodicity().toggleAll(true);

    defineInitialDensityAtCenter(lattice);

    // Add the data processor which implements the Shan/Chen interaction potential.
    plint processorLevel = 1;
    integrateProcessingFunctional (
            new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
                G, new interparticlePotential::PsiShanChen94<T>(psi0,rho0) ),
            lattice.getBoundingBox(),
            lattice,
            processorLevel );

    lattice.initialize();

    pcout << "Starting simulation" << endl;
    for (int iT=0; iT<maxIter; ++iT) {
        if (iT%statIter==0) {
            auto_ptr<MultiScalarField2D<T> > rho( computeDensity(lattice) );
            pcout << iT << ": Average rho fluid one = " << computeAverage(*rho) << endl;
            pcout << "Minimum density: " << computeMin(*rho) << endl;
            pcout << "Maximum density: " << computeMax(*rho) << endl;
        }
        if (iT%saveIter == 0) {
            ImageWriter<T>("leeloo.map").writeScaledGif (
                    createFileName("rho", iT, 6), *computeDensity(lattice) );
            writeVTK(lattice, iT);
        }

        lattice.collideAndStream();
    }
}
#2

computeAverage will compute the spatial average of a quantity (the density in your case). This will result in giving you a number. ANd therefore trying to write it in a vtk file will not make much sense. In a vtk file you usually want to write a ScalarField or a TensorField…

#3

Thank you, Malaspin.

Are there any other options for outputting that number into some sort of data file?

Thanks,

Dan

#4

If you want to write the information into a text file you can use plb_ofstream (which is similar to the std::ofstream of c++).