LatticeModel.tmpl.h 18 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
//======================================================================================================================
//
//  This file is part of waLBerla. waLBerla is free software: you can
//  redistribute it and/or modify it under the terms of the GNU General Public
//  License as published by the Free Software Foundation, either version 3 of
//  the License, or (at your option) any later version.
//
//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
//  for more details.
//
//  You should have received a copy of the GNU General Public License along
//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \\author Martin Bauer <martin.bauer@fau.de>
//
//======================================================================================================================


#include "core/DataTypes.h"
#include "core/logging/Logging.h"

#include "field/GhostLayerField.h"
#include "field/SwapableCompare.h"
#include "domain_decomposition/BlockDataID.h"
#include "domain_decomposition/IBlock.h"
Martin Bauer's avatar
Martin Bauer committed
28
#include "stencil/{{stencil_name}}.h"
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

#include "lbm/lattice_model/EquilibriumDistribution.h"
#include "lbm/field/Density.h"
#include "lbm/field/DensityAndMomentumDensity.h"
#include "lbm/field/DensityAndVelocity.h"
#include "lbm/field/PressureTensor.h"
#include "lbm/field/ShearRate.h"

#include <set>

#ifdef __GNUC__
#define RESTRICT __restrict__
#elif _MSC_VER
#define RESTRICT __restrict
#else
#define RESTRICT
#endif

#ifdef WALBERLA_CXX_COMPILER_IS_GNU
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#endif

53
54
55
56
57
#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wunused-variable"
#pragma clang diagnostic ignored "-Wunused-parameter"
#endif
58

59
{% set lmIgnores = ('pdfs', 'pdfs_tmp') %}
60
61
62
63
64


// Forward declarations
namespace walberla{
namespace {{namespace}} {
Martin Bauer's avatar
Martin Bauer committed
65
   class {{class_name}};
66
67
68
}}
namespace walberla {
namespace mpi {
Martin Bauer's avatar
Martin Bauer committed
69
70
    mpi::SendBuffer & operator<< (mpi::SendBuffer & buf, const ::walberla::{{namespace}}::{{class_name}} & lm);
    mpi::RecvBuffer & operator>> (mpi::RecvBuffer & buf,       ::walberla::{{namespace}}::{{class_name}} & lm);
71
72
73
74
75
76
77
78
79
80
}}




namespace walberla {
namespace {{namespace}} {


/**
Martin Bauer's avatar
Martin Bauer committed
81
{{class_name}} was generated with lbmpy. Do not edit this file directly. Instead modify {{class_name}}.py.
82
83
84
85
86
87
88
89
For details see documentation of lbmpy.

Usage:
    - Create an instance of this lattice model class: the constructor parameters vary depending on the configure
      lattice model. A model with constant force needs a single force vector, while a model with variable forces needs
      a force field. All constructor parameters are ordered alphabetically.
    - Create a PDFField with the lattice model as template argument to store the particle distribution functions.
      Use the PDFField to get and modify macroscopic values.
Martin Bauer's avatar
Martin Bauer committed
90
    - The internal class {{class_name}}::Sweep is a functor to execute one LB time step.
91
92
93
      Stream, collide steps can be executed separately, or together in an optimized stream-pull-collide scheme

*/
Martin Bauer's avatar
Martin Bauer committed
94
class {{class_name}}
95
96
97
{

public:
Martin Bauer's avatar
Martin Bauer committed
98
99
    typedef stencil::{{stencil_name}} Stencil;
    typedef stencil::{{stencil_name}} CommunicationStencil;
100
101
102
    static const real_t w[{{Q}}];
    static const real_t wInv[{{Q}}];

103
    static const bool compressible = {% if compressible %}true{% else %}false{% endif %};
104
105
106
107
108
109
110
111
112

    class Sweep
    {
    public:
        Sweep( BlockDataID _pdfsID ) : pdfsID(_pdfsID) {};

        //void stream       ( IBlock * const block, const uint_t numberOfGhostLayersToInclude = uint_t(0) );
        void collide      ( IBlock * const block, const uint_t numberOfGhostLayersToInclude = uint_t(0) );
        void streamCollide( IBlock * const block, const uint_t numberOfGhostLayersToInclude = uint_t(0) );
113
        void stream       ( IBlock * const block, const uint_t numberOfGhostLayersToInclude = uint_t(0) );
114
115
116
117
118
119
120

        void operator() ( IBlock * const block, const uint_t numberOfGhostLayersToInclude = uint_t(0) )
        {
            streamCollide( block, numberOfGhostLayersToInclude );
        }

    private:
121
        {{stream_collide_kernel|generate_members(only_fields=True)|indent(8)}}
122
123
    };

Martin Bauer's avatar
Martin Bauer committed
124
125
    {{class_name}}( {{stream_collide_kernel|generate_constructor_parameters(lmIgnores) }} )
        : {{ stream_collide_kernel|generate_constructor_initializer_list(lmIgnores) }}, currentLevel(0)
126
127
128
129
    {};

    void configure( IBlock & block, StructuredBlockStorage &)  { configureBlock( &block ); }

130
131
132
    // Parameters:
    {{stream_collide_kernel|generate_members(lmIgnores)|indent(4)}}

133
134
135
private:
    void configureBlock(IBlock * block)
    {
Martin Bauer's avatar
Martin Bauer committed
136
        {{stream_collide_kernel|generate_block_data_to_field_extraction(lmIgnores, no_declarations=True)|indent(8)}}
137

Martin Bauer's avatar
Martin Bauer committed
138
        {% if refinement_scaling -%}
139
140
141
142
        const uint_t targetLevel = block->getBlockStorage().getLevel(*block);

        if( targetLevel != currentLevel )
        {
Jan Hönig's avatar
Jan Hönig committed
143
            real_t level_scale_factor(1);
144
            if( currentLevel < targetLevel )
Jan Hönig's avatar
Jan Hönig committed
145
               level_scale_factor = real_c( uint_t(1) << ( targetLevel - currentLevel ) );
146
            else // currentLevel > targetLevel
Jan Hönig's avatar
Jan Hönig committed
147
               level_scale_factor = real_t(1) / real_c( uint_t(1) << ( currentLevel - targetLevel ) );
148

Martin Bauer's avatar
Martin Bauer committed
149
            {% for scalingType, name, expression in refinement_scaling.scaling_info -%}
150
151
            {% if scalingType == 'normal' %}
            {{name}} = {{expression}};
Martin Bauer's avatar
Martin Bauer committed
152
153
            {% elif scalingType in ('field_with_f', 'field_xyz') %}
            auto it = {{name}}->{% if scalingType == 'field_with_f'%} beginWithGhostLayer(){% else %}beginWithGhostLayerXYZ(){% endif %};
154
155
156
157
158
            for( ; it != {{name}}->end(); ++it )
            {
                 auto x = it.x();
                 auto y = it.y();
                 auto z = it.z();
Martin Bauer's avatar
Martin Bauer committed
159
                 {% if scalingType == 'field_with_f' -%}
160
161
162
163
164
165
166
167
168
169
170
171
172
                 auto f = it.f();
                 {% endif -%}
                 *it = {{expression}};
            }
            {% endif -%}
            {% endfor -%}

            currentLevel = targetLevel;
        }
        {% endif -%}
    }

    // Updated by configureBlock:
Martin Bauer's avatar
Martin Bauer committed
173
    {{stream_collide_kernel|generate_block_data_to_field_extraction(lmIgnores, declarations_only=True)|indent(4)}}
174
175
176
    uint_t currentLevel;

    // Backend classes can access private members:
Martin Bauer's avatar
Martin Bauer committed
177
    friend class {{class_name}}::Sweep;
178
179
180
181
182
183
184
185
    template<class LM, class Enable> friend class  EquilibriumDistribution;
    template<class LM, class Enable> friend struct Equilibrium;
    template<class LM, class Enable> friend struct internal::AdaptVelocityToForce;
    template<class LM, class Enable> friend struct Density;
    template<class LM>               friend struct DensityAndVelocity;
    template<class LM, class Enable> friend struct DensityAndMomentumDensity;
    template<class LM, class Enable> friend struct MomentumDensity;
    template<class LM, class It, class Enable> friend struct DensityAndVelocityRange;
186

Martin Bauer's avatar
Martin Bauer committed
187
188
    friend mpi::SendBuffer & ::walberla::mpi::operator<< (mpi::SendBuffer & , const {{class_name}} & );
    friend mpi::RecvBuffer & ::walberla::mpi::operator>> (mpi::RecvBuffer & ,       {{class_name}} & );
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203

};




//======================================================================================================================
//
//  Implementation of macroscopic value backend
//
//======================================================================================================================



template<>
Martin Bauer's avatar
Martin Bauer committed
204
class EquilibriumDistribution< {{class_name}}, void>
205
{
206
public:
Martin Bauer's avatar
Martin Bauer committed
207
   typedef typename {{class_name}}::Stencil Stencil;
208
209
210
211
212

   static real_t get( const stencil::Direction direction,
                      const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ),
                      real_t rho = real_t(1.0) )
   {
213
        {% if not compressible %}
214
215
        rho -= real_t(1.0);
        {% endif %}
Martin Bauer's avatar
Martin Bauer committed
216
        {{equilibrium_from_direction}}
217
218
219
220
221
222
   }

   static real_t getSymmetricPart( const stencil::Direction direction,
                                   const Vector3<real_t> & u = Vector3< real_t >(real_t(0.0)),
                                   real_t rho = real_t(1.0) )
   {
223
        {% if not compressible %}
224
225
        rho -= real_t(1.0);
        {% endif %}
226
        {{symmetric_equilibrium_from_direction}}
227
228
229
230
231
232
   }

   static real_t getAsymmetricPart( const stencil::Direction direction,
                                    const Vector3< real_t > & u = Vector3<real_t>( real_t(0.0) ),
                                    real_t rho = real_t(1.0) )
   {
233
        {% if not compressible %}
234
235
        rho -= real_t(1.0);
        {% endif %}
Martin Bauer's avatar
Martin Bauer committed
236
        {{asymmetric_equilibrium_from_direction}}
237
238
239
240
241
   }

   static std::vector< real_t > get( const Vector3< real_t > & u = Vector3<real_t>( real_t(0.0) ),
                                     real_t rho = real_t(1.0) )
   {
242
      {% if not compressible %}
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
      rho -= real_t(1.0);
      {% endif %}

      std::vector< real_t > equilibrium( Stencil::Size );
      for( auto d = Stencil::begin(); d != Stencil::end(); ++d )
      {
         equilibrium[d.toIdx()] = get(*d, u, rho);
      }
      return equilibrium;
   }
};


namespace internal {

template<>
Martin Bauer's avatar
Martin Bauer committed
259
struct AdaptVelocityToForce<{{class_name}}, void>
260
261
{
   template< typename FieldPtrOrIterator >
Martin Bauer's avatar
Martin Bauer committed
262
   static Vector3<real_t> get( FieldPtrOrIterator & it, const {{class_name}} & lm,
263
264
265
266
267
                               const Vector3< real_t > & velocity, const real_t rho )
   {
      auto x = it.x();
      auto y = it.y();
      auto z = it.z();
Martin Bauer's avatar
Martin Bauer committed
268
269
      {% if macroscopic_velocity_shift %}
      return velocity - Vector3<real_t>({{macroscopic_velocity_shift | join(",") }} {% if D == 2 %}, 0.0 {%endif %} );
270
271
272
273
274
      {% else %}
      return velocity;
      {% endif %}
   }

Martin Bauer's avatar
Martin Bauer committed
275
   static Vector3<real_t> get( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const {{class_name}} & lm,
276
277
                               const Vector3< real_t > & velocity, const real_t rho )
   {
Martin Bauer's avatar
Martin Bauer committed
278
      {% if macroscopic_velocity_shift %}
279

Martin Bauer's avatar
Martin Bauer committed
280
      return velocity - Vector3<real_t>({{macroscopic_velocity_shift | join(",") }} {% if D == 2 %}, 0.0 {%endif %} );
281
282
283
284
285
286
287
288
289
290
      {% else %}
      return velocity;
      {% endif %}
   }
};
} // namespace internal



template<>
Martin Bauer's avatar
Martin Bauer committed
291
struct Equilibrium< {{class_name}}, void >
292
293
294
295
{

   template< typename FieldPtrOrIterator >
   static void set( FieldPtrOrIterator & it,
296
                    const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), real_t rho = real_t(1.0) )
297
   {
298
        {%if not compressible %}
299
300
301
        rho -= real_t(1.0);
        {%endif %}

302
       {% for eqTerm in equilibrium -%}
303
       it[{{loop.index0 }}] = {{eqTerm}};
304
305
306
307
308
       {% endfor -%}
   }

   template< typename PdfField_T >
   static void set( PdfField_T & pdf, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z,
309
                    const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), real_t rho = real_t(1.0) )
310
   {
311
      {%if not compressible %}
312
313
314
      rho -= real_t(1.0);
      {%endif %}

315
316
      real_t & xyz0 = pdf(x,y,z,0);
      {% for eqTerm in equilibrium -%}
317
      pdf.getF( &xyz0, {{loop.index0 }})= {{eqTerm}};
318
319
320
321
322
323
      {% endfor -%}
   }
};


template<>
Martin Bauer's avatar
Martin Bauer committed
324
struct Density<{{class_name}}, void>
325
326
{
   template< typename FieldPtrOrIterator >
Martin Bauer's avatar
Martin Bauer committed
327
   static inline real_t get( const {{class_name}} & , const FieldPtrOrIterator & it )
328
329
330
331
   {
        {% for i in range(Q) -%}
            const real_t f_{{i}} = it[{{i}}];
        {% endfor -%}
Martin Bauer's avatar
Martin Bauer committed
332
        {{density_getters | indent(8)}}
333
334
335
336
        return rho;
   }

   template< typename PdfField_T >
Martin Bauer's avatar
Martin Bauer committed
337
   static inline real_t get( const {{class_name}} & ,
338
339
340
341
342
343
                             const PdfField_T & pdf, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
   {
        const real_t & xyz0 = pdf(x,y,z,0);
        {% for i in range(Q) -%}
            const real_t f_{{i}} = pdf.getF( &xyz0, {{i}});
        {% endfor -%}
Martin Bauer's avatar
Martin Bauer committed
344
        {{density_getters | indent(8)}}
345
346
347
348
349
350
        return rho;
   }
};


template<>
Martin Bauer's avatar
Martin Bauer committed
351
struct DensityAndVelocity<{{class_name}}>
352
353
{
    template< typename FieldPtrOrIterator >
Martin Bauer's avatar
Martin Bauer committed
354
    static void set( FieldPtrOrIterator & it, const {{class_name}} & lm,
Jan Hönig's avatar
Jan Hönig committed
355
                     const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), const real_t rho_in = real_t(1.0) )
356
357
358
359
360
    {
        auto x = it.x();
        auto y = it.y();
        auto z = it.z();

Martin Bauer's avatar
Martin Bauer committed
361
        {{density_velocity_setter_macroscopic_values | indent(8)}}
362
363
364
365
        {% if D == 2 -%}
        const real_t u_2(0.0);
        {% endif %}

366
        Equilibrium<{{class_name}}>::set(it, Vector3<real_t>(u_0, u_1, u_2), rho{%if not compressible %} + real_t(1) {%endif%});
367
368
369
    }

    template< typename PdfField_T >
Martin Bauer's avatar
Martin Bauer committed
370
    static void set( PdfField_T & pdf, const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const {{class_name}} & lm,
Jan Hönig's avatar
Jan Hönig committed
371
                     const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), const real_t rho_in = real_t(1.0) )
372
    {
Martin Bauer's avatar
Martin Bauer committed
373
        {{density_velocity_setter_macroscopic_values | indent(8)}}
374
375
376
377
        {% if D == 2 -%}
        const real_t u_2(0.0);
        {% endif %}

378
        Equilibrium<{{class_name}}>::set(pdf, x, y, z, Vector3<real_t>(u_0, u_1, u_2), rho {%if not compressible %} + real_t(1) {%endif%});
379
380
381
382
383
    }
};


template<typename FieldIteratorXYZ >
Martin Bauer's avatar
Martin Bauer committed
384
struct DensityAndVelocityRange<{{class_name}}, FieldIteratorXYZ>
385
386
{

Martin Bauer's avatar
Martin Bauer committed
387
   static void set( FieldIteratorXYZ & begin, const FieldIteratorXYZ & end, const {{class_name}} & lm,
Jan Hönig's avatar
Jan Hönig committed
388
                    const Vector3< real_t > & u = Vector3< real_t >( real_t(0.0) ), const real_t rho_in = real_t(1.0) )
389
390
391
392
393
394
   {
        for( auto cellIt = begin; cellIt != end; ++cellIt )
        {
            const auto x = cellIt.x();
            const auto y = cellIt.y();
            const auto z = cellIt.z();
Martin Bauer's avatar
Martin Bauer committed
395
            {{density_velocity_setter_macroscopic_values | indent(12)}}
396
397
398
399
            {% if D == 2 -%}
            const real_t u_2(0.0);
            {% endif %}

400
            Equilibrium<{{class_name}}>::set(cellIt, Vector3<real_t>(u_0, u_1, u_2), rho{%if not compressible %} + real_t(1) {%endif%});
401
402
403
404
405
406
407
        }
   }
};



template<>
Martin Bauer's avatar
Martin Bauer committed
408
struct DensityAndMomentumDensity<{{class_name}}>
409
410
{
   template< typename FieldPtrOrIterator >
Martin Bauer's avatar
Martin Bauer committed
411
   static real_t get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm,
412
413
414
415
416
417
418
419
420
421
                      const FieldPtrOrIterator & it )
   {
        const auto x = it.x();
        const auto y = it.y();
        const auto z = it.z();

        {% for i in range(Q) -%}
            const real_t f_{{i}} = it[{{i}}];
        {% endfor -%}

Martin Bauer's avatar
Martin Bauer committed
422
        {{momentum_density_getter | indent(8) }}
423
424
425
426
427
428
429
        {% for i in range(D) -%}
            momentumDensity[{{i}}] = md_{{i}};
        {% endfor %}
        return rho;
   }

   template< typename PdfField_T >
Martin Bauer's avatar
Martin Bauer committed
430
   static real_t get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm, const PdfField_T & pdf,
431
432
433
434
435
436
437
                      const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
   {
        const real_t & xyz0 = pdf(x,y,z,0);
        {% for i in range(Q) -%}
            const real_t f_{{i}} = pdf.getF( &xyz0, {{i}});
        {% endfor -%}

Martin Bauer's avatar
Martin Bauer committed
438
        {{momentum_density_getter | indent(8) }}
439
440
441
442
443
444
445
446
447
        {% for i in range(D) -%}
            momentumDensity[{{i}}] = md_{{i}};
        {% endfor %}
       return rho;
   }
};


template<>
Martin Bauer's avatar
Martin Bauer committed
448
struct MomentumDensity< {{class_name}}>
449
450
{
   template< typename FieldPtrOrIterator >
Martin Bauer's avatar
Martin Bauer committed
451
   static void get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm, const FieldPtrOrIterator & it )
452
453
454
455
456
457
458
459
460
   {
        const auto x = it.x();
        const auto y = it.y();
        const auto z = it.z();

        {% for i in range(Q) -%}
            const real_t f_{{i}} = it[{{i}}];
        {% endfor -%}

Martin Bauer's avatar
Martin Bauer committed
461
        {{momentum_density_getter | indent(8) }}
462
463
464
465
466
467
        {% for i in range(D) -%}
            momentumDensity[{{i}}] = md_{{i}};
        {% endfor %}
   }

   template< typename PdfField_T >
Martin Bauer's avatar
Martin Bauer committed
468
   static void get( Vector3< real_t > & momentumDensity, const {{class_name}} & lm, const PdfField_T & pdf,
469
470
471
472
473
474
475
                    const cell_idx_t x, const cell_idx_t y, const cell_idx_t z )
   {
        const real_t & xyz0 = pdf(x,y,z,0);
        {% for i in range(Q) -%}
            const real_t f_{{i}} = pdf.getF( &xyz0, {{i}});
        {% endfor -%}

Martin Bauer's avatar
Martin Bauer committed
476
        {{momentum_density_getter | indent(8) }}
477
478
479
480
481
482
483
484
        {% for i in range(D) -%}
            momentumDensity[{{i}}] = md_{{i}};
        {% endfor %}
   }
};


template<>
Martin Bauer's avatar
Martin Bauer committed
485
struct PressureTensor<{{class_name}}>
486
487
{
   template< typename FieldPtrOrIterator >
Martin Bauer's avatar
Martin Bauer committed
488
   static void get( Matrix3< real_t > & /* pressureTensor */, const {{class_name}} & /* latticeModel */, const FieldPtrOrIterator & /* it */ )
489
490
491
492
493
   {
       WALBERLA_ABORT("Not implemented");
   }

   template< typename PdfField_T >
Martin Bauer's avatar
Martin Bauer committed
494
   static void get( Matrix3< real_t > & /* pressureTensor */, const {{class_name}} & /* latticeModel */, const PdfField_T & /* pdf */,
495
                    const cell_idx_t /* x */, const cell_idx_t /* y */, const cell_idx_t /* z */ )
496
497
498
499
500
501
502
   {
       WALBERLA_ABORT("Not implemented");
   }
};


template<>
Martin Bauer's avatar
Martin Bauer committed
503
struct ShearRate<{{class_name}}>
504
505
{
   template< typename FieldPtrOrIterator >
Martin Bauer's avatar
Martin Bauer committed
506
   static inline real_t get( const {{class_name}} & /* latticeModel */, const FieldPtrOrIterator & /* it */,
507
                             const Vector3< real_t > & /* velocity */, const real_t /* rho */)
508
509
510
511
512
513
   {
       WALBERLA_ABORT("Not implemented");
       return real_t(0.0);
   }

   template< typename PdfField_T >
Martin Bauer's avatar
Martin Bauer committed
514
   static inline real_t get( const {{class_name}} & latticeModel,
515
516
                             const PdfField_T & /* pdf */, const cell_idx_t /* x */, const cell_idx_t /* y */, const cell_idx_t /* z */,
                             const Vector3< real_t > & /* velocity */, const real_t /* rho */ )
517
518
519
520
521
   {
       WALBERLA_ABORT("Not implemented");
       return real_t(0.0);
   }

522
523
   static inline real_t get( const std::vector< real_t > & /* nonEquilibrium */, const real_t /* relaxationParam */,
                             const real_t /* rho */ = real_t(1) )
524
525
526
527
528
529
530
531
532
533
534
535
536
537
   {
       WALBERLA_ABORT("Not implemented");
       return real_t(0.0);
   }
};


} // namespace {{namespace}}
} // namespace walberla



#ifdef WALBERLA_CXX_COMPILER_IS_GNU
#pragma GCC diagnostic pop
538
539
540
541
542
#endif

#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
#pragma clang diagnostic pop
#endif