Hello Everyone,
I am trying to change the Womersley example to a 3D case. When I compile the cpp file I get an error: Assertion failed: vectorStartsAt+Descriptor::d <= Descriptor::ExternalField::numScalars
I would be really thankful if someone help with that.
Here is my code for Womersley 3D.
#include "palabos3D.h"
#include "palabos3D.hh"
#include <cstdlib>
#include <iostream>
#include <complex>
using namespace plb;
using namespace std;
typedef double T;
const T pi = (T)4.0 * std::atan((T)1.0);
#define DESCRIPTOR descriptors::D3Q19Descriptor
#define DYNAMICS GuoExternalForceBGKdynamics<T, DESCRIPTOR>(omega)
// #define DYNAMICS HeExternalForceBGKdynamics<T, NSDESCRIPTOR>(omega)
// #define DYNAMICS ShanExternalForceBGKdynamics<T, NSDESCRIPTOR>(omega)
/// Velocity on the parabolic Womersley profile
T womersleyVelocity(plint iY, plint iZ, T t, T A, T omega, T alpha, T ly, IncomprFlowParam<T> const& parameters)
{
const complex<T> I(0., 1.);
T y = (T)iY / parameters.getResolution();
T z = (T)iZ / parameters.getResolution();
return (std::exp(I * omega * t)*((T)1./(alpha*alpha) + A*std::cosh(alpha*z) * A*std::cos(pi*y) *
(std::cosh(std::sqrt(pi*pi+alpha*alpha)*z))/(std::cosh(std::sqrt(pi*pi+alpha*alpha)*ly)))).real();
}
/// Time dependent but space-independent Womersley force.
T womersleyForce(T t, T A, T omega, IncomprFlowParam<T> const& parameters)
{
return A * std::cos(omega * t);
}
/// A functional, used to initialize the velocity for the boundary conditions
template<typename T>
class WomersleyVelocity {
public:
WomersleyVelocity(IncomprFlowParam<T> parameters_, T alpha_, T t_)
: parameters(parameters_), alpha(alpha_),
omega((T)4 * util::sqr(alpha) * parameters.getLatticeNu() / (T)(util::sqr(parameters.getResolution()))),
A(8. * parameters.getLatticeNu() * parameters.getLatticeU() / ((T)(util::sqr(parameters.getResolution())))),
t(t_)
{ }
void operator()(plint iX, plint iY, plint iZ, Array<T, 3>& u) const {
u[0] = womersleyVelocity(iY, iZ, t, A, omega, alpha, ly, parameters);
u[1] = T();
u[2] = T();
}
private:
IncomprFlowParam<T> parameters;
T alpha, omega, A;
T t;
};
void channelSetup(MultiBlockLattice3D<T, DESCRIPTOR>& lattice,
IncomprFlowParam<T> const& parameters,
OnLatticeBoundaryCondition3D<T, DESCRIPTOR>& boundaryCondition,
T alpha, T frequency, T amplitude)
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();
Box3D top = Box3D(0, nx - 1, ny - 1, ny - 1, 0, nz - 1);
Box3D bottom = Box3D(0, nx - 1, 0, 0, 0, nz - 1);
Box3D left = Box3D(0, 0, 1, ny - 2, 0, nz - 1);
Box3D right = Box3D(nx - 1, nx - 1, 1, ny - 2, 0, nz - 1);
boundaryCondition.addVelocityBoundary1N(bottom, lattice);
boundaryCondition.addPressureBoundary1P(top, lattice);
boundaryCondition.addVelocityBoundary2P(right, lattice);
boundaryCondition.addPressureBoundary2N(left, lattice);
Array<T, 3> u((T)0., (T)0., (T)0.);
setBoundaryVelocity(lattice, lattice.getBoundingBox(), u);
initializeAtEquilibrium(lattice, lattice.getBoundingBox(), (T)1.0, u);
Array<T, DESCRIPTOR<T>::d> force(womersleyForce((T)0, amplitude, frequency, parameters), 0., 0.);
setExternalVector(lattice, lattice.getBoundingBox(), DESCRIPTOR<T>::ExternalField::forceBeginsAt, force);
lattice.initialize();
}
template<class BlockLatticeT>
void writeVTK(BlockLatticeT& lattice,
IncomprFlowParam<T> const& parameters, plint iter)
{
T dx = parameters.getDeltaX();
T dt = parameters.getDeltaT();
VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx);
vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx / dt);
}
template<class BlockLatticeT>
void writeGif(BlockLatticeT& lattice, plint iT)
{
const plint imSize = 600;
ImageWriter<T> imageWriter("leeloo.map");
imageWriter.writeScaledGif(createFileName("u", iT, 6),
*computeVelocityNorm(lattice),
imSize, imSize);
}
T computeRMSerror(MultiBlockLattice3D<T, DESCRIPTOR>& lattice,
IncomprFlowParam<T> const& parameters,
T alpha, plint iT, bool createImage = false)
{
MultiTensorField3D<T, 3> analyticalVelocity(lattice);
setToFunction(analyticalVelocity, analyticalVelocity.getBoundingBox(),
WomersleyVelocity<T>(parameters, alpha, (T)iT));
MultiTensorField3D<T, 3> numericalVelocity(lattice);
computeVelocity(lattice, numericalVelocity, lattice.getBoundingBox());
// Divide by lattice velocity to normalize the error
return 1. / parameters.getLatticeU() *
// Compute RMS difference between analytical and numerical solution
std::sqrt(computeAverage(*computeNormSqr(
*subtract(analyticalVelocity, numericalVelocity)
)));
}
int main(int argc, char* argv[])
{
plbInit(&argc, &argv);
if (argc != 2)
{
pcout << "Error : Wrong parameters specified." << endl;
pcout << "1 : N." << endl;
exit(1);
}
const plint N = atoi(argv[1]);
const T Re = 1.0;
const T alpha = 1.0; // womersley number
const plint Nref = 10;
const T uMaxRef = 0.01;
const T uMax = uMaxRef / (T)N * (T)Nref; // needed to avoid compressibility errors.
const T lx = 10.0;
const T ly = 2.0;
const T lz = 2.0;
pcout << "uMaxRef=" << uMaxRef << std::endl;
pcout << "uMax=" << uMax << std::endl;
global::directories().setOutputDir("C:/Users/Asia/Documents/0Magisterka/palabos/examples/showCases/womersley/tmp/wom3d/");
IncomprFlowParam<T> parameters(uMax, Re, N, lx, ly, lz);
// The frequency of the force (lattice units)
T frequency = (T)4 * alpha * alpha * parameters.getLatticeNu()
/ (T)(parameters.getResolution() * parameters.getResolution());
// The amplitude of the forcing term (lattice units)
T amplitude = 8. * parameters.getLatticeNu() * parameters.getLatticeU()
/ ((T)(parameters.getResolution() * parameters.getResolution()));
// Period of the force (lattice units)
plint tPeriod = (plint)((T)2 * pi / frequency + 0.5);
writeLogFile(parameters, "palabos.log");
plint nx = parameters.getNx();
plint ny = parameters.getNy();
plint nz = parameters.getNz();
T omega = parameters.getOmega();
MultiBlockLattice3D<T, DESCRIPTOR> lattice(
nx, ny, nz, new DYNAMICS);
OnLatticeBoundaryCondition3D<T, DESCRIPTOR>*
boundaryCondition = createLocalBoundaryCondition3D<T, DESCRIPTOR>();
lattice.periodicity().toggle(0, true);
channelSetup(lattice, parameters, *boundaryCondition, alpha, frequency, amplitude);
pcout << "Starting simulation" << endl;
const plint maxIter = tPeriod * 100;
const plint tSave = tPeriod / 24;
const plint vtkSave = tPeriod;
T error = T();
lattice.resetTime(1);
pcout << "Omega = " << omega << ", it period = " << tPeriod << endl;
util::ValueTracer<T> converge(uMax, N, 1.0e-3);
plint iT = 0;
for (iT = 0; iT < maxIter; ++iT) {
// Updating the force in the whole domain
Array<T, DESCRIPTOR<T>::d> force(womersleyForce((T)iT, amplitude, frequency, parameters), 0., 0.);
setExternalVector(lattice, lattice.getBoundingBox(),
DESCRIPTOR<T>::ExternalField::forceBeginsAt, force);
T errorTemp = computeRMSerror(lattice, parameters, alpha, iT);
error += errorTemp;
/*
if (iT % tSave == 0) {
pcout << "Writing Gif at time : " << iT << std::endl;
writeGif(lattice,iT);
}*/
if (iT % vtkSave == 0 && iT > 0) {
pcout << "Saving VTK file at time : " << iT << std::endl;
writeVTK(lattice, parameters, iT);
}
if (iT % tPeriod == 0)
{
// The error is averaged over one period
error /= (T)(tPeriod);
pcout << "For N = " << N << ", Error = " << error << endl;
converge.takeValue(error, true);
if (converge.hasConverged())
{
cout << "Simulation converged!\n";
break;
}
error = T();
}
lattice.collideAndStream();
}
pcout << "For N = " << N << ", Error = " << computeRMSerror(lattice, parameters, alpha, iT, true) << endl;
}