~registry/+junk/palabos-packaging : revision 1

1

#!/bin/sh

2

# autopkgtest check

3

4

5

set -e

6

7

8

WORKDIR=$(mktemp -d)

9

trap "rm -rf $WORKDIR" 0 INT QUIT ABRT PIPE TERM

10

cd $WORKDIR

11

mkdir src

12

cd src

13

14

cat <<EOF > CMakeLists.txt

15

cmake_minimum_required(VERSION 2.6)

16

project(Test)

17

find_package(Palabos REQUIRED)

18

include_directories(\${PALABOS_INCLUDE_DIR})

19

20

add_executable(demo demo.cpp)

21

target_link_libraries(demo \${PALABOS_LIBRARIES})

22

install(TARGETS demo DESTINATION bin)

23

24

EOF

25

26

cat <<EOF > demo.cpp

27

/* This file is part of the Palabos library.

28

*

29

30

* Route d'Oron 2

31

* 1010 Lausanne, Switzerland

32

* E-mail contact: contact@flowkit.com

33

*

34

* The most recent release of Palabos can be downloaded at

35

* <http://www.palabos.org/>

36

*

37

* The library Palabos is free software: you can redistribute it and/or

38

* modify it under the terms of the GNU Affero General Public License as

39

* published by the Free Software Foundation, either version 3 of the

40

* License, or (at your option) any later version.

41

*

42

* The library is distributed in the hope that it will be useful,

43

* but WITHOUT ANY WARRANTY; without even the implied warranty of

44

* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

45

* GNU Affero General Public License for more details.

46

*

47

* You should have received a copy of the GNU Affero General Public License

48

* along with this program. If not, see <http://www.gnu.org/licenses/>.

49

*/

50

51

/** \file

52

* Flow in a lid-driven 3D cavity. The cavity is square and has no-slip walls,

53

* except for the top wall which is diagonally driven with a constant

54

* velocity. The benchmark is challenging because of the velocity

55

* discontinuities on corner nodes.

56

**/

57

58

#include "palabos3D.h"

59

#include "palabos3D.hh"

60

#include <vector>

61

#include <cmath>

62

#include <iostream>

63

#include <fstream>

64

65

using namespace plb;

66

using namespace std;

67

68

typedef double T;

69

#define DESCRIPTOR descriptors::D3Q19Descriptor

70

71

void cavitySetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,

72

IncomprFlowParam<T> const& parameters,

73

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition )

74

{

75

const plint nx = parameters.getNx();

76

const plint ny = parameters.getNy();

77

const plint nz = parameters.getNz();

78

Box3D topLid = Box3D(0, nx-1, ny-1, ny-1, 0, nz-1);

79

Box3D everythingButTopLid = Box3D(0, nx-1, 0, ny-2, 0, nz-1);

80

81

// All walls implement a Dirichlet velocity condition.

82

boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

83

84

T u = std::sqrt((T)2)/(T)2 * parameters.getLatticeU();

85

initializeAtEquilibrium(lattice, everythingButTopLid, (T)1., Array<T,3>((T)0.,(T)0.,(T)0.) );

86

initializeAtEquilibrium(lattice, topLid, (T)1., Array<T,3>(u,(T)0.,u) );

87

setBoundaryVelocity(lattice, topLid, Array<T,3>(u,(T)0.,u) );

88

89

lattice.initialize();

90

}

91

92

template<class BlockLatticeT>

93

void writeGifs(BlockLatticeT& lattice,

94

IncomprFlowParam<T> const& parameters, plint iter)

95

{

96

const plint imSize = 600;

97

const plint nx = parameters.getNx();

98

const plint ny = parameters.getNy();

99

const plint nz = parameters.getNz();

100

const plint zComponent = 2;

101

102

Box3D slice(0, nx-1, 0, ny-1, nz/2, nz/2);

103

ImageWriter<T> imageWriter("leeloo");

104

105

imageWriter.writeScaledGif( createFileName("uz", iter, 6),

106

*computeVelocityComponent (lattice, slice, zComponent),

107

imSize, imSize );

108

109

imageWriter.writeScaledGif( createFileName("uNorm", iter, 6),

110

*computeVelocityNorm (lattice, slice),

111

imSize, imSize );

112

imageWriter.writeScaledGif( createFileName("omega", iter, 6),

113

*computeNorm(*computeVorticity (

114

*computeVelocity(lattice) ), slice ),

115

imSize, imSize );

116

}

117

118

template<class BlockLatticeT>

119

void writeVTK(BlockLatticeT& lattice,

120

IncomprFlowParam<T> const& parameters, plint iter)

121

{

122

T dx = parameters.getDeltaX();

123

T dt = parameters.getDeltaT();

124

VtkImageOutput3D<T> vtkOut(createFileName("vtk", iter, 6), dx);

125

vtkOut.writeData<float>(*computeVelocityNorm(lattice), "velocityNorm", dx/dt);

126

vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", dx/dt);

127

vtkOut.writeData<3,float>(*computeVorticity(*computeVelocity(lattice)), "vorticity", 1./dt);

128

}

129

130

131

SparseBlockStructure3D createRegularDistribution3D (

132

std::vector<plint> const& xVal, std::vector<plint> const& yVal, std::vector<plint> const& zVal )

133

{

134

PLB_ASSERT(xVal.size()>=2);

135

PLB_ASSERT(yVal.size()>=2);

136

PLB_ASSERT(zVal.size()>=2);

137

SparseBlockStructure3D dataGeometry (

138

Box3D(xVal[0], xVal.back()-1, yVal[0], yVal.back()-1, zVal[0], zVal.back()-1) );

139

for (plint iX=0; iX<(plint)xVal.size()-1; ++iX) {

140

for (plint iY=0; iY<(plint)yVal.size()-1; ++iY) {

141

for (plint iZ=0; iZ<(plint)zVal.size()-1; ++iZ) {

142

dataGeometry.addBlock (

143

Box3D( xVal[iX], xVal[iX+1]-1, yVal[iY],

144

yVal[iY+1]-1, zVal[iZ], zVal[iZ+1]-1 ),

145

dataGeometry.nextIncrementalId() );

146

}

147

}

148

}

149

return dataGeometry;

150

}

151

152

SparseBlockStructure3D createCavityDistribution3D(plint nx, plint ny, plint nz)

153

{

154

static const plint numval=4;

155

plint x[numval] = {0, nx/4, 3*nx/4, nx};

156

plint y[numval] = {0, ny/4, 3*ny/4, ny};

157

plint z[numval] = {0, nz/4, 3*nz/4, nz};

158

std::vector<plint> xVal(x, x+numval);

159

std::vector<plint> yVal(y, y+numval);

160

std::vector<plint> zVal(z, z+numval);

161

return createRegularDistribution3D(xVal, yVal, zVal);

162

}

163

164

165

int main(int argc, char* argv[]) {

166

167

plbInit(&argc, &argv);

168

global::directories().setOutputDir("./tmp/");

169

170

IncomprFlowParam<T> parameters(

171

(T) 1e-4, // uMax

172

(T) 10., // Re

173

30, // N

174

1., // lx

175

1., // ly

176

1. // lz

177

);

178

const T logT = (T)1/(T)1000;

179

const T imSave = (T)1/(T)40;

180

const T vtkSave = (T)1;

181

const T maxT = (T)10.1;

182

183

pcout << "omega= " << parameters.getOmega() << std::endl;

184

writeLogFile(parameters, "3D diagonal cavity");

185

186

// @Tomasz: STEP 1

187

// Instead of simply creating a MultiBlockLattice3D as usual, the internal

188

// structure and parallelization of the block are created manually. The block

189

// is covered by 5x5x5 sub-domains. Each of these 125 sub-domains, if it

190

// contains no boundary area (i.e. only pure BGK), will then be off-loaded to

191

// the co-processor.

192

193

// Here the 5x5x5 cover-up is instantiated.

194

plint numBlocksX = 3;

195

plint numBlocksY = 3;

196

plint numBlocksZ = 3;

197

plint numBlocks = numBlocksX*numBlocksY*numBlocksZ;

198

plint envelopeWidth = 1;

199

SparseBlockStructure3D blockStructure (

200

createCavityDistribution3D (

201

parameters.getNx(), parameters.getNy(), parameters.getNz() ) );

202

203

// In case of MPI parallelism, the blocks are explicitly assigned to processors,

204

// with equal load.

205

ExplicitThreadAttribution* threadAttribution = new ExplicitThreadAttribution;

206

std::vector<std::pair<plint,plint> > ranges;

207

plint numRanges = std::min(numBlocks, (plint)global::mpi().getSize());

208

util::linearRepartition(0, numBlocks-1, numRanges, ranges);

209

for (pluint iProc=0; iProc<ranges.size(); ++iProc) {

210

for (plint blockId=ranges[iProc].first; blockId<=ranges[iProc].second; ++blockId) {

211

threadAttribution -> addBlock(blockId, iProc);

212

}

213

}

214

215

216

// Create a lattice with the above specified internal structure.

217

MultiBlockLattice3D<T, DESCRIPTOR> lattice (

218

MultiBlockManagement3D (

219

blockStructure, threadAttribution, envelopeWidth ),

220

defaultMultiBlockPolicy3D().getBlockCommunicator(),

221

defaultMultiBlockPolicy3D().getCombinedStatistics(),

222

defaultMultiBlockPolicy3D().getMultiCellAccess<T,DESCRIPTOR>(),

223

new BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()) );

224

225

226

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition

227

= createInterpBoundaryCondition3D<T,DESCRIPTOR>();

228

cavitySetup(lattice, parameters, *boundaryCondition);

229

230

// @Tomasz: STEP 2

231

// Here the co-processor is informed about the domains for which it will need to do computations.

232

bool printInfo=true;

233

initiateCoProcessors(lattice, BGKdynamics<T,DESCRIPTOR>(parameters.getOmega()).getId(), printInfo);

234

235

T previousIterationTime = T();

236

// Loop over main time iteration.

237

for (plint iT=0; iT<parameters.nStep(maxT)/1000; ++iT) {

238

global::timer("mainLoop").restart();

239

240

if (iT%parameters.nStep(imSave)==0) {

241

pcout << "Writing Gif ..." << endl;

242

writeGifs(lattice, parameters, iT);

243

}

244

245

if (iT%parameters.nStep(vtkSave)==0 && iT>0) {

246

pcout << "Saving VTK file ..." << endl;

247

writeVTK(lattice, parameters, iT);

248

}

249

250

if (iT%parameters.nStep(logT)==0) {

251

pcout << "step " << iT

252

<< "; t=" << iT*parameters.getDeltaT();

253

}

254

255

// @Tomasz: STEP 3

256

// For now, all the data assigned to co-processors is transferred from CPU memory

257

// to co-processor memory before the collision-streaming step, and back to CPU

258

// memory after the collision-streaming step. This is not efficient at all: we

259

// will replace this by communicating outer boundary layers only. It is however

260

// simpler for now, for our first proof of concept.

261

262

// Execute a time iteration.

263

transferToCoProcessors(lattice);

264

lattice.collideAndStream();

265

transferFromCoProcessors(lattice);

266

// After transferring back from co-processor to CPU memory, data in the envelopes

267

// of the CPU blocks must be synchronized.

268

lattice.duplicateOverlaps(modif::staticVariables);

269

270

271

// Access averages from internal statistics ( their value is defined

272

// only after the call to lattice.collideAndStream() )

273

if (iT%parameters.nStep(logT)==0) {

274

pcout << "; av energy="

275

<< setprecision(10) << computeAverageEnergy<T>(lattice)

276

<< "; av rho="

277

<< setprecision(10) << computeAverageDensity<T>(lattice) << endl;

278

pcout << "Time spent during previous iteration: "

279

<< previousIterationTime << endl;

280

}

281

282

previousIterationTime = global::timer("mainLoop").stop();

283

}

284

285

delete boundaryCondition;

286

}

287

288

EOF

289

290

cd ..

291

mkdir build

292

cd build

293

cmake -DCMAKE_VERBOSE_MAKEFILE=ON -DCMAKE_INSTALL_PREFIX=./../inst ./../src

294

make

295

make install

296

echo "build: OK"

297

[ -x demo ]

298

./demo

299

echo "run: OK"