An issue when running parallel simulations

[avrachan1]

Dear Community,

I define two lattices from the same geometry.

SuperLattice<T,DESCRIPTOR> A(superGeometry);
SuperLattice<T,DESCRIPTOR> B(superGeometry);

I define the dynamics for both of the lattices.
A.defineDynamics<BGKdynamics>(bulkIndicator);
B.defineDynamics<BGKdynamics>(bulkIndicator);

where bulkIndicater indicates the whole region – which is a simple cuboid with all regions as one material.

I define *ONE* of the lattice’s density and velocity using some AnalyticalF3D function.

phaseField.defineRhoU(bulkIndicator,tanhs,uF);
phaseField.iniEquilibrium(bulkIndicator,tanhs,uF);

I initialize
A.initialize();

get the populations on A
SuperLatticeField3D <T,DESCRIPTOR,POPULATION> populationA(A);
AnalyticalFfromSuperF3D <T> analyticalFpopA(populationA,true,true);

Use this to define the population of B,

B.definePopulations(bulkIndicator,analyticalFpopA);

However when I inspect the populations of B, it is different depending on the number of
cores used.

For example I get what I expect when I use “mpirun -np 1 …”

But the results change when I use 2 or 4 cores etc.

What could be the reason ?

#################### CODE TO REPRODUCE ISSUE ##############################

#include “olb3D.h”
#include “olb3D.hh”

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

using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = D3Q15<>;

//const T flow_vel = 4.0;
//const T radiusCylinder = 0.05;
//const T viscosity = 0.2e-3; // This is Kappa for heat equation

const int N = 10;
//const T char_phys_l = 1.0;

const T W0 = 0.01;
const T tau0 = 1e-4;

const T dx= 0.01;
const T dt= 1.5*1e-5;

const T lengthX = 51*dx;
const T lengthY = 51*dx;
const T lengthZ = 51*dx;

//const T eta_theta = 0.535; // This is eta_theta
//const T eta_phi = 0.55; // this is eta_phi

const T kappa = 0.7;
const T lambda = 10;
const T eta_theta = descriptors::invCs2<T,DESCRIPTOR>()*kappa*(dt/(dx*dx))+1./2.; // This is eta_theta
const T eta_phi = descriptors::invCs2<T,DESCRIPTOR>()*W0*W0/tau0*(dt/(dx*dx))+1./2.;//0.55; // this is eta_phi

template <typename T, typename S, typename DESCRIPTOR>
class TanhSphere : public AnalyticalF3D<T, S>
{
private:
T centerx,centery,centerz,radius, ls; // ls is the width of the border .
//UnitConverter<T, DESCRIPTOR> const &converter;
//int iT, tempCase;

public:
TanhSphere(T centerx_, T centery_,T centerz_ , T radius_, T ls_) : AnalyticalF3D<T, S>(3),
centerx(centerx_), centery(centery_), centerz(centerz_), radius(radius_), ls(ls_)
{
this->getName() = “TanhSphere”;
};

bool operator()(T output[1], const S x[3])
{
//T res = converter.getResolution();
T real_distx = x[0];//*converter.getPhysDeltaX();
T real_disty = x[1];//*converter.getPhysDeltaX();
T real_distz = x[2];//*converter.getPhysDeltaX();
T distance = sqrt((real_distz-centerz)*(real_distz-centerz)+(real_distx-centerx)*(real_distx-centerx)+(real_disty-centery)*(real_disty-centery));
T argument = (radius – distance)/sqrt(ls);

output[0] = util::tanh(argument);
//std::cout<<x[0]<<“\t”<<x[1]<<“\t”<<distance<<“\t”<<output[0]<<“\n”;

return true;
};
};
int main( int argc, char* argv[] )
{
olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./AB_error_check/” );
OstreamManager clout(std::cout,”main”);

clout<<“Simulation start..”<<std::endl;
UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converter(
int {1}, // resolution: number of voxels per charPhysL
(T) eta_phi, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) dx*1, // charPhysLength: reference length of simulation geometry
(T) 0.02*(dx/dt), // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) W0*W0/tau0, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 1.0 // physDensity: physical density in __kg / m^3__
);
converter.print();

T L = converter.getPhysDeltaX();
Vector<T,3> extend(lengthX,lengthY,lengthZ);
Vector<T,3> center(lengthX/2.,lengthY/2.,lengthZ/2.);
Vector<T,3> origin;

IndicatorCuboid3D<T> region(extend,origin);

#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = singleton::mpi().getSize();
#else
const int noOfCuboids = 1;
#endif
CuboidGeometry3D<T> regionGeometry(region,L,noOfCuboids);
regionGeometry.setPeriodicity(true,true,true);

HeuristicLoadBalancer<T> loadBalancer( regionGeometry );
SuperGeometry<T,3> superGeometry(regionGeometry,loadBalancer);

superGeometry.rename(0,1);

superGeometry.clean(true,{1});

superGeometry.checkForErrors();
superGeometry.print();

SuperLattice<T,DESCRIPTOR> A(superGeometry);
SuperLattice<T,DESCRIPTOR> B(superGeometry);

auto bulkIndicator = superGeometry.getMaterialIndicator({1});

TanhSphere <T,T,DESCRIPTOR> tanhs(center[0],center[1],center[2],8*dx,2.*W0);

A.defineDynamics<BGKdynamics>(bulkIndicator);
B.defineDynamics<BGKdynamics>(bulkIndicator);

std::vector<T> velocity(3,T(0.));

AnalyticalConst3D<T,T> uF(velocity);

A.defineRhoU(bulkIndicator,tanhs,uF);
A.iniEquilibrium(bulkIndicator,tanhs,uF);

A.initialize();
B.initialize();

SuperLatticeField3D <T,DESCRIPTOR,POPULATION> populationA(A);

AnalyticalFfromSuperF3D <T> analyticalFpopA(populationA,true,true);

SuperLatticeField3D <T,DESCRIPTOR,POPULATION> populationB(B);

SuperVTMwriter3D<T> vtmWriterPF(“A”);
SuperVTMwriter3D<T> vtmWriterPFm1(“B”);

vtmWriterPF.createMasterFile();
vtmWriterPFm1.createMasterFile();

vtmWriterPF.addFunctor(populationA);
vtmWriterPFm1.addFunctor(populationB);

B.definePopulations(bulkIndicator,analyticalFpopA);

vtmWriterPF.write(-1);
vtmWriterPFm1.write(-1);

return 0;
}

[avrachan1]

I think I have identified the issue to be

AnalyticalFfromSuperF3D <T> analyticalFpopA(populationA,true,true);

I get the correct results when I used

AnalyticalFfromSuperF3D <T> analyticalFpopA(populationA);

[Adrian]

Yes, this issue was caused by enabling the “communicate to all”-mode with the second constructor argument. Happy to hear that the default arguments worked for you.

In any case: Be aware that this functor will interpolate the values, this is not necessary when sharing the same discretization (as you do here due to constructing both lattices from the same super geometry).

[avrachan1]

Hi Adrian,

Thanks for the reply.

I am aware that it interpolates values, but I could not find an alternate way to define populations on the new lattice. It would be great if you can describe it!

Also if you can elaborate on what does “communicate-to-all” mode does here?

Thanks!

[Adrian]

You can use the SuperF-accepting overload of SuperLattice::definePopulations instead. If the discretization is identical (as it needs to be for coupling two lattices) you can directly pass your populationA functor.

Communicate-to-all communicates the results of all evaluations to all processes. e.g. this can be used if all processes need to iterate over the whole simulation domain (but not all processes hold all data, i.e. parallelization is used).

[avrachan1]

Dear Adrian,

Thanks!

When I replace

B.definePopulations(bulkIndicator,analyticalFpopA);
with
B.definePopulations(bulkIndicator,populationPF);
I am getting errors, what is the correct way to define populations directly from SuperF ?

####################

1 olb-1.6r0/src/core/blockLattice.hh:142:9: error: no matching function for call to object of type ‘BlockF<double, D2Q5<>::d>’ (aka ‘olb::BlockF2D<double>’)
2 popF(pop, loc);
3 ^~~~
4 olb-1.6r0/src/core/superLattice.hh:321:17: note: in instantiation of member function ‘olb::BlockLattice<double, olb::descriptors::D2Q5<>>::definePopulations’ requested here
5 _block[iC]->definePopulations(indicator->getBlockIndicatorF(iC), Pop.getBlockF(iC));
6 ^
7 error.cpp:175:4: note: in instantiation of member function ‘olb::SuperLattice<double, olb::descriptors::D2Q5<>>::definePopulations’ requested here
8 B.definePopulations(bulkIndicator,populationPF);
9 ^
10 olb-1.6r0/src/functors/genericF.h:81:16: note: candidate function not viable: no known conversion from ‘LatticeR<D2Q5<>::d>’ (aka ‘Vector<int, 2U>’) to ‘const int *’ for 2nd argument
11 virtual bool operator() (T output[], const S input[])=0;
12 ^
13 olb-1.6r0/src/functors/genericF.h:87:8: note: candidate function not viable: no known conversion from ‘LatticeR<D2Q5<>::d>’ (aka ‘Vector<int, 2U>’) to ‘int’ for 2nd argument
14 bool operator() (T output[], S input0);
15 ^
16 olb-1.6r0/src/functors/genericF.h:86:8: note: candidate function not viable: requires single argument ‘output’, but 2 arguments were provided
17 bool operator() (T output[]);
18 ^
19 olb-1.6r0/src/functors/genericF.h:88:8: note: candidate function not viable: requires 3 arguments, but 2 were provided
20 bool operator() (T output[], S input0, S input1);
21 ^
22 olb-1.6r0/src/functors/genericF.h:89:8: note: candidate function not viable: requires 4 arguments, but 2 were provided
23 bool operator() (T output[], S input0, S input1, S input2);
24 ^
25 olb-1.6r0/src/functors/genericF.h:90:8: note: candidate function not viable: requires 5 arguments, but 2 were provided
26 bool operator() (T output[], S input0, S input1, S input2, S input3);
27 ^
28 2 warnings and 1 error generated.
29 make: *** [error.o] Error 1
30

[Adrian]

What is the type of populationPF?

[avrachan1]

Sorry my mistake,

SuperLatticeField3D <T,DESCRIPTOR,POPULATION> populationPF(A);

followed by

B.definePopulations(bulkIndicator,populationPF);

leads to the error.

[Adrian]

Ok, this is indeed a bug in the code. You should be able to fix it by changing line 142 of src/core/blockLattice.hh to:

        popF(pop, loc.data());

or using SuperLattice::defineField<POPULATION>

Sorry about this and thanks for discovering this issue!

I also pushed a fix to the public branch.

[Author]

Posts