Hi community,
I have a programm at palabos, but I chanhge condition boundary, I had this wrong
*** Process received signal ***
[machine04:80260] Signal: Segmentation fault (11)
[machine04:80260] Signal code: (128)
*******************************************
the change was
boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, top,VelocityWall<T>parameters);
to
boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, top, boundary::freeslip);
Anyone can help me
Hi mfduqued,
I think that setVelocityConditionOnBlockBoundaries takes a boundary type as an argument. So boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, top, boundary::freeslip); seems more correct to me than [...],VelocityWall<T>parameters);. But probably if this was the problem, you would get a compilation error. Can you post a larger section of your code?
Cheers
Mars, Thanks for your time and your response.
At this moment, I want to simulate a Poiseuille flow with slip (by mean populations), that flow has external force,
#define DESCRIPTOR descriptors::ForcedD2Q9Descriptor
template <typename T>
class channelForce {
public:
channelForce(IncomprFlowParam<T> const& parameters_, T& forCe_)
: parameters(parameters_), forCe(forCe_)
{ }
void operator()(plint iX, plint iY, Array<T,2>& force) const {
T Force = forCe;
force[0] = Force;
force[1] = (T)0.0;
}
private:
IncomprFlowParam<T> parameters;
T forCe;
};
template <typename T>
class VelocityWall {
public:
VelocityWall(IncomprFlowParam<T> const& parameters_)
: parameters(parameters_)
{ }
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = (T)0.0;
u[1] = (T)0.0;
}
private:
IncomprFlowParam<T> parameters;
};
template <typename T>
class channelDensityAndVelocityInitial {
public:
channelDensityAndVelocityInitial(IncomprFlowParam<T> const& parameters_)
: parameters(parameters_)
{ }
void operator()(plint iX, plint iY, T &rho, Array<T,2>& u) const {
rho = (T)1.;
u[0] = (T)0.0;
u[1] = (T)0.0;
}
private:
IncomprFlowParam<T> parameters;
};
template<typename T, template<typename U> class Descriptor>
class TopBoundaryFunctional2D
: public BoxProcessingFunctional2D_L<T,Descriptor>
{
public:
TopBoundaryFunctional2D(T& parameterR_)
: parameterR(parameterR_)
{ }
virtual void process(Box2D domain, BlockLattice2D<T,Descriptor>& lattice){
for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
Cell<T,Descriptor>& cell = lattice.get(iX,iY);
Cell<T,Descriptor>& cell1 = lattice.get(iX-1,iY);
Cell<T,Descriptor>& cell2 = lattice.get(iX+1,iY);
Cell<T,Descriptor>& cell3 = lattice.get(iX,iY-1);
cell[3] = ((cell1[7]) * parameterR) + ((cell2[1])*((T)1.0-parameterR));
cell[5] = ((cell2[1]) * parameterR) + ((cell1[7])*((T)1.0-parameterR));
cell[4] = cell3[8];
}
}
}
virtual TopBoundaryFunctional2D<T,Descriptor>* clone() const
{
return new TopBoundaryFunctional2D<T,Descriptor>(*this);
}
virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const {
modified[0] = modif::dataStructure;
}
private:
T parameterR;
};
template<typename T, template<typename U> class Descriptor>
class BottomBoundaryFunctional2D
: public BoxProcessingFunctional2D_L<T,Descriptor>
{
public:
BottomBoundaryFunctional2D(T& parameterR_)
: parameterR(parameterR_)
{ }
virtual void process(Box2D domain, BlockLattice2D<T,Descriptor>& lattice){
for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
Cell<T,Descriptor>& cell = lattice.get(iX,iY);
Cell<T,Descriptor>& cell1 = lattice.get(iX-1,iY);
Cell<T,Descriptor>& cell2 = lattice.get(iX+1,iY);
Cell<T,Descriptor>& cell3 = lattice.get(iX,iY+1);
cell[1] = ((cell1[5]) * parameterR) + ((cell2[3])*((T)1.0-parameterR));
cell[7] = ((cell2[3]) * parameterR) + ((cell1[5])*((T)1.0-parameterR));
cell[8] = cell3[4];
}
}
}
virtual BottomBoundaryFunctional2D<T,Descriptor>* clone() const
{
return new BottomBoundaryFunctional2D<T,Descriptor>(*this);
}
virtual void getTypeOfModification(std::vector<modif::ModifT>& modified)
const {
modified[0] = modif::dataStructure;
}
private:
T parameterR;
};
void createBoundariesAtBottom(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, Box2D domain, T parameterR)
{
integrateProcessingFunctional (new BottomBoundaryFunctional2D<T,DESCRIPTOR>(parameterR), domain, lattice,0);
}
void createBoundariesAtTop(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, Box2D domain, T parameterR)
{
integrateProcessingFunctional (new TopBoundaryFunctional2D<T,DESCRIPTOR>(parameterR), domain, lattice, 0);
};
void channelSetup(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, IncomprFlowParam<T> const& parameters, OnLatticeBoundaryCondition2D<T,DESCRIPTOR>& boundaryCondition, T force, T parameterR)
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
//===================================================
Box2D top = Box2D(0, nx-1, ny-1, ny-1);//ok
Box2D bottom = Box2D(0, nx-1, 0, 0 );//ok
Box2D inlet = Box2D(0, 0, 1, ny-2);//ok
Box2D outlet = Box2D(nx-1, nx-1, 1, ny-2);//ok
//===================================================
Box2D topSlip = Box2D(0, nx-1, ny-2, ny-2);//ok
Box2D bottomSlip = Box2D(0, nx-1, 1, 1 );//ok
//Boundary Consitions
//====================================================
//top
boundaryCondition.setVelocityConditionOnBlockBoundaries ( lattice, outlet, boundary::normalOutflow);
//boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice, top);
// setBoundaryVelocity(lattice, bottom, VelocityWall<T>(parameters));
initializeAtEquilibrium(lattice, lattice.getBoundingBox(),channelDensityAndVelocityInitial<T>(parameters));
//Boundary population topSlilp
createBoundariesAtTop(lattice, top, parameterR);
//bottom
createBoundariesAtBottom(lattice, bottom, parameterR);
//Force
setExternalVector(lattice,Box2D(0,nx-1,1,ny-2),DESCRIPTOR<T>::ExternalField::forceBeginsAt,channelForce<T>(parameters,force));
lattice.initialize();
}
void getParametersFromParamFile(T& Uphysical, T& NuPhysical, T& Tau, T& ULB, T& LenghtPhys, T& lx, T& ly, T& deltaP, T& density, T& convergencia, T& Force, T& maxTime, T& parameterR, plint& parallel, plint& N) {
plb_ifstream ifile("parameters.dat");
ifile >> Uphysical >> NuPhysical >> Tau >> ULB >> LenghtPhys >> lx >> ly >> deltaP >> density >> convergencia >> Force >> maxTime >> parameterR >> parallel >> N ;
global::mpi().bCast(&Uphysical, 1);
global::mpi().bCast(&NuPhysical, 1);
global::mpi().bCast(&Tau, 1);
global::mpi().bCast(&ULB, 1);
global::mpi().bCast(&LenghtPhys, 1);
global::mpi().bCast(&lx, 1);
global::mpi().bCast(&ly, 1);
global::mpi().bCast(&deltaP, 1);
global::mpi().bCast(&density, 1);
global::mpi().bCast(&convergencia, 1);
global::mpi().bCast(&Force, 1);
global::mpi().bCast(&maxTime, 1);
global::mpi().bCast(¶meterR, 1);
global::mpi().bCast(¶llel, 1);
global::mpi().bCast(&N, 1);
}
T computeVmedia (
MultiBlockLattice2D<T,DESCRIPTOR>& lattice, Box2D domain ){
return computeAverage ( *computeVelocityComponent (lattice, domain,0) );
}
T unit_velo(IncomprFlowParam<T> const& parameters, T Uphysical){
return Uphysical * parameters.getDeltaX() / parameters.getDeltaT(); //->ok
}
int main(int argc, char* argv[]) {
plbInit(&argc, &argv);
T Uphysical; //Velocity at physical's units
T NuPhysical; //Viscosity Air 300K, units physical
T ULB; //velocity units LB
T LenghtPhys; // lenght caractheristic
T lx; //times length caractheristic at axis X
T ly; //times length caractheristic at axis Y
T lz; //times length caractheristic at axis Z
T deltaP; //delta P
T density;//density
T convergencia;//convergence
T Force; //Force
T Tau; //Relaxation time
T maxTime; //Maximum iteraction at each node
T parameterR; // fraction de bounce-back
plint parallel;
plint N; //resolution
getParametersFromParamFile(Uphysical, Kn, NuPhysical, Tau, ULB, LenghtPhys, lx, ly, deltaP, density, convergencia, Force, maxTime, parameterR, alpha, parallel, N);
T Re = Uphysical*LenghtPhys/NuPhysical; //Physical Velocity from Reynolds' Number
IncomprFlowParam<T> parameters(
ULB,
Re,
N,
lx,
ly
);
const T maxT = (T)maxTime;
const T imSave = (T)maxT/(T)1000;
const T promedio = (T)maxT/(T)100;
const T checkPointIter = 200;
T iniT, endT;
global::directories().setOutputDir("./tmp/");
global::IOpolicy().activateParallelIO(true);
writeLogFile(parameters, "2D channel Forced");
plint nx = parameters.getNx();
plint ny = parameters.getNy();
T Omega0 = parameter.getOmega();
#define Dynamics GuoExternalForceBGKdynamics<T,DESCRIPTOR>(Omega0)
// Parallelization in x-direction.
plint procNumber = parallel - (T)1.0;
plint yTiles = 1;
plint xTiles = procNumber;
std::auto_ptr<MultiBlockLattice2D<T,DESCRIPTOR> > lattice (
new MultiBlockLattice2D<T,DESCRIPTOR> (
MultiBlockManagement2D (
// Parallelize in x-direction
createRegularDistribution2D (
nx, ny, xTiles, yTiles),
defaultMultiBlockPolicy2D().getThreadAttribution(), 1 ),
defaultMultiBlockPolicy2D().getBlockCommunicator(),
defaultMultiBlockPolicy2D().getCombinedStatistics(),
defaultMultiBlockPolicy2D().getMultiCellAccess<T,DESCRIPTOR>(),
new Dynamics));
// Use periodic boundary conditions - axis x.
lattice->periodicity().toggle(0,true);
// Use periodic boundary conditions - axis y.
lattice->periodicity().toggle(1,false);
OnLatticeBoundaryCondition2D<T,DESCRIPTOR>* boundaryCondition=createLocalBoundaryCondition2D<T,DESCRIPTOR>();
channelSetup(*lattice, parameters, *boundaryCondition, Force, parameterR);
plb_ofstream Results("./tmp/Results.dat");
plb_ofstream successiveProfilesFinal1("./tmp/Profile_Velocity.dat");
plb_ofstream UPvsT("./tmp/VvsT.dat");
Box2D profileSectionInicial(0, 0, 1, ny-2);
Box2D profileSectionMity(0, nx-1, ny/2, ny/2);
Box2D profileSectionFinal1(nx - 1, nx - 1, 1, ny -2);
// Convergence.
util::ValueTracer<T> converge(parameters.getLatticeU(), parameters.getResolution(), convergencia);
T previousIterationTime = T();
global::timer("total").restart();
// Loop over main time iteration.
for (plint iT=0; iT<parameters.nStep(maxT); ++iT) {
global::timer("mainLoop").restart();
counter += 1;
double TimE = iT*parameters.getDeltaT()*unit_time;
UPvsT << setprecision(16)
<< TimE << " " << *computeVelocityNorm (*lattice, profileSectionFinal1) << std::endl;
// Converged
converge.takeValue(getStoredAverageEnergy(*lattice),false);
if (converge.hasConverged())
{
break;
}
// Execute a time iteration.
lattice->collideAndStream();
previousIterationTime = global::timer("mainLoop").stop();
}
for (plint iy=0; iy<ny; iy++){
Box2D profileSectionFinal1(nx - 1, nx - 1, iy, iy);
successiveProfilesFin01 << setprecision(16)
<< iy << " "<< *computeVelocityNorm(*lattice, profileSectionFinal1) << std::endl;
}
Pressure << setprecision(16) << *multiply(*computeDensity (*lattice, profileSectionMity), (T)1.0/(T)3.0) << endl;
delete boundaryCondition;
T U_meanI = computeVmedia(*lattice,profileSectionInicial) ;
T U_meanF = computeVmedia(*lattice,profileSectionFinal1);
Results << setprecision(16)
<< "mean Velocity outlet =" <<" "<< U_meanF << '\n'
<< "mean Velocity inlet =" <<" "<< U_meanI << '\n'
}
This programm worked very good when I used boundaryCondition.setVelocityConditionOnBlockBoundaries, but this option doesn’t show a slip on boundary.
Thanks Mars for your help
I could to solve this problem. It was the cell3 at the code. But, at this moment I have another problem, it is that in the TopBoundaryFunctional2D class, when the program runs it does not depend on the parameter R that is in the class.
I change this parameter, and I always get the same result.
This problem is the cell3, at this moment I don’t know the reason
Hello,
I think that when you do cell[3] = ... you are trying to write in a location that is outside of your accessible region.
Your domain is indeed a plane located at y = ny-1 but you are trying to write in y-1. Now, usually, this does not cause a segmentation fault, because around the domain you have an envelope, but writing in the envelope won’t actually change the value in your simulation. In this case, it is probably this the problem, try to enlarge the dataprocessor domain…
Hi mars,
Thanks for your response.I have another problem as,
template<typename T, template<typename U> class Descriptor>
class TopBoundaryFunctional2D
: public BoxProcessingFunctional2D_L<T,Descriptor>
{
public:
TopBoundaryFunctional2D(T parameterR_)
: parameterR(parameterR_)
{ }
virtual void process(Box2D domain, BlockLattice2D<T,Descriptor>& lattice){
for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
Cell<T,Descriptor>& cell = lattice.get(iX,iY);
Cell<T,Descriptor>& cell1 = lattice.get(iX-1,iY);
Cell<T,Descriptor>& cell2 = lattice.get(iX+1,iY);
Cell<T,Descriptor>& cell3 = lattice.get(iX,iY-1);
cell[3] = ((cell1[7]) * parameter) + ((cell2[1])*((T)1.0-paramete$
cell[5] = ((cell2[1]) * parameter) + ((cell1[7])*((T)1.0-paramete$
cell[4] = cell[8];
}
}
}
virtual TopBoundaryFunctional2D<T,Descriptor>* clone() const
{
return new TopBoundaryFunctional2D<T,Descriptor>(*this);
}
virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) cons$
modified[0] = modif::dataStructure;
}
private:
T parameterR;
};
void createBoundariesAtTop(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, Box2D domain, T& parameterR)
{
integrateProcessingFunctional (new TopBoundaryFunctional2D<T,DESCRIPTOR>(parameterR), domain, lattice,0);
}
I read the pararameterR from a file, but it doesn’t affect the above class and I obtain the same results for different values of parameterR
Thanks for your help
Hello,
first. When you encounter such problems do not hesitate either to use a debugger (like gdb) or a sanitizer (see here https://stackoverflow.com/questions/37970758/how-to-use-addresssanitizer-in-gcc). It is usually of great help.
Second. from what I read here you are trying to read from the points iX-1, iX+1, iY-1, iY+1. Are you sure these points are in the domain and that you are not trying to read form points that are located outside of the computational domain?