Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 121-buffersystem-receiver-info-without-sender-ranks
  • 3-stable
  • 4-stable
  • AddaptTypeSystem
  • CMakeCodeGenPartTwo
  • ChannelFlow
  • CoVortex
  • CodegenForRefinement
  • CommunicationGPUBenchmark
  • ComnbinedGPUPackinfo
  • ExportCudaDeviceSelection
  • FixSinglePrecisionProblems
  • FlagFieldExample
  • FlowAroundSphere
  • FreeSurface
  • GPURefineTest
  • GPURefinement
  • GPURefinementImprovement
  • HRR
  • HydroPressure
  • IBC
  • InterpolationBC
  • Italy
  • LDC
  • Lagoon
  • LeesEdwards
  • ListLBM
  • NewChannelBenchmark
  • Remove_fSize_from_templates
  • SphereMovie
  • TGA
  • TaylorBubble
  • TurbulentChannel
  • UpgradePystencils
  • VTKUnstructured
  • clang11
  • develop
  • develop2
  • fluidizedbed_showcase
  • master
  • phaseField
  • phasefield-drop
  • porous
  • porousHeat
  • remiPorous
  • s2a
  • setup_walberla_codegen
  • vbondmodel_integrated
  • vbondmodel_isotropic
  • v3.1
  • v3.2
  • v3.3
  • v4.0dev
  • v4.1
  • v4.2
  • v5.0dev
56 results

Target

Select target project
  • castellsc/walberla
  • ravi.k.ayyala/walberla
  • em73etav/walberla
  • hoenig/walberla
  • le45zyci/walberla
  • sudesh.rathnayake/walberla
  • el38efib/walberla
  • rahil.doshi/walberla
  • Bindgen/walberla
  • ob28imeq/walberla
  • ArashPartow/walberla
  • jarmatz/walberla
  • jbadwaik/walberla
  • ec93ujoh/walberla
  • walberla/walberla
  • ProjectPhysX/walberla
  • shellshocked2003/walberla
  • stewart/walberla
  • behzad.safaei/walberla
  • schruff/walberla
  • loreson/walberla
  • Novermars/walberla
  • itischler/walberla
  • holzer/walberla
  • da15siwa/walberla
  • he66coqe/walberla
  • jngrad/walberla
  • uq60ifih/walberla
  • ostanin/walberla
  • bauer/walberla
  • zy79zopo/walberla
  • jonas_schmitt/walberla
  • po60nani/walberla
  • ro36vugi/walberla
  • fweik/walberla
  • ab04unyc/walberla
  • yw25ynew/walberla
  • ig38otak/walberla
  • RudolfWeeber/walberla
39 results
Select Git revision
  • 121-buffersystem-receiver-info-without-sender-ranks
  • 128-some-tests-are-not-active
  • 146-cuda-gcc-config-warning
  • 3-stable
  • 4-stable
  • 5-stable
  • 6-stable
  • 7-stable
  • 727-refactor-sqlExport
  • AtomicAdd_for_CUDA_compute_capabilities<6.0
  • ChargedParticles
  • CodegenForRefinement
  • GeneratedOutflowBC
  • RayleighBernardConvection
  • Remove_fSize_from_templates
  • UpdateGPUBenchmark
  • UpdatePhaseField
  • angersbach/coding-day-01-09
  • antidunes-visualization
  • bam_piping_erosion
  • benchmark_sqlite_modify
  • change-default-layout-fzyx
  • clang-tidy
  • clang11
  • clang_tidy2
  • cmake_cleanup
  • cnt_app
  • codegen-update
  • coding-day-01-09-mesh
  • coupling_tutorial
  • doshi/coding-day-01-09
  • externalize_dependencies
  • fix_nvcc_compiler_warnings
  • fluidizedbed_showcase
  • hip-ShiftedPeriodicity
  • kajol/coding-day
  • kemmler/particle_coupling_GPU
  • lbmpy-kernel-comparison
  • master
  • plewinski/fix-Guo-force-model-TRT-MRT
  • pystencils2.0-adoption
  • rangersbach/doxygen_style
  • ravi/coding-day
  • ravi/material_transport
  • setup_walberla_codegen
  • suction_bucket
  • suffa/NorthWind
  • suffa/NorthWind_refined
  • suffa/SYCL
  • suffa/Sparse
  • suffa/compact_interpolation
  • suffa/fix_force_on_boundary
  • suffa/integrate_moving_geo
  • suffa/psm_lbm_package
  • thermalFreeSurfaceLBM
  • thoennes/cusotm-mpi-reduce-function
  • use-correct-codegen-data-type
  • viscLDCwithFSLBM
  • v3.1
  • v3.2
  • v3.3
  • v4.0dev
  • v4.1
  • v4.2
  • v5.0dev
  • v5.1
  • v6.0dev
  • v6.1
  • v7.0dev
  • v7.1
70 results
Show changes
Showing
with 4324 additions and 0 deletions
/Users/holzer/walberla/cmake-build-release/apps/showcases/Channel/Channel Channel.prm
[0][INFO ]------(0.000 sec) Creating the block structure ...
[0][INFO ]------(0.000 sec) SetupBlockForest created successfully:
[0] - AABB: [ <0,0,0>, <8,8,4> ]
[0] - initial decomposition: 2 x 2 x 1 (= forest size)
[0] - periodicity: true x false x true
[0] - number of blocks discarded from the initial grid: 0 (= 0 %)
[0] - number of levels: 2
[0] - tree ID digits: 3 (-> block ID bytes = 1)
[0] - total number of blocks: 18
[0] - number of processes: 1 (1 worker process(es) / 0 empty buffer process(es))
[0] - buffer processes are inserted into the process network: no
[0] - process ID bytes: 0
[0] - blocks/memory/workload per process:
[0] + blocks:
[0] - min = 18
[0] - max = 18
[0] - avg = 18
[0] - stdDev = 0
[0] - relStdDev = 0
[0] + memory:
[0] - min = 4.35917e+06
[0] - max = 4.35917e+06
[0] - avg = 4.35917e+06
[0] - stdDev = 0
[0] - relStdDev = 0
[0] + workload:
[0] - min = 34
[0] - max = 34
[0] - avg = 34
[0] - stdDev = 0
[0] - relStdDev = 0
[0] - distribution of space/memory/work to different grid levels:
[0] + level 0:
[0] - 2 blocks ...
[0] - ... cover 50 % of the total simulation space
[0] - ... account for 11.1111 % of the total memory foot print
[0] - ... generate 5.88235 % of the total work load
[0] - blocks per process:
[0] + min = 2
[0] + max = 2
[0] + avg = 2
[0] + stdDev = 0
[0] + relStdDev = 0
[0] - memory per process:
[0] + min = 484352
[0] + max = 484352
[0] + avg = 484352
[0] + stdDev = 0
[0] + relStdDev = 0
[0] - workload per process:
[0] + min = 2
[0] + max = 2
[0] + avg = 2
[0] + stdDev = 0
[0] + relStdDev = 0
[0] + level 1:
[0] - 16 blocks ...
[0] - ... cover 50 % of the total simulation space
[0] - ... account for 88.8889 % of the total memory foot print
[0] - ... generate 94.1176 % of the total work load
[0] - blocks per process:
[0] + min = 16
[0] + max = 16
[0] + avg = 16
[0] + stdDev = 0
[0] + relStdDev = 0
[0] - memory per process:
[0] + min = 3.87482e+06
[0] + max = 3.87482e+06
[0] + avg = 3.87482e+06
[0] + stdDev = 0
[0] + relStdDev = 0
[0] - workload per process:
[0] + min = 32
[0] + max = 32
[0] + avg = 32
[0] + stdDev = 0
[0] + relStdDev = 0
[0] - using a uniform decomposition with a resolution equal to the finest level, one would ...
[0] + ... need 1.77778 times the memory
[0] + ... generate 1.88235 times the workload
[0][INFO ]------(0.003 sec) Setting up communication...
[0][INFO ]------(0.008 sec) Blocks created: 18
[0][INFO ]------(0.008 sec) Level 0 Blocks: 2
[0][INFO ]------(0.008 sec) Level 1 Blocks: 16
[0][INFO ]------(0.008 sec) Starting Simulation
[0][INFO ]------(0.013 sec)[ 0] Evaluation of accuracy:
[0] - L1: 0.00479167
[0] - L2: 0.005469
[0] - Lmax: 0.0090625
[0][INFO ]------(0.435 sec)[ 100] Evaluation of accuracy:
[0] - L1: 0.00134704
[0] - L2: 0.00149049
[0] - Lmax: 0.002099
[0][INFO ]------(0.851 sec)[ 200] Evaluation of accuracy:
[0] - L1: 0.000448633
[0] - L2: 0.000496303
[0] - Lmax: 0.000696362
[0][INFO ]------(1.262 sec)[ 300] Evaluation of accuracy:
[0] - L1: 0.000149393
[0] - L2: 0.000165266
[0] - Lmax: 0.00023185
[0][INFO ]------(1.673 sec)[ 400] Evaluation of accuracy:
[0] - L1: 4.97466e-05
[0] - L2: 5.50323e-05
[0] - Lmax: 7.7204e-05
[0][INFO ]------(2.085 sec)[ 500] Evaluation of accuracy:
[0] - L1: 1.65652e-05
[0] - L2: 1.83253e-05
[0] - Lmax: 2.57083e-05
[0][INFO ]------(2.507 sec)[ 600] Evaluation of accuracy:
[0] - L1: 5.5161e-06
[0] - L2: 6.1022e-06
[0] - Lmax: 8.56068e-06
[0][INFO ]------(2.919 sec)[ 700] Evaluation of accuracy:
[0] - L1: 1.83682e-06
[0] - L2: 2.03199e-06
[0] - Lmax: 2.85064e-06
[0][INFO ]------(3.330 sec)[ 800] Evaluation of accuracy:
[0] - L1: 6.11647e-07
[0] - L2: 6.76636e-07
[0] - Lmax: 9.49242e-07
[0][INFO ]------(3.741 sec)[ 900] Evaluation of accuracy:
[0] - L1: 2.03674e-07
[0] - L2: 2.25315e-07
[0] - Lmax: 3.1609e-07
[0][INFO ]------(4.153 sec)[1000] Evaluation of accuracy:
[0] - L1: 6.78218e-08
[0] - L2: 7.5028e-08
[0] - Lmax: 1.05256e-07
[0][INFO ]------(4.162 sec) Simulation finished
[0][RESULT ]------(4.162 sec) Simulation performance:
[0] - processes: 1
[0] - threads: 1 (threads per process = 1, threads per core = n/a *))
[0] - cores: n/a *)
[0] - time steps: 1001 (on the coarsest grid, 2002 on the finest grid)
[0] - time: 4.1544 sec
[0] - cells: 1152 (2048 if everything were fine -> data reduction by factor of 1.77778)
[0] - fluid cells: 1152 (100 % of all cells)
[0] - distribution of cells to different grid levels:
[0] + level 0: 128 cells (128 fluid cells = 100 % of all cells on this level)
[0] + level 1: 1024 cells (1024 fluid cells = 100 % of all cells on this level)
[0] - performance: 0.524305 MLUPS (million lattice cell updates per second)
[0] 0.524305 MLUPS / process
[0] n/a *) MLUPS / core
[0] 0.524305 MFLUPS (million fluid lattice cell updates per second)
[0] 0.524305 MFLUPS / process
[0] n/a *) MFLUPS / core
[0] 240.949 time steps / second
[0] - 'virtual' performance (if everything were fine): 0.986928 MLUPS (million lattice cell updates per second)
[0] 0.986928 MLUPS / process
[0] n/a *) MLUPS / core
[0] 0.986928 MFLUPS (million fluid lattice cell updates per second)
[0] 0.986928 MFLUPS / process
[0] n/a *) MFLUPS / core
[0] 481.898 fine time steps / second
[0] - build / run information:
[0] + host machine: gazleu.cerfacs.fr
[0] + build machine: gazleu.cerfacs.fr
[0] + git SHA1: d5dba67210e6cc8141af6ce4b87c70822dca19e9
[0] + build type: Release
[0] + compiler flags: -O3 -DNDEBUG -Wall -Wconversion -Wshadow -Wno-c++11-extensions -Qunused-arguments -ffast-math -pthread
[0]
[0] *) only available if environment variable 'THREADS_PER_CORE' is set
[0][RESULT ]------(4.162 sec) Time loop timing:
[0] Timer | % | Total| Average| Count| Min| Max| Variance|
[0] -----------------------------------------------------------------------------------------------------------------------------------
[0] Refinement Cycle | 99.24% | 4.12| 0.004| 1001| 0.004| 0.005| 0.000|
[0] VTK Output | 0.55% | 0.02| 0.000| 1001| 0.000| 0.009| 0.000|
[0] evaluation | 0.20% | 0.01| 0.000| 1001| 0.000| 0.001| 0.000|
[0] remaining time logger | 0.00% | 0.00| 0.000| 1001| 0.000| 0.000| 0.000|
[0]
[0][WARNING ]------(4.162 sec) Getting memory statistics is currently not supported on non-Linux systems!
[0][INFO ]------(4.162 sec) resident memory: Sample has 1 values in [0.000000, 0.000000], sum = 0.000000, mean = 0.000000, med = 0.000000, stddev = 0.000000 (relative: nan), mad = 0.000000
Process finished with exit code 0
\ No newline at end of file
waLBerla_link_files_to_builddir( *.py )
waLBerla_link_files_to_builddir( *.prm )
waLBerla_link_files_to_builddir( *.png )
waLBerla_link_files_to_builddir( *.obj )
waLBerla_link_files_to_builddir( *.stl )
waLBerla_link_files_to_builddir( *.mtl )
waLBerla_generate_target_from_python(NAME FlowAroundCylinderCodeGen
FILE FlowAroundCylinder.py
OUT_FILES FlowAroundCylinderSweepCollection.h FlowAroundCylinderSweepCollection.${CODEGEN_FILE_SUFFIX}
FlowAroundCylinderStorageSpecification.h FlowAroundCylinderStorageSpecification.${CODEGEN_FILE_SUFFIX}
FreeSlip.h FreeSlip.${CODEGEN_FILE_SUFFIX}
NoSlip.h NoSlip.${CODEGEN_FILE_SUFFIX}
Obstacle.h Obstacle.${CODEGEN_FILE_SUFFIX}
UBB.h UBB.${CODEGEN_FILE_SUFFIX}
Outflow.h Outflow.${CODEGEN_FILE_SUFFIX}
FlowAroundCylinderBoundaryCollection.h
FlowAroundCylinderInfoHeader.h
FlowAroundCylinderStaticDefines.h)
if (WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
waLBerla_add_executable ( NAME FlowAroundCylinder
FILES FlowAroundCylinder.cpp Cylinder.cpp Evaluation.cpp
DEPENDS blockforest boundary core field gpu lbm_generated stencil timeloop vtk FlowAroundCylinderCodeGen )
else()
waLBerla_add_executable ( NAME FlowAroundCylinder
FILES FlowAroundCylinder.cpp Cylinder.cpp Evaluation.cpp
DEPENDS blockforest boundary core field lbm_generated stencil timeloop vtk FlowAroundCylinderCodeGen )
endif(WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
//======================================================================================================================
//
// 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/>.
//
//! \file Cylinder.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "Cylinder.h"
namespace walberla
{
bool Cylinder::contains(const Vector3< real_t >& point) const
{
const real_t px = setup_.cylinderXPosition;
const real_t py = setup_.cylinderYPosition;
const real_t H = setup_.H;
const real_t r = setup_.cylinderRadius;
if (setup_.circularCrossSection)
{
const real_t xd = point[0] - px;
const real_t yd = point[1] - py;
const real_t d = xd * xd + yd * yd;
return point[2] >= real_t(0.0) && point[2] <= H && d <= (r * r);
}
else
{
const AABB cylinder(px - r, py - r, real_t(0), px + r, py + r, H);
return cylinder.contains(point);
}
}
bool Cylinder::contains(const AABB& aabb) const
{
if (setup_.circularCrossSection)
{
Vector3< real_t > p[8];
p[0].set(aabb.xMin(), aabb.yMin(), aabb.zMin());
p[1].set(aabb.xMax(), aabb.yMin(), aabb.zMin());
p[2].set(aabb.xMin(), aabb.yMax(), aabb.zMin());
p[3].set(aabb.xMax(), aabb.yMax(), aabb.zMin());
p[4].set(aabb.xMin(), aabb.yMin(), aabb.zMax());
p[5].set(aabb.xMax(), aabb.yMin(), aabb.zMax());
p[6].set(aabb.xMin(), aabb.yMax(), aabb.zMax());
p[7].set(aabb.xMax(), aabb.yMax(), aabb.zMax());
return contains(p[0]) && contains(p[1]) && contains(p[2]) && contains(p[3]) && contains(p[4]) && contains(p[5]) &&
contains(p[6]) && contains(p[7]);
}
else { return contains(aabb.min()) && contains(aabb.max()); }
}
bool Cylinder::intersects(const AABB& aabb, const real_t bufferDistance) const
{
const real_t px = setup_.cylinderXPosition;
const real_t py = setup_.cylinderYPosition;
const real_t r = setup_.cylinderRadius;
if (setup_.circularCrossSection)
{
Vector3< real_t > p(px, py, real_t(0));
if (p[0] < aabb.xMin())
p[0] = aabb.xMin();
else if (p[0] > aabb.xMax())
p[0] = aabb.xMax();
if (p[1] < aabb.yMin())
p[1] = aabb.yMin();
else if (p[1] > aabb.yMax())
p[1] = aabb.yMax();
const real_t xd = p[0] - px;
const real_t yd = p[1] - py;
const real_t d = xd * xd + yd * yd;
const real_t rr = (r + bufferDistance) * (r + bufferDistance);
return d <= rr;
}
else
{
const AABB cylinder(px - r, py - r, real_t(0), px + r, py + r, setup_.H);
return cylinder.intersects(aabb, bufferDistance);
}
}
real_t Cylinder::delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const
{
WALBERLA_CHECK(!contains(fluid), "Point ( " << fluid << " ) is not outside the cylinder!")
WALBERLA_CHECK(contains(boundary), "Point ( " << boundary << " ) is not inside the cylinder!")
const real_t px = setup_.cylinderXPosition;
const real_t py = setup_.cylinderYPosition;
const real_t r = setup_.cylinderRadius;
if (setup_.circularCrossSection)
{
// http://devmag.org.za/2009/04/17/basic-collision-detection-in-2d-part-2/
const Vector3< real_t > circle(px, py, real_t(0));
const Vector3< real_t > f = fluid - circle;
const Vector3< real_t > d = (boundary - circle) - f;
const real_t a = d[0] * d[0] + d[1] * d[1];
const real_t b = real_t(2) * (d[0] * f[0] + d[1] * f[1]);
const real_t c = f[0] * f[0] + f[1] * f[1] - r * r;
const real_t bb4ac = b * b - (real_t(4) * a * c);
WALBERLA_CHECK_GREATER_EQUAL(bb4ac, real_t(0))
const real_t sqrtbb4ac = std::sqrt(bb4ac);
const real_t alpha = std::min((-b + sqrtbb4ac) / (real_t(2) * a), (-b - sqrtbb4ac) / (real_t(2) * a));
WALBERLA_CHECK_GREATER_EQUAL(alpha, real_t(0))
WALBERLA_CHECK_LESS_EQUAL(alpha, real_t(1))
return alpha;
}
const AABB cylinder(px - r, py - r, real_t(0), px + r, py + r, setup_.H);
if (fluid[0] <= cylinder.xMin())
{
const real_t xdiff = cylinder.xMin() - fluid[0];
if (fluid[1] <= cylinder.yMin())
{
const real_t ydiff = cylinder.yMin() - fluid[1];
if (xdiff >= ydiff)
{
WALBERLA_CHECK_LESS_EQUAL(fluid[0], boundary[0])
return xdiff / (boundary[0] - fluid[0]);
}
WALBERLA_CHECK_LESS_EQUAL(fluid[1], boundary[1])
return ydiff / (boundary[1] - fluid[1]);
}
else if (fluid[1] >= cylinder.yMax())
{
const real_t ydiff = fluid[1] - cylinder.yMax();
if (xdiff >= ydiff)
{
WALBERLA_CHECK_LESS_EQUAL(fluid[0], boundary[0])
return xdiff / (boundary[0] - fluid[0]);
}
WALBERLA_CHECK_GREATER_EQUAL(fluid[1], boundary[1])
return ydiff / (fluid[1] - boundary[1]);
}
WALBERLA_CHECK_LESS_EQUAL(fluid[0], boundary[0])
return xdiff / (boundary[0] - fluid[0]);
}
else if (fluid[0] >= cylinder.xMax())
{
const real_t xdiff = fluid[0] - cylinder.xMax();
if (fluid[1] <= cylinder.yMin())
{
const real_t ydiff = cylinder.yMin() - fluid[1];
if (xdiff >= ydiff)
{
WALBERLA_CHECK_GREATER_EQUAL(fluid[0], boundary[0])
return xdiff / (fluid[0] - boundary[0]);
}
WALBERLA_CHECK_LESS_EQUAL(fluid[1], boundary[1])
return ydiff / (boundary[1] - fluid[1]);
}
else if (fluid[1] >= cylinder.yMax())
{
const real_t ydiff = fluid[1] - cylinder.yMax();
if (xdiff >= ydiff)
{
WALBERLA_CHECK_GREATER_EQUAL(fluid[0], boundary[0])
return xdiff / (fluid[0] - boundary[0]);
}
WALBERLA_CHECK_GREATER_EQUAL(fluid[1], boundary[1])
return ydiff / (fluid[1] - boundary[1]);
}
WALBERLA_CHECK_GREATER_EQUAL(fluid[0], boundary[0])
return xdiff / (fluid[0] - boundary[0]);
}
if (fluid[1] <= cylinder.yMin())
{
WALBERLA_CHECK_LESS_EQUAL(fluid[1], boundary[1])
return (cylinder.yMin() - fluid[1]) / (boundary[1] - fluid[1]);
}
WALBERLA_CHECK_GREATER_EQUAL(fluid[1], cylinder.yMax())
WALBERLA_CHECK_GREATER_EQUAL(fluid[1], boundary[1])
return (fluid[1] - cylinder.yMax()) / (fluid[1] - boundary[1]);
}
real_t wallDistance::operator()(const Cell& fluidCell, const Cell& boundaryCell,
const shared_ptr< StructuredBlockForest >& SbF, IBlock& block) const
{
Vector3< real_t > boundary = SbF->getBlockLocalCellCenter( block, boundaryCell );
Vector3< real_t > fluid = SbF->getBlockLocalCellCenter( block, fluidCell );
SbF->mapToPeriodicDomain(boundary);
SbF->mapToPeriodicDomain(fluid);
WALBERLA_ASSERT(cylinder_.contains(boundary), "Boundary cell must be inside the cylinder!")
WALBERLA_ASSERT(!cylinder_.contains(fluid), "Fluid cell must not be inside the cylinder!")
return cylinder_.delta( fluid, boundary );
}
} // namespace walberla
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file Cylinder.h
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#pragma once
#include "blockforest/SetupBlock.h"
#include "blockforest/SetupBlockForest.h"
#include "blockforest/StructuredBlockForest.h"
#include "domain_decomposition/IBlock.h"
#include "core/DataTypes.h"
#include "core/math/AABB.h"
#include "core/math/Vector3.h"
#include "core/cell/Cell.h"
#include "Setup.h"
namespace walberla
{
class Cylinder
{
public:
Cylinder(const Setup& setup) : setup_(setup) {}
bool operator()(const Vector3< real_t >& point) const { return contains(point); }
bool contains(const Vector3< real_t >& point) const;
bool contains(const AABB& aabb) const;
bool intersects(const AABB& aabb, const real_t bufferDistance = real_t(0)) const;
real_t delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const;
private:
Setup setup_;
}; // class Cylinder
class CylinderRefinementSelection
{
public:
CylinderRefinementSelection(const Cylinder& cylinder, const uint_t level, const real_t bufferDistance)
: cylinder_(cylinder), level_(level), bufferDistance_(bufferDistance)
{}
void operator()(SetupBlockForest& forest)
{
for (auto block = forest.begin(); block != forest.end(); ++block)
{
const AABB& aabb = block->getAABB();
if (block->getLevel() < level_ && cylinder_.intersects(aabb, bufferDistance_) && !cylinder_.contains(aabb))
block->setMarker(true);
}
}
private:
Cylinder cylinder_;
uint_t level_;
real_t bufferDistance_;
}; // class CylinderRefinementSelection
class CylinderBlockExclusion
{
public:
CylinderBlockExclusion(const Cylinder& cylinder) : cylinder_(cylinder) {}
bool operator()(const blockforest::SetupBlock& block)
{
const AABB aabb = block.getAABB();
return static_cast< bool >(cylinder_.contains(aabb));
}
private:
Cylinder cylinder_;
}; // class CylinderBlockExclusion
class wallDistance
{
public:
wallDistance(const Cylinder& cylinder) : cylinder_(cylinder) {}
real_t operator()(const Cell& fluidCell, const Cell& boundaryCell, const shared_ptr< StructuredBlockForest >& SbF,
IBlock& block) const;
private:
Cylinder cylinder_;
}; // class wallDistance
} // namespace walberla
\ No newline at end of file
Parameters
{
kinViscosity 0.005148; // [m^2/s]
rho 1; // [kg/m^3]
inflowVelocity 20.0772; // [m/s]
maxLatticeVelocity 0.0337749907475931;
timesteps 150001;
diameterCylinder 1.0; // On the finest mesh
cylinderXPosition 16;
cylinderyPosition 15;
coarseMeshSize 0.05859375;
circularCrossSection true;
inflowProfile 2; // 1: Parabolic (SchaeferTurek), 2: uniform (Jacob et al)
evaluatePressure false;
pAlpha < 0.45, 0.2, 0.205 >; // points for evaluating
pOmega < 0.55, 0.2, 0.205 >; // the pressure difference
evaluateStrouhal true;
pStrouhal < 30, 23, 3 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
cylinderRefinementBuffer 0;
processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>;
// GPU related Parameters, only used if GPU is enabled
sendDirectlyFromGPU true;
gpuBlockSize <128, 1, 1>;
}
DomainSetup
{
cellsPerBlock < 64, 64, 64 >;
domainSize < 60, 30, 7.5 >;
periodic < false, false, true>;
refinementLevels 2;
numberProcesses 1; // This is for load balancing, overwritten if more than one proc is used
}
Boundaries
{
Border { direction N; walldistance -1; flag FreeSlip; }
Border { direction S; walldistance -1; flag FreeSlip; }
Border { direction W; walldistance -1; flag UBB; }
Border { direction E; walldistance -1; flag Outflow; }
}
StabilityChecker
{
checkFrequency 0;
streamOutput false;
vtkOutput true;
}
VTKWriter
{
vtkWriteFrequency 25000;
velocity true;
density false;
averageFields true;
flag false;
writeOnlySlice false;
amrFileFormat false;
oneFilePerProcess false;
}
Logging
{
logLevel info; // info progress detail tracing
writeSetupForestAndReturn true;
readSetupFromFile false;
remainingTimeLoggerFrequency 60; // in seconds
}
Evaluation
{
evaluationCheckFrequency 2000;
rampUpTime 50000;
logToStream true;
logToFile true;
filename CylinderRe3900.txt;
}
AABBRefinementSelection
{
// coordinates of Regions are always in [0,1] -> the entire simulation space is mapped to [ <0,0,0>, <1,1,1> ]
Region
{
level 2;
region [ <0.25, 0.49, 0>, <0.35, 0.51, 1> ];
}
/*Region
{
level 2;
region [ <0.26, 0.41, 0>, <0.7, 0.59, 1> ];
}
Region
{
level 3;
region [ <0.26, 0.41, 0>, <0.5, 0.59, 1> ];
}
Region
{
level 4;
region [ <0.26, 0.41, 0>, <0.4, 0.59, 1> ];
}
Region
{
level 5;
region [ <0.26, 0.41, 0>, <0.33, 0.59, 1> ];
}*/
}
\ No newline at end of file
Parameters
{
kinViscosity 0.005148; // [m^2/s]
rho 1; // [kg/m^3]
inflowVelocity 20.0772; // [m/s]
maxLatticeVelocity 0.0337749907475931;
timesteps 150001;
diameterCylinder 25; // On the finest mesh
cylinderXPosition 30;
cylinderyPosition 30;
coarseMeshSize 1;
circularCrossSection true;
inflowProfile 2; // 1: Parabolic (SchaeferTurek), 2: uniform (Jacob et al)
evaluatePressure false;
pAlpha < 0.45, 0.2, 0.205 >; // points for evaluating
pOmega < 0.55, 0.2, 0.205 >; // the pressure difference
evaluateStrouhal true;
pStrouhal < 30, 23, 3 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
cylinderRefinementBuffer 0;
processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>;
// GPU related Parameters, only used if GPU is enabled
sendDirectlyFromGPU true;
gpuBlockSize <128, 1, 1>;
}
DomainSetup
{
cellsPerBlock < 30, 30, 30 >;
domainSize < 90, 60, 30 >;
periodic < false, false, true>;
refinementLevels 2;
numberProcesses 3; // This is for load balancing, overwritten if more than one proc is used
}
Boundaries
{
Border { direction N; walldistance -1; flag FreeSlip; }
Border { direction S; walldistance -1; flag FreeSlip; }
Border { direction W; walldistance -1; flag UBB; }
Border { direction E; walldistance -1; flag Outflow; }
}
StabilityChecker
{
checkFrequency 0;
streamOutput false;
vtkOutput true;
}
VTKWriter
{
vtkWriteFrequency 25000;
velocity true;
density false;
averageFields true;
flag false;
writeOnlySlice false;
amrFileFormat false;
oneFilePerProcess false;
}
Logging
{
logLevel info; // info progress detail tracing
writeSetupForestAndReturn true;
readSetupFromFile false;
remainingTimeLoggerFrequency 60; // in seconds
}
Evaluation
{
evaluationCheckFrequency 2000;
rampUpTime 50000;
logToStream true;
logToFile true;
filename CylinderRe3900.txt;
}
AABBRefinementSelection
{
// coordinates of Regions are always in [0,1] -> the entire simulation space is mapped to [ <0,0,0>, <1,1,1> ]
/*Region
{
level 2;
region [ <0.25, 0.49, 0>, <0.35, 0.51, 1> ];
}
Region
{
level 2;
region [ <0.26, 0.41, 0>, <0.7, 0.59, 1> ];
}
Region
{
level 3;
region [ <0.26, 0.41, 0>, <0.5, 0.59, 1> ];
}
Region
{
level 4;
region [ <0.26, 0.41, 0>, <0.4, 0.59, 1> ];
}
Region
{
level 5;
region [ <0.26, 0.41, 0>, <0.33, 0.59, 1> ];
}*/
}
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file Evaluation.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "Evaluation.h"
namespace walberla
{
void Evaluation::writeMacroFields()
{
if ((executionCounter_ - uint_c(1)) % checkFrequency_ != 0) return;
auto blocks = blocks_.lock();
WALBERLA_CHECK_NOT_NULLPTR(blocks)
getMacroFields_();
const real_t Cs2 = ((setup_.dxC / setup_.dt) * (setup_.dxC / setup_.dt) / real_c(3.0)); // speed of sound squared
for (auto& block : *blocks)
{
FlagField_T* flagField = block.getData<FlagField_T>(ids_.flagField);
VelocityField_T* velocityField = block.getData<VelocityField_T>(ids_.velocityField);
ScalarField_T* densityField = block.getData<ScalarField_T>(ids_.densityField);
VelocityField_T* avgVelField = block.getData<VelocityField_T>(ids_.avgVelField);
VelocityField_T* avgVelSqrField = block.getData<VelocityField_T>(ids_.avgVelSqrField);
ScalarField_T* avgPressureField = block.getData<ScalarField_T>(ids_.avgPressureField);
auto fluid = flagField->getFlag(fluid_);
// iterate all cells of the current block
for (auto cell = flagField->beginXYZ(); cell != flagField->end(); ++cell)
{
const cell_idx_t x = cell.x();
const cell_idx_t y = cell.y();
const cell_idx_t z = cell.z();
// If the current cell is marked as being fluid ...
if (flagField->isFlagSet(x, y, z, fluid)) {
const real_t density = densityField->get(x, y, z) * setup_.rho;
const real_t pressure = (density - 1.0) * Cs2; // ρ' times speed of sound squared
avgVelField->get(x, y, z, 0) += velocityField->get(x, y, z, 0) / real_c(checkFrequency_);
avgVelSqrField->get(x, y, z, 0) += velocityField->get(x, y, z, 0) * velocityField->get(x, y, z, 0) / real_c(checkFrequency_);
avgVelField->get(x, y, z, 1) += velocityField->get(x, y, z, 1) / real_c(checkFrequency_);
avgVelSqrField->get(x, y, z, 1) += velocityField->get(x, y, z, 1) * velocityField->get(x, y, z, 1) / real_c(checkFrequency_);
avgVelField->get(x, y, z, 2) += velocityField->get(x, y, z, 2) / real_c(checkFrequency_);
avgVelSqrField->get(x, y, z, 2) += velocityField->get(x, y, z, 2) * velocityField->get(x, y, z, 2) / real_c(checkFrequency_);
avgPressureField->get(x, y, z) += pressure / real_c(checkFrequency_);
}
}
}
}
void Evaluation::operator()()
{
if (checkFrequency_ == uint_t(0)) return;
++executionCounter_;
writeMacroFields();
if (rampUpTime_ > executionCounter_) return;
real_t cDRealArea(real_c(0.0));
real_t cLRealArea(real_c(0.0));
real_t cDDiscreteArea(real_c(0.0));
real_t cLDiscreteArea(real_c(0.0));
real_t pressureDifference_L(real_c(0.0));
real_t pressureDifference(real_c(0.0));
evaluate(cDRealArea, cLRealArea, cDDiscreteArea, cLDiscreteArea, pressureDifference_L, pressureDifference);
if ((executionCounter_ - uint_c(1)) % checkFrequency_ != 0) return;
auto blocks = blocks_.lock();
WALBERLA_CHECK_NOT_NULLPTR(blocks)
// Strouhal number (needs vortex shedding frequency)
real_t vortexVelocity(real_c(0.0));
if (setup_.evaluateStrouhal)
{
auto block = blocks->getBlock(setup_.pStrouhal);
if (block != nullptr)
{
const VelocityField_T* const velocityField = block->template getData< VelocityField_T >(ids_.velocityField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pStrouhal);
WALBERLA_ASSERT(velocityField->xyzSize().contains(cell))
vortexVelocity += velocityField->get(cell.x(), cell.y(), cell.z(), cell_idx_c(0));
}
mpi::reduceInplace(vortexVelocity, mpi::SUM);
}
WALBERLA_ROOT_SECTION()
{
coefficients_[0].push_back(cDRealArea);
coefficients_[1].push_back(cLRealArea);
coefficients_[2].push_back(cDDiscreteArea);
coefficients_[3].push_back(cLDiscreteArea);
if (coefficients_[0].size() > setup_.nbrOfEvaluationPointsForCoefficientExtremas)
{
for (uint_t i = uint_t(0); i < uint_t(4); ++i)
coefficients_[i].pop_front();
}
for (uint_t i = uint_t(0); i < uint_t(4); ++i)
{
coefficientExtremas_[i] = std::make_pair(*(coefficients_[i].begin()), *(coefficients_[i].begin()));
for (auto v = coefficients_[i].begin(); v != coefficients_[i].end(); ++v)
{
coefficientExtremas_[i].first = std::min(coefficientExtremas_[i].first, *v);
coefficientExtremas_[i].second = std::max(coefficientExtremas_[i].second, *v);
}
}
std::ostringstream oss;
if (logToStream_)
{
WALBERLA_LOG_RESULT_ON_ROOT(
"force acting on cylinder (in dimensionless lattice units of the coarsest grid - evaluated in time step "
<< executionCounter_ - uint_c(1) << "):\n " << force_ << oss.str()
<< "\ndrag and lift coefficients (including extrema of last " << (coefficients_[0].size() * checkFrequency_)
<< " time steps):"
"\n \"real\" area:"
"\n c_D: "
<< cDRealArea << " (min = " << coefficientExtremas_[0].first << ", max = " << coefficientExtremas_[0].second
<< ")"
<< "\n c_L: " << cLRealArea << " (min = " << coefficientExtremas_[1].first
<< ", max = " << coefficientExtremas_[1].second << ")"
<< "\n discrete area:"
"\n c_D: "
<< cDDiscreteArea << " (min = " << coefficientExtremas_[2].first
<< ", max = " << coefficientExtremas_[2].second << ")"
<< "\n c_L: " << cLDiscreteArea << " (min = " << coefficientExtremas_[3].first
<< ", max = " << coefficientExtremas_[3].second << ")")
}
if (setup_.evaluatePressure && logToStream_)
{
WALBERLA_LOG_RESULT_ON_ROOT("pressure:\n difference: " << pressureDifference << " Pa ("
<< pressureDifference_L << ")")
}
if (setup_.evaluateStrouhal)
{
// We evaluate the derivative (-> strouhalRising_) in order to find the local minima and maxima of the velocity
// over time. If we know the time between a local minimum and a local maximum, we can calculate the frequency.
// -> "Smooth noise-robust differentiators"
// (http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/)
if (strouhalVelocities_.size() < uint_t(11)) { strouhalVelocities_.push_back(vortexVelocity); }
else
{
for (uint_t i = uint_t(0); i < uint_t(10); ++i)
strouhalVelocities_[i] = strouhalVelocities_[i + 1];
strouhalVelocities_[10] = vortexVelocity;
const real_t f1 = strouhalVelocities_[6] - strouhalVelocities_[4];
const real_t f2 = strouhalVelocities_[7] - strouhalVelocities_[3];
const real_t f3 = strouhalVelocities_[8] - strouhalVelocities_[2];
const real_t f4 = strouhalVelocities_[9] - strouhalVelocities_[1];
const real_t f5 = strouhalVelocities_[10] - strouhalVelocities_[0];
const real_t diff =
(real_c(322) * f1 + real_c(256) * f2 + real_c(39) * f3 - real_c(32) * f4 - real_c(11) * f5) /
real_c(1536);
if ((diff > real_t(0)) != strouhalRising_)
{
strouhalRising_ = (diff > real_t(0));
if (strouhalTimeStep_.size() < uint_t(3)) { strouhalTimeStep_.push_back(executionCounter_); }
else
{
strouhalTimeStep_[0] = strouhalTimeStep_[1];
strouhalTimeStep_[1] = strouhalTimeStep_[2];
strouhalTimeStep_[2] = executionCounter_;
}
}
}
if (strouhalTimeStep_.size() == uint_t(3))
{
strouhalNumberRealD_ = cylinderDiameter_ / (setup_.uMean * real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]));
strouhalNumberDiscreteD_ = real_c(D_) / (setup_.uMean * real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]));
strouhalEvaluationExecutionCount_ = executionCounter_ - uint_t(1);
if (logToStream_)
{
WALBERLA_LOG_RESULT_ON_ROOT(
"Strouhal number (evaluated in time step "
<< strouhalEvaluationExecutionCount_
<< "):"
"\n D/U (in lattice units): "
<< (cylinderDiameter_ / setup_.uMean) << " (\"real\" D), " << (real_c(D_) / setup_.uMean)
<< " (discrete D)"
"\n T: "
<< (real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]) * setup_.dt) << " s ("
<< real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0])
<< ")"
"\n f: "
<< (real_t(1) / (real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]) * setup_.dt)) << " Hz ("
<< (real_t(1) / real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]))
<< ")"
"\n St (\"real\" D): "
<< strouhalNumberRealD_ << "\n St (discrete D): " << strouhalNumberDiscreteD_)
}
}
}
if (logToFile_)
{
std::ofstream file(filename_.c_str(), std::ios_base::app);
file << executionCounter_ - uint_t(1) << " " << force_[0] << " " << force_[1] << " " << force_[2] << " "
<< cDRealArea << " " << cLRealArea << " " << cDDiscreteArea << " " << cLDiscreteArea << " "
<< pressureDifference_L << " " << pressureDifference << " " << vortexVelocity << " "
<< strouhalNumberRealD_ << " " << strouhalNumberDiscreteD_ << '\n';
file.close();
}
}
// WALBERLA_MPI_WORLD_BARRIER();
}
void Evaluation::resetForce()
{
// Supposed to be used after the timestep; after the evaluation is written
if (!initialized_) refresh();
}
void Evaluation::forceCalculation(IBlock* block, const uint_t /*level*/)
{
// Supposed to be used as a post boundary handling function on every level
if (rampUpTime_ > executionCounter_) return;
force_ += boundaryCollection_.ObstacleObject->getForce(block);
}
void Evaluation::prepareResultsForSQL()
{
real_t cDRealArea(real_t(0));
real_t cLRealArea(real_t(0));
real_t cDDiscreteArea(real_t(0));
real_t cLDiscreteArea(real_t(0));
real_t pressureDifference_L(real_t(0));
real_t pressureDifference(real_t(0));
evaluate(cDRealArea, cLRealArea, cDDiscreteArea, cLDiscreteArea, pressureDifference_L, pressureDifference);
WALBERLA_ROOT_SECTION()
{
sqlResults_["forceX_L"] = double_c(force_[0]);
sqlResults_["forceY_L"] = double_c(force_[1]);
sqlResults_["forceZ_L"] = double_c(force_[2]);
sqlResults_["forceXMin_L"] = 0.0;
sqlResults_["forceXMax_L"] = 0.0;
sqlResults_["forceXAvg_L"] = 0.0;
sqlResults_["forceXMedian_L"] = 0.0;
sqlResults_["forceXStdDev_L"] = 0.0;
sqlResults_["forceYMin_L"] = 0.0;
sqlResults_["forceYMax_L"] = 0.0;
sqlResults_["forceYAvg_L"] = 0.0;
sqlResults_["forceYMedian_L"] = 0.0;
sqlResults_["forceYStdDev_L"] = 0.0;
sqlResults_["cDRealArea"] = double_c(cDRealArea);
sqlResults_["cLRealArea"] = double_c(cLRealArea);
sqlResults_["cDDiscreteArea"] = double_c(cDDiscreteArea);
sqlResults_["cLDiscreteArea"] = double_c(cLDiscreteArea);
sqlResults_["cDRealAreaMin"] = double_c(coefficientExtremas_[0].first);
sqlResults_["cDRealAreaMax"] = double_c(coefficientExtremas_[0].second);
sqlResults_["cLRealAreaMin"] = double_c(coefficientExtremas_[1].first);
sqlResults_["cLRealAreaMax"] = double_c(coefficientExtremas_[1].second);
sqlResults_["cDDiscreteAreaMin"] = double_c(coefficientExtremas_[2].first);
sqlResults_["cDDiscreteAreaMax"] = double_c(coefficientExtremas_[2].second);
sqlResults_["cLDiscreteAreaMin"] = double_c(coefficientExtremas_[3].first);
sqlResults_["cLDiscreteAreaMax"] = double_c(coefficientExtremas_[3].second);
sqlResults_["pressureDifference_L"] = double_c(pressureDifference_L);
sqlResults_["pressureDifference"] = double_c(pressureDifference);
sqlResults_["strouhalNumberRealD"] = double_c(strouhalNumberRealD_);
sqlResults_["strouhalNumberDiscreteD"] = double_c(strouhalNumberDiscreteD_);
}
}
void Evaluation::getResultsForSQLOnRoot(std::map< std::string, double >& realProperties,
std::map< std::string, int >& integerProperties) const
{
WALBERLA_ROOT_SECTION()
{
for (auto result = sqlResults_.begin(); result != sqlResults_.end(); ++result)
realProperties[result->first] = result->second;
integerProperties["forceEvaluationTimeStep"] = int_c(executionCounter_ - uint_c(1));
integerProperties["strouhalEvaluationTimeStep"] = int_c(strouhalEvaluationExecutionCount_);
}
}
void Evaluation::refresh()
{
auto blocks = blocks_.lock();
WALBERLA_CHECK_NOT_NULLPTR(blocks)
// Calculate obstacle surface areas required for evaluating drag and lift force
real_t yMin(setup_.H);
real_t yMax(real_t(0));
real_t AD(real_t(0));
real_t AL(real_t(0));
for (auto block = blocks->begin(); block != blocks->end(); ++block)
{
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
auto fluid = flagField->getFlag(fluid_);
auto obstacle = flagField->getFlag(obstacle_);
const real_t area = real_c(1.0);
auto xyzSize = flagField->xyzSize();
for (auto z = xyzSize.zMin(); z <= xyzSize.zMax(); ++z)
{
for (auto y = xyzSize.yMin(); y <= xyzSize.yMax(); ++y)
{
for (auto x = xyzSize.xMin(); x <= xyzSize.xMax(); ++x)
{
if (flagField->isFlagSet(x, y, z, fluid))
{
for (auto it = Stencil_T::beginNoCenter(); it != Stencil_T::end(); ++it)
{
auto nx = x + cell_idx_c(it.cx());
auto ny = y + cell_idx_c(it.cy());
auto nz = z + cell_idx_c(it.cz());
if (flagField->isFlagSet(nx, ny, nz, obstacle))
{
const Vector3< real_t > p = blocks->getBlockLocalCellCenter(*block, Cell(nx, ny, nz));
if (it.cx() == 1 && it.cy() == 0 && it.cz() == 0) { AD += area; }
else if (it.cx() == 0 && it.cz() == 0)
{
if (it.cy() == -1) { yMax = std::max(yMax, p[1]); }
else if (it.cy() == 1)
{
yMin = std::min(yMin, p[1]);
AL += area;
}
}
}
}
}
}
}
}
}
mpi::reduceInplace(yMin, mpi::MIN);
mpi::reduceInplace(yMax, mpi::MAX);
mpi::reduceInplace(AD, mpi::SUM);
mpi::reduceInplace(AL, mpi::SUM);
WALBERLA_ROOT_SECTION()
{
const Cell yMinCell = blocks->getCell(real_t(0), yMin, real_t(0));
const Cell yMaxCell = blocks->getCell(real_t(0), yMax, real_t(0));
D_ = uint_c(yMaxCell[1] - yMinCell[1]) + uint_t(1);
AD_ = AD;
AL_ = AL;
}
// Check if points alpha and omega (required for evaluating the pressure difference) are located in fluid cells
if (setup_.evaluatePressure)
{
int alpha(0);
int omega(0);
auto block = blocks->getBlock(setup_.pAlpha);
if (block != nullptr)
{
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pAlpha);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(fluid_);
if (!flagField->isFlagSet(cell, fluid))
{
const auto aabb = blocks->getBlockLocalCellAABB(*block, cell);
const Vector3< real_t > pAlpha = setup_.pAlpha;
setup_.pAlpha[0] = aabb.xMin() - aabb.xSize() / real_t(2);
WALBERLA_LOG_WARNING("Cell for evaluating pressure difference at point alpha "
<< pAlpha
<< " is not a fluid cell!"
"\nChanging point alpha to "
<< setup_.pAlpha << " ...")
}
else { alpha = 1; }
}
block = blocks->getBlock(setup_.pOmega);
if (block != nullptr)
{
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pOmega);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(fluid_);
if (!flagField->isFlagSet(cell, fluid))
{
const auto aabb = blocks->getBlockLocalCellAABB(*block, cell);
const Vector3< real_t > pOmega = setup_.pOmega;
setup_.pOmega[0] = aabb.xMax() + aabb.xSize() / real_t(2);
WALBERLA_LOG_WARNING("Cell for evaluating pressure difference at point omega "
<< pOmega
<< " is not a fluid cell!"
"\nChanging point omega to "
<< setup_.pOmega << " ...")
}
else { omega = 1; }
}
mpi::allReduceInplace(alpha, mpi::SUM);
mpi::allReduceInplace(omega, mpi::SUM);
if (alpha == 0)
{
block = blocks->getBlock(setup_.pAlpha);
if (block != nullptr)
{
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pAlpha);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(fluid_);
if (!flagField->isFlagSet(cell, fluid))
WALBERLA_ABORT("Cell for evaluating pressure difference at point alpha "
<< setup_.pAlpha << " is still not a fluid cell!")
alpha = 1;
}
mpi::reduceInplace(alpha, mpi::SUM);
}
if (omega == 0)
{
block = blocks->getBlock(setup_.pOmega);
if (block != nullptr)
{
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pOmega);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(fluid_);
if (!flagField->isFlagSet(cell, fluid))
WALBERLA_ABORT("Cell for evaluating pressure difference at point omega "
<< setup_.pOmega << " is still not a fluid cell!")
omega = 1;
}
mpi::reduceInplace(omega, mpi::SUM);
}
WALBERLA_ROOT_SECTION()
{
if (alpha == 0)
WALBERLA_ABORT(
"Point alpha "
<< setup_.pAlpha
<< " (required for evaluating the pressure difference) is not located inside the fluid domain!")
WALBERLA_ASSERT_EQUAL(alpha, 1)
if (omega == 0)
WALBERLA_ABORT(
"Point omega "
<< setup_.pOmega
<< " (required for evaluating the pressure difference) is not located inside the fluid domain!")
WALBERLA_ASSERT_EQUAL(omega, 1)
}
}
// Check if point for evaluating the Strouhal number is located inside a fluid cell
if (setup_.evaluateStrouhal)
{
int strouhal(0);
auto block = blocks->getBlock(setup_.pStrouhal);
if (block != nullptr)
{
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pStrouhal);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(fluid_);
if (!flagField->isFlagSet(cell, fluid))
WALBERLA_ABORT("Cell for evaluating the Strouhal number at point " << setup_.pStrouhal
<< " is not a fluid cell!")
strouhal = 1;
}
mpi::reduceInplace(strouhal, mpi::SUM);
WALBERLA_ROOT_SECTION()
{
if (strouhal == 0)
WALBERLA_ABORT(
"Point " << setup_.pStrouhal
<< " (required for evaluating the Strouhal number) is not located inside the fluid domain!")
WALBERLA_ASSERT_EQUAL(strouhal, 1)
}
}
initialized_ = true;
}
void Evaluation::evaluate(real_t& cDRealArea, real_t& cLRealArea, real_t& cDDiscreteArea, real_t& cLDiscreteArea,
real_t& pressureDifference_L, real_t& pressureDifference)
{
if (!initialized_) refresh();
// force on obstacle
mpi::reduceInplace(force_, mpi::SUM);
// pressure difference
real_t pAlpha(real_c(0.0));
real_t pOmega(real_c(0.0));
auto blocks = blocks_.lock();
WALBERLA_CHECK_NOT_NULLPTR(blocks)
if (setup_.evaluatePressure)
{
auto block = blocks->getBlock(setup_.pAlpha);
if (block != nullptr)
{
const ScalarField_T* const densityField = block->template getData< ScalarField_T >(ids_.densityField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pAlpha);
WALBERLA_ASSERT(densityField->xyzSize().contains(cell))
pAlpha += densityField->get(cell) / real_c(3.0);
}
block = blocks->getBlock(setup_.pOmega);
if (block != nullptr)
{
const ScalarField_T* const densityField = block->template getData< ScalarField_T >(ids_.densityField);
const auto cell = blocks->getBlockLocalCell(*block, setup_.pOmega);
WALBERLA_ASSERT(densityField->xyzSize().contains(cell))
pOmega += densityField->get(cell) / real_c(3.0);
}
mpi::reduceInplace(pAlpha, mpi::SUM);
mpi::reduceInplace(pOmega, mpi::SUM);
}
WALBERLA_ROOT_SECTION()
{
cDRealArea = (real_c(2.0) * force_[0]) / (setup_.uMean * setup_.uMean * cylinderDiameter_ * cylinderHeight_);
cLRealArea = (real_c(2.0) * force_[1]) / (setup_.uMean * setup_.uMean * cylinderDiameter_ * cylinderHeight_);
cDDiscreteArea = (real_c(2.0) * force_[0]) / (setup_.uMean * setup_.uMean * AD_);
cLDiscreteArea = (real_c(2.0) * force_[1]) / (setup_.uMean * setup_.uMean * AL_);
DragCoefficient DC(executionCounter_ - uint_c(1), force_, cDRealArea, cLRealArea, cDDiscreteArea, cLDiscreteArea);
dragResults.push_back(DC);
pressureDifference_L = pAlpha - pOmega;
pressureDifference = (pressureDifference_L * setup_.rho * setup_.dxC * setup_.dxC) / (setup_.dt * setup_.dt);
}
force_[0] = real_c(0.0);
force_[1] = real_c(0.0);
force_[2] = real_c(0.0);
}
}
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file Evaluation.h
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
# pragma once
#include "core/config/Config.h"
#include "core/math/Sample.h"
#include "lbm_generated/field/PdfField.h"
#include "sqlite/SQLite.h"
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <functional>
#include <iostream>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
#include "Setup.h"
#include "Types.h"
#include "FlowAroundCylinderInfoHeader.h"
using namespace walberla;
using StorageSpecification_T = lbm::FlowAroundCylinderStorageSpecification;
using Stencil_T = StorageSpecification_T::Stencil;
using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
using FlagField_T = FlagField< uint8_t >;
using BoundaryCollection_T = lbm::FlowAroundCylinderBoundaryCollection< FlagField_T >;
using VoidFunction = std::function<void ()>;
namespace walberla
{
struct DragCoefficient {
uint_t timestep;
double Fx;
double Fy;
double Fz;
double cDRealArea;
double cLRealArea;
double cDDiscreteArea;
double cLDiscreteArea;
DragCoefficient(uint_t t, Vector3<double> f, double cdR, double clR, double cdD, double clD) : timestep(t), Fx(f[0]), Fy(f[1]), Fz(f[2]), cDRealArea(cdR), cLRealArea(clR), cDDiscreteArea(cdD), cLDiscreteArea(clD) {}
};
class Evaluation
{
public:
Evaluation(const weak_ptr< StructuredBlockStorage >& blocks, const uint_t checkFrequency, const uint_t rampUpTime,
VoidFunction& getMacroFields, BoundaryCollection_T & boundaryCollection,
const IDs& ids, const FlagUID& fluid, const FlagUID& obstacle, const Setup& setup,
const bool logToStream = true, const bool logToFile = true,
const std::string& filename = std::string("FlowAroundCylinder.txt"))
: blocks_(blocks), executionCounter_(uint_t(0)), checkFrequency_(checkFrequency), rampUpTime_(rampUpTime),
getMacroFields_(getMacroFields), boundaryCollection_(boundaryCollection), ids_(ids), fluid_(fluid), obstacle_(obstacle), setup_(setup),
strouhalNumberRealD_(real_t(0)), strouhalEvaluationExecutionCount_(uint_t(0)),
logToStream_(logToStream), logToFile_(logToFile), filename_(filename)
{
coefficients_.resize(uint_t(4));
coefficientExtremas_.resize(uint_t(4));
WALBERLA_ROOT_SECTION()
{
if (logToFile_)
{
std::ofstream file(filename_.c_str());
file << "# time step [1], force (x) [2], force (y) [3], force (z) [4], "
"cD (real area) [5], cL (real area) [6], cD (discrete area) [7], cL (discrete area) [8], "
"pressure difference (in lattice units) [9], pressure difference (in Pa) [10], vortex velocity (in "
"lattice units) [11], "
"Strouhal number (real D) [12], Strouhal number (discrete D) [13]"
<< '\n';
if (!setup_.evaluatePressure)
file << "# ATTENTION: pressure was not evaluated, pressure difference is set to zero!" << '\n';
if (!setup_.evaluateStrouhal)
file << "# ATTENTION: vortex velocities were not evaluated, Strouhal number is set to zero!"
<< '\n';
file.close();
}
}
refresh();
cylinderDiameter_ = real_c(2.0) * setup_.cylinderRadius / setup_.dxF;
cylinderHeight_ = setup_.H / setup_.dxF;
WALBERLA_LOG_INFO_ON_ROOT(
"Evaluation initialised:"
"\n + Cylinder real area: " << cylinderDiameter_ * cylinderHeight_
<< "\n + Cylinder real diameter: " << cylinderDiameter_
<< "\n + Cylinder real height: " << cylinderHeight_
<< "\n + Cylinder discrete area: " << AD_
)
}
void operator()();
void forceCalculation(IBlock* block, const uint_t level); // for calculating the force
void resetForce();
void writeMacroFields();
void writeDrag()
{
filesystem::path dataFile( "dragCoefficient.csv" );
if( filesystem::exists( dataFile ) )
std::remove( dataFile.string().c_str() );
std::ofstream outfile( "dragCoefficient.csv" );
outfile << "SEP=," << "\n";
outfile << "timestep," << "Fx," << "Fy," << "Fz," << "cDRealArea," << "cLRealArea," << "cDDiscreteArea," << "cLDiscreteArea" << "\n";
for(auto it = dragResults.begin(); it != dragResults.end(); ++it)
{
outfile << it->timestep << ",";
outfile << it->Fx << "," << it->Fy << "," << it->Fz << ",";
outfile << it->cDRealArea << ",";
outfile << it->cLRealArea << ",";
outfile << it->cDDiscreteArea << ",";
outfile << it->cLDiscreteArea;
outfile << "\n";
}
outfile.close();
};
std::function<void (IBlock *, const uint_t)> forceCalculationFunctor()
{
return [this](IBlock* block, uint_t level) { forceCalculation(block, level); };
}
std::function<void()> resetForceFunctor()
{
return [this]() { resetForce(); };
}
void prepareResultsForSQL();
void getResultsForSQLOnRoot(std::map< std::string, double >& realProperties,
std::map< std::string, int >& integerProperties) const;
void refresh();
protected:
void evaluate(real_t& cDRealArea, real_t& cLRealArea, real_t& cDDiscreteArea, real_t& cLDiscreteArea,
real_t& pressureDifference_L, real_t& pressureDifference);
bool initialized_{false};
weak_ptr< StructuredBlockStorage > blocks_;
uint_t executionCounter_;
uint_t checkFrequency_;
uint_t rampUpTime_;
VoidFunction & getMacroFields_;
BoundaryCollection_T & boundaryCollection_;
IDs ids_;
FlagUID fluid_;
FlagUID obstacle_;
Setup setup_;
uint_t D_;
real_t AD_;
real_t AL_;
std::vector<DragCoefficient> dragResults;
real_t cylinderDiameter_;
real_t cylinderHeight_;
Vector3< real_t > force_;
std::vector< std::deque< real_t > > coefficients_;
std::vector< std::pair< real_t, real_t > > coefficientExtremas_;
std::vector< real_t > strouhalVelocities_;
std::vector< uint_t > strouhalTimeStep_;
bool strouhalRising_;
real_t strouhalNumberRealD_;
real_t strouhalNumberDiscreteD_;
uint_t strouhalEvaluationExecutionCount_;
bool logToStream_;
bool logToFile_;
std::string filename_;
std::map< std::string, double > sqlResults_;
}; // class Evaluation
class EvaluationRefresh
{
public:
EvaluationRefresh(const weak_ptr< Evaluation >& evaluation) : evaluation_(evaluation) {}
void operator()()
{
auto evaluation = evaluation_.lock();
WALBERLA_CHECK_NOT_NULLPTR(evaluation)
evaluation->refresh();
}
void operator()(BlockForest&, const PhantomBlockForest&)
{
auto evaluation = evaluation_.lock();
WALBERLA_CHECK_NOT_NULLPTR(evaluation)
evaluation->refresh();
}
private:
weak_ptr< Evaluation > evaluation_;
};
}
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file FlowAroundCylinder.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "blockforest/AABBRefinementSelection.h"
#include "blockforest/SetupBlockForest.h"
#include "blockforest/StructuredBlockForest.h"
#include "blockforest/communication/NonUniformBufferedScheme.h"
#include "blockforest/loadbalancing/StaticCurve.h"
#include "core/Abort.h"
#include "core/DataTypes.h"
#include "core/MemoryUsage.h"
#include "core/SharedFunctor.h"
#include "core/debug/CheckFunctions.h"
#include "core/logging/Initialization.h"
#include "core/logging/Logging.h"
#include "core/math/Vector3.h"
#include "core/mpi/Environment.h"
#include "core/mpi/MPIManager.h"
#include "core/mpi/Reduce.h"
#include "core/timing/RemainingTimeLogger.h"
#include "core/timing/TimingPool.h"
#include "field/AddToStorage.h"
#include "field/CellCounter.h"
#include "field/FlagField.h"
#include "field/StabilityChecker.h"
#include "field/adaptors/AdaptorCreators.h"
#include "field/vtk/FlagFieldCellFilter.h"
#include "field/vtk/VTKWriter.h"
#include "geometry/InitBoundaryHandling.h"
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
# include "gpu/AddGPUFieldToStorage.h"
# include "gpu/DeviceSelectMPI.h"
# include "gpu/ErrorChecking.h"
# include "gpu/FieldCopy.h"
# include "gpu/HostFieldAllocator.h"
# include "gpu/communication/NonUniformGPUScheme.h"
# include "gpu/communication/UniformGPUScheme.h"
#endif
#include "lbm_generated/communication/NonuniformGeneratedPdfPackInfo.h"
#include "lbm_generated/communication/UniformGeneratedPdfPackInfo.h"
#include "lbm_generated/evaluation/PerformanceEvaluation.h"
#include "lbm_generated/field/AddToStorage.h"
#include "lbm_generated/field/PdfField.h"
#include "lbm_generated/refinement/BasicRecursiveTimeStep.h"
#include "lbm_generated/refinement/RefinementScaling.h"
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
# include "lbm_generated/gpu/AddToStorage.h"
# include "lbm_generated/gpu/BasicRecursiveTimeStepGPU.h"
# include "lbm_generated/gpu/GPUPdfField.h"
# include "lbm_generated/gpu/NonuniformGeneratedGPUPdfPackInfo.h"
# include "lbm_generated/gpu/UniformGeneratedGPUPdfPackInfo.h"
#endif
#include "stencil/D3Q27.h"
#include "timeloop/SweepTimeloop.h"
#include "vtk/VTKOutput.h"
#include <cmath>
#include <cstdlib>
#include <functional>
#include <iostream>
#include <memory>
#include <vector>
#include "Cylinder.h"
#include "Evaluation.h"
#include "FlowAroundCylinderInfoHeader.h"
#include "FlowAroundCylinderStaticDefines.h"
#include "Setup.h"
#include "Types.h"
using namespace walberla;
using StorageSpecification_T = lbm::FlowAroundCylinderStorageSpecification;
using Stencil_T = StorageSpecification_T::Stencil;
using CommunicationStencil_T = StorageSpecification_T::CommunicationStencil;
using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
using FlagField_T = FlagField< uint8_t >;
using BoundaryCollection_T = lbm::FlowAroundCylinderBoundaryCollection< FlagField_T >;
using SweepCollection_T = lbm::FlowAroundCylinderSweepCollection;
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
using GPUPdfField_T = lbm_generated::GPUPdfField< StorageSpecification_T >;
using gpu::communication::NonUniformGPUScheme;
using gpu::communication::UniformGPUScheme;
using lbm_generated::NonuniformGeneratedGPUPdfPackInfo;
using lbm_generated::UniformGeneratedGPUPdfPackInfo;
#else
using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
using blockforest::communication::NonUniformBufferedScheme;
using blockforest::communication::UniformBufferedScheme;
using lbm_generated::NonuniformGeneratedPdfPackInfo;
using lbm_generated::UniformGeneratedPdfPackInfo;
#endif
using RefinementSelectionFunctor = SetupBlockForest::RefinementSelectionFunction;
using BlockExclusionFunction = SetupBlockForest::BlockExclusionFunction;
//////////////////////
// Parameter Struct //
//////////////////////
namespace
{
void workloadMemoryAndSUIDAssignment(SetupBlockForest& forest, const memory_t memoryPerBlock, const Setup& setup)
{
for (auto block = forest.begin(); block != forest.end(); ++block)
{
block->setWorkload(numeric_cast< workload_t >(uint_t(1) << block->getLevel()));
block->setMemory(memoryPerBlock);
}
}
shared_ptr< SetupBlockForest >
createSetupBlockForest(const blockforest::RefinementSelectionFunctions& refinementSelectionFunctions,
const BlockExclusionFunction& blockExclusionFunction, const Setup& setup,
uint_t numberOfProcesses, const memory_t memoryPerCell, const memory_t processMemoryLimit,
const bool outputSetupForest)
{
shared_ptr< SetupBlockForest > forest = make_shared< SetupBlockForest >();
const memory_t memoryPerBlock = numeric_cast< memory_t >((setup.cellsPerBlock[0] + uint_t(2) * setup.numGhostLayers) *
(setup.cellsPerBlock[1] + uint_t(2) * setup.numGhostLayers) *
(setup.cellsPerBlock[2] + uint_t(2) * setup.numGhostLayers)) * memoryPerCell;
forest->addRefinementSelectionFunction(refinementSelectionFunctions);
forest->addBlockExclusionFunction(blockExclusionFunction);
forest->addWorkloadMemorySUIDAssignmentFunction(
std::bind(workloadMemoryAndSUIDAssignment, std::placeholders::_1, memoryPerBlock, std::cref(setup)));
forest->init(AABB(real_c(0), real_c(0), real_c(0), real_c(setup.domainSize[0]), real_c(setup.domainSize[1]), real_c(setup.domainSize[2])),
setup.xBlocks, setup.yBlocks, setup.zBlocks, setup.periodic[0], setup.periodic[1], setup.periodic[2]);
MPIManager::instance()->useWorldComm();
forest->balanceLoad(blockforest::StaticLevelwiseCurveBalanceWeighted(), numberOfProcesses);
// forest->balanceLoad(blockforest::StaticLevelwiseCurveBalanceWeighted(true),
// numberOfProcesses, real_t(0), processMemoryLimit, true, false);
if (outputSetupForest)
{
forest->writeVTKOutput("domain_decomposition");
forest->writeCSV("process_distribution");
}
WALBERLA_LOG_INFO_ON_ROOT("SetupBlockForest created successfully:\n" << *forest);
return forest;
}
shared_ptr< blockforest::StructuredBlockForest >
createStructuredBlockForest(const blockforest::RefinementSelectionFunctions& refinementSelectionFunctions,
const BlockExclusionFunction& blockExclusionFunction, const Setup& setup,
const memory_t memoryPerCell, const memory_t processMemoryLimit)
{
// if (configBlock.isDefined("sbffile"))
// {
// std::string sbffile = configBlock.getParameter< std::string >("sbffile");
//
// WALBERLA_LOG_INFO_ON_ROOT("Creating the block structure: loading from file \'" << sbffile << "\' ...");
//
// MPIManager::instance()->useWorldComm();
//
// auto bf = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->rank()), sbffile.c_str(), true,
// false);
//
// auto sbf = std::make_shared< StructuredBlockForest >(bf, setup.xCells, setup.yCells, setup.zCells);
// sbf->createCellBoundingBoxes();
//
// return sbf;
// }
WALBERLA_LOG_INFO_ON_ROOT("Creating the block structure ...");
shared_ptr< SetupBlockForest > sforest =
createSetupBlockForest(refinementSelectionFunctions, blockExclusionFunction, setup,
uint_c(MPIManager::instance()->numProcesses()), memoryPerCell, processMemoryLimit, false);
auto bf = std::make_shared< blockforest::BlockForest >(uint_c(MPIManager::instance()->rank()), *sforest, false);
auto sbf = std::make_shared< blockforest::StructuredBlockForest >(bf, setup.cellsPerBlock[0], setup.cellsPerBlock[1],
setup.cellsPerBlock[2]);
sbf->createCellBoundingBoxes();
return sbf;
}
void consistentlySetBoundary(const std::shared_ptr< StructuredBlockForest >& blocks, Block& block,
FlagField_T* flagField, const uint8_t flag,
const std::function< bool(const Vector3< real_t >&) >& isBoundary)
{
const uint_t level = blocks->getLevel(block);
int ghostLayers = int_c(flagField->nrOfGhostLayers());
CellInterval cells = flagField->xyzSize();
cells.expand(cell_idx_c(ghostLayers));
std::vector< CellInterval > coarseRegions;
for (auto dir = stencil::D3Q27::beginNoCenter(); dir != stencil::D3Q27::end(); ++dir)
{
const auto index = blockforest::getBlockNeighborhoodSectionIndex(dir.cx(), dir.cy(), dir.cz());
if (block.neighborhoodSectionHasLargerBlock(index))
{
CellInterval coarseRegion(cells);
for (uint_t i = 0; i != 3; ++i)
{
const auto c = stencil::c[i][*dir];
if (c == -1)
coarseRegion.max()[i] = coarseRegion.min()[i] + cell_idx_c(2 * ghostLayers - 1);
else if (c == 1)
coarseRegion.min()[i] = coarseRegion.max()[i] - cell_idx_c(2 * ghostLayers - 1);
}
coarseRegions.push_back(coarseRegion);
}
}
for (auto cell = cells.begin(); cell != cells.end(); ++cell)
{
bool inCoarseRegion(false);
for (auto region = coarseRegions.begin(); region != coarseRegions.end() && !inCoarseRegion; ++region)
inCoarseRegion = region->contains(*cell);
if (!inCoarseRegion)
{
Vector3< real_t > center;
blocks->getBlockLocalCellCenter(block, *cell, center);
blocks->mapToPeriodicDomain(center);
if (isBoundary(center)) { flagField->addFlag(cell->x(), cell->y(), cell->z(), flag); }
}
else
{
Cell globalCell(*cell);
blocks->transformBlockLocalToGlobalCell(globalCell, block);
Cell coarseCell(globalCell);
for (uint_t i = 0; i < 3; ++i)
{
if (coarseCell[i] < cell_idx_t(0)) { coarseCell[i] = -((cell_idx_t(1) - coarseCell[i]) >> 1); }
else { coarseCell[i] >>= 1; }
}
Vector3< real_t > coarseCenter;
blocks->getCellCenter(coarseCenter, coarseCell, level - uint_t(1));
blocks->mapToPeriodicDomain(coarseCenter);
if (isBoundary(coarseCenter)) { flagField->addFlag(cell->x(), cell->y(), cell->z(), flag); }
}
}
}
void setupBoundaryCylinder(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID,
const FlagUID& obstacleFlagUID, const std::function< bool(const Vector3< real_t >&) >& isObstacleBoundary)
{
for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
{
Block& b = dynamic_cast< Block& >(*bIt);
auto flagField = b.getData< FlagField_T >(flagFieldID);
uint8_t obstacleFlag = flagField->registerFlag(obstacleFlagUID);
consistentlySetBoundary(sbfs, b, flagField, obstacleFlag, isObstacleBoundary);
}
}
class InflowProfile
{
public:
InflowProfile( const real_t velocity, const real_t H, const uint_t inflowProfile ) : velocity_(velocity), H_( H ), inflowProfile_( inflowProfile )
{
uTerm_ = ( real_c(16.0) * velocity );
HTerm_ = ( real_c(1.0) / ( H * H * H * H ) );
tConstTerm_ = uTerm_ * HTerm_;
}
Vector3< real_t > operator()( const Cell& pos, const shared_ptr< StructuredBlockForest >& SbF, IBlock& block ) const;
private:
real_t velocity_;
real_t H_;
real_t uTerm_;
real_t HTerm_;
real_t tConstTerm_;
uint_t inflowProfile_;
}; // class InflowProfile
Vector3< real_t > InflowProfile::operator()( const Cell& pos, const shared_ptr< StructuredBlockForest >& SbF, IBlock& block ) const
{
if (inflowProfile_ == 1)
{
Cell globalCell;
real_t x;
real_t y;
real_t z;
const uint_t level = SbF->getLevel(block);
SbF->transformBlockLocalToGlobalCell(globalCell, block, pos);
SbF->getCellCenter(x, y, z, globalCell, level);
return Vector3< real_t >(tConstTerm_ * y * z * ( H_ - y ) * ( H_ - z ), real_c(0.0), real_c(0.0) );
}
else if (inflowProfile_ == 2)
{
return Vector3< real_t >(velocity_, real_c(0.0), real_c(0.0) );
}
else
{
WALBERLA_ABORT("Inflow profile not implemented")
}
}
} // namespace
//////////////////////
// Parameter Struct //
//////////////////////
int main(int argc, char** argv)
{
mpi::Environment env(argc, argv);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::selectDeviceBasedOnMpiRank();
#endif
shared_ptr< Config > config = make_shared< Config >();
config->readParameterFile(argv[1]);
logging::configureLogging(config);
///////////////////////
/// PARAMETER INPUT ///
///////////////////////
// read general simulation parameters
auto parameters = config->getOneBlock("Parameters");
const real_t kinViscosity = parameters.getParameter< real_t >("kinViscosity");
const real_t rho = parameters.getParameter< real_t >("rho");
const real_t referenceVelocity = parameters.getParameter< real_t >("inflowVelocity");
const real_t maxLatticeVelocity = parameters.getParameter< real_t >("maxLatticeVelocity");
const uint_t inflowProfile = parameters.getParameter< uint_t >("inflowProfile");
const real_t diameterCylinder = parameters.getParameter< real_t >("diameterCylinder");
const real_t coarseMeshSize = parameters.getParameter< real_t >("coarseMeshSize");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps");
const real_t dt = maxLatticeVelocity / referenceVelocity * coarseMeshSize;
const real_t latticeViscosity = kinViscosity / coarseMeshSize / coarseMeshSize * dt;
const real_t omega = real_c(real_c(1.0) / (real_c(3.0) * latticeViscosity + real_c(0.5)));
const real_t uMean = inflowProfile == 1 ? real_c(4.0) / real_c(9.0) * maxLatticeVelocity : maxLatticeVelocity;
const real_t reynoldsNumber = (uMean * (diameterCylinder / coarseMeshSize)) / latticeViscosity;
// read domain parameters
auto domainParameters = config->getOneBlock("DomainSetup");
const uint_t refinementLevels = domainParameters.getParameter< uint_t >("refinementLevels");
const Vector3< uint_t > cellsPerBlock = domainParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
const Vector3< real_t > domainSize = domainParameters.getParameter< Vector3< real_t > >("domainSize");
const Vector3< bool > periodic = domainParameters.getParameter< Vector3< bool > >("periodic");
const uint_t numGhostLayers = uint_c(2);
const real_t fineMeshSize = real_c(coarseMeshSize) / real_c(1 << refinementLevels);
WALBERLA_LOG_INFO_ON_ROOT("Diameter of the Cylinder is resolved with " << diameterCylinder / fineMeshSize << " lattice cells.")
auto loggingParameters = config->getOneBlock("Logging");
bool writeSetupForestAndReturn = loggingParameters.getParameter< bool >("writeSetupForestAndReturn", false);
if (uint_c(MPIManager::instance()->numProcesses()) > 1) writeSetupForestAndReturn = false;
for (uint_t i = 0; i < 3; ++i)
{
auto cells = uint_c(std::ceil((domainSize[i] / coarseMeshSize) / real_c(cellsPerBlock[i]))) * cellsPerBlock[i];
if (cells > uint_c(std::ceil(domainSize[i] / coarseMeshSize)))
{
WALBERLA_LOG_WARNING_ON_ROOT("Total cells in direction "
<< i << ": " << uint_c(domainSize[i] / coarseMeshSize)
<< " is not dividable by cellsPerBlock in that direction. Domain was extended")
}
}
const uint_t xBlocks = uint_c(std::ceil((domainSize[0] / coarseMeshSize) / real_c(cellsPerBlock[0])));
const uint_t yBlocks = uint_c(std::ceil((domainSize[1] / coarseMeshSize) / real_c(cellsPerBlock[1])));
const uint_t zBlocks = uint_c(std::ceil((domainSize[2] / coarseMeshSize) / real_c(cellsPerBlock[2])));
Setup setup;
setup.xBlocks = xBlocks;
setup.yBlocks = yBlocks;
setup.zBlocks = zBlocks;
setup.xCells = xBlocks * cellsPerBlock[0];
setup.yCells = yBlocks * cellsPerBlock[1];
setup.zCells = zBlocks * cellsPerBlock[2];
setup.cellsPerBlock = cellsPerBlock;
setup.domainSize = domainSize;
setup.periodic = periodic;
setup.numGhostLayers = numGhostLayers;
setup.H = domainSize[2];
setup.L = domainSize[0];
setup.cylinderXPosition = parameters.getParameter< real_t >("cylinderXPosition");
setup.cylinderYPosition = parameters.getParameter< real_t >("cylinderYPosition");
setup.cylinderRadius = diameterCylinder / real_c(2.0);
setup.circularCrossSection = parameters.getParameter< bool >("circularCrossSection");
setup.evaluateForceComponents = false;
setup.nbrOfEvaluationPointsForCoefficientExtremas = 100;
setup.evaluatePressure = parameters.getParameter< bool >("evaluatePressure");
setup.pAlpha = parameters.getParameter< Vector3< real_t > >("pAlpha", Vector3< real_t >(real_c(0.45), real_c(0.2), real_c(0.205)));
setup.pOmega = parameters.getParameter< Vector3< real_t > >("pOmega", Vector3< real_t >(real_c(0.55), real_c(0.2), real_c(0.205)));
setup.evaluateStrouhal = parameters.getParameter< bool >("evaluateStrouhal");
setup.pStrouhal = parameters.getParameter< Vector3< real_t > >("pStrouhal", Vector3< real_t >(real_c(1), real_c(0.325), real_c(0.205)));
setup.kinViscosity = kinViscosity;
setup.rho = rho;
setup.inflowVelocity = referenceVelocity;
setup.uMean = uMean;
setup.dxC = coarseMeshSize;
setup.dxF = fineMeshSize;
setup.dt = dt;
const uint_t valuesPerCell = (Stencil_T::Q + VelocityField_T::F_SIZE + uint_c(2) * ScalarField_T::F_SIZE);
const uint_t sizePerValue = sizeof(PdfField_T::value_type);
const memory_t memoryPerCell = memory_t(valuesPerCell * sizePerValue + uint_c(1));
const memory_t processMemoryLimit = parameters.getParameter< memory_t >("processMemoryLimit", memory_t(512)) * memory_t(1024 * 1024);
///////////////////////////
/// Refinement ///
///////////////////////////
blockforest::RefinementSelectionFunctions refinementSelectionFunctions;
blockforest::AABBRefinementSelection aabbRefinementSelection(config->getOneBlock("AABBRefinementSelection"));
refinementSelectionFunctions.add(aabbRefinementSelection);
const real_t cylinderRefinementBuffer = parameters.getParameter< real_t >("cylinderRefinementBuffer", real_t(0));
Cylinder cylinder(setup);
CylinderBlockExclusion cylinderBlockExclusion(cylinder);
CylinderRefinementSelection cylinderRefinementSelection(cylinder, refinementLevels, cylinderRefinementBuffer);
refinementSelectionFunctions.add(cylinderRefinementSelection);
///////////////////////////
/// CREATE BLOCK FOREST ///
///////////////////////////
if (writeSetupForestAndReturn)
{
std::string sbffile = "sbfFlowAroundCylinder.bfs";
std::ostringstream infoString;
infoString << "You have selected the option of just creating the block structure (= domain decomposition) and "
"saving the result to file\n"
"by specifying the output file name \'"
<< sbffile << "\' AND also specifying \'saveToFile\'.\n";
if (MPIManager::instance()->numProcesses() > 1)
WALBERLA_ABORT(infoString.str() << "In this mode you need to start " << argv[0] << " with just one process!")
WALBERLA_LOG_INFO_ON_ROOT(infoString.str() << "Creating the block structure ...")
const uint_t numberProcesses = domainParameters.getParameter< uint_t >("numberProcesses");
shared_ptr< SetupBlockForest > sforest =
createSetupBlockForest(refinementSelectionFunctions, cylinderBlockExclusion, setup, numberProcesses,
memoryPerCell, processMemoryLimit, true);
sforest->saveToFile(sbffile.c_str());
WALBERLA_LOG_INFO_ON_ROOT("Benchmark run data:"
"\n- simulation parameters:"
"\n + collision model: "
<< infoCollisionOperator << "\n + stencil: " << infoStencil
<< "\n + streaming: " << infoStreamingPattern
<< "\n + compressible: " << (StorageSpecification_T::compressible ? "yes" : "no")
<< "\n + mesh levels: " << refinementLevels + uint_c(1)
<< "\n + resolution: " << coarseMeshSize << " - on the coarsest grid"
<< "\n + resolution: " << fineMeshSize << " - on the finest grid"
<< "\n- simulation properties:"
"\n + H: "
<< setup.H << " [m]"
<< "\n + L: " << setup.L << " [m]"
<< "\n + cylinder pos.(x): " << setup.cylinderXPosition << " [m]"
<< "\n + cylinder pos.(y): " << setup.cylinderYPosition << " [m]"
<< "\n + cylinder radius: " << setup.cylinderRadius << " [m]"
<< "\n + circular profile: " << (setup.circularCrossSection ? "yes" : "no (= box)")
<< "\n + kin. viscosity: " << setup.kinViscosity << " [m^2/s] (" << setup.kinViscosity
<< " - on the coarsest grid)"
<< "\n + omega: " << omega << " on the coarsest grid"
<< "\n + rho: " << setup.rho << " [kg/m^3]"
<< "\n + inflow velocity: " << setup.inflowVelocity << " [m/s]"
<< "\n + lattice velocity: " << maxLatticeVelocity
<< "\n + Reynolds number: " << reynoldsNumber
<< "\n + dt (coarsest grid): " << setup.dt << " [s]"
<< "\n + #time steps: " << (real_t(1) / setup.dt) << " (for 1s of real time)")
logging::Logging::printFooterOnStream();
return EXIT_SUCCESS;
}
auto blocks = createStructuredBlockForest(refinementSelectionFunctions, cylinderBlockExclusion, setup, memoryPerCell,
processMemoryLimit);
////////////////////////////////////
/// CREATE AND INITIALIZE FIELDS ///
////////////////////////////////////
// create fields
const StorageSpecification_T StorageSpec = StorageSpecification_T();
IDs ids;
ids.avgVelField = field::addToStorage<VelocityField_T>(blocks, "average velocity", real_t(0.0), field::fzyx, numGhostLayers);
ids.avgVelSqrField = field::addToStorage<VelocityField_T>(blocks, "average velocity squared", real_t(0.0), field::fzyx, numGhostLayers);
ids.avgPressureField = field::addToStorage<ScalarField_T>(blocks, "average pressure", real_t(0.0), field::fzyx, numGhostLayers);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
auto allocator = make_shared< gpu::HostFieldAllocator< real_t > >();
ids.pdfField = lbm_generated::addPdfFieldToStorage(blocks, "pdfs", StorageSpec, numGhostLayers, field::fzyx);
ids.velocityField = field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx, numGhostLayers);
ids.densityField = field::addToStorage< ScalarField_T >(blocks, "density", setup.rho, field::fzyx, numGhostLayers);
ids.flagField = field::addFlagFieldToStorage< FlagField_T >(blocks, "Boundary Flag Field", uint_c(3));
ids.pdfFieldGPU = lbm_generated::addGPUPdfFieldToStorage< PdfField_T >(blocks, ids.pdfField, StorageSpec, "pdfs on GPU", true);
ids.velocityFieldGPU = gpu::addGPUFieldToStorage< VelocityField_T >(blocks, ids.velocityField, "velocity on GPU", true);
ids.densityFieldGPU = gpu::addGPUFieldToStorage< ScalarField_T >(blocks, ids.densityField, "density on GPU", true);
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#else
ids.pdfField = lbm_generated::addPdfFieldToStorage(blocks, "pdfs", StorageSpec, numGhostLayers, field::fzyx);
ids.velocityField = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx, numGhostLayers);
ids.densityField = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx, numGhostLayers);
ids.flagField = field::addFlagFieldToStorage< FlagField_T >(blocks, "Boundary Flag Field", uint_c(3));
#endif
WALBERLA_MPI_BARRIER()
const Cell innerOuterSplit =
Cell(parameters.getParameter< Vector3< cell_idx_t > >("innerOuterSplit", Vector3< cell_idx_t >(1, 1, 1)));
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
const Vector3< int64_t > gpuBlockSize = parameters.getParameter< Vector3< int64_t > >("gpuBlockSize");
SweepCollection_T sweepCollection(blocks, ids.pdfFieldGPU, ids.densityFieldGPU, ids.velocityFieldGPU, gpuBlockSize[0],
gpuBlockSize[1], gpuBlockSize[2], omega, innerOuterSplit);
for (auto& block : *blocks)
{
sweepCollection.initialise(&block, cell_idx_c(1));
}
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
#else
SweepCollection_T sweepCollection(blocks, ids.pdfField, ids.densityField, ids.velocityField, omega, innerOuterSplit);
for (auto& block : *blocks)
{
sweepCollection.initialise(&block, cell_idx_c(1));
}
#endif
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< PdfField_T, gpu::GPUField< real_t > >(blocks, ids.pdfField, ids.pdfFieldGPU);
#endif
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
WALBERLA_LOG_INFO_ON_ROOT("Setting up communication")
std::shared_ptr< NonUniformGPUScheme< CommunicationStencil_T > > nonUniformCommunication =
std::make_shared< NonUniformGPUScheme< CommunicationStencil_T > >(blocks);
std::shared_ptr< NonuniformGeneratedGPUPdfPackInfo< GPUPdfField_T > > nonUniformPackInfo =
lbm_generated::setupNonuniformGPUPdfCommunication< GPUPdfField_T >(blocks, ids.pdfFieldGPU);
nonUniformCommunication->addPackInfo(nonUniformPackInfo);
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#else
WALBERLA_LOG_INFO_ON_ROOT("Setting up communication...")
auto nonUniformCommunication = std::make_shared< NonUniformBufferedScheme< CommunicationStencil_T > >(blocks);
auto nonUniformPackInfo = lbm_generated::setupNonuniformPdfCommunication< PdfField_T >(blocks, ids.pdfField);
nonUniformCommunication->addPackInfo(nonUniformPackInfo);
#endif
WALBERLA_MPI_BARRIER()
WALBERLA_LOG_INFO_ON_ROOT("Setting up communication done")
/////////////////////////
/// BOUNDARY HANDLING ///
/////////////////////////
WALBERLA_LOG_INFO_ON_ROOT("Start BOUNDARY HANDLING")
// create and initialize boundary handling
const FlagUID fluidFlagUID("Fluid");
const FlagUID obstacleFlagUID("Obstacle");
auto boundariesConfig = config->getBlock("Boundaries");
geometry::initBoundaryHandling< FlagField_T >(*blocks, ids.flagField, boundariesConfig);
setupBoundaryCylinder(blocks, ids.flagField, obstacleFlagUID, cylinder);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, fluidFlagUID, cell_idx_c(0));
std::function< real_t(const Cell&, const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
wallDistanceFunctor = wallDistance(cylinder);
InflowProfile const velocityCallback{maxLatticeVelocity, setup.H, inflowProfile};
std::function< Vector3< real_t >(const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
velocity_initialisation = velocityCallback;
const real_t omegaFinestLevel = lbm_generated::relaxationRateScaling(omega, refinementLevels);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, fluidFlagUID, omegaFinestLevel,
wallDistanceFunctor, velocity_initialisation, ids.pdfField);
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#else
BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfField, fluidFlagUID, omegaFinestLevel,
wallDistanceFunctor, velocity_initialisation);
#endif
WALBERLA_MPI_BARRIER()
WALBERLA_LOG_INFO_ON_ROOT("BOUNDARY HANDLING done")
//////////////////////////////////
/// SET UP SWEEPS AND TIMELOOP ///
//////////////////////////////////
WALBERLA_LOG_INFO_ON_ROOT("Start SWEEPS AND TIMELOOP")
// flow evaluation
auto EvaluationParameters = config->getOneBlock("Evaluation");
const uint_t evaluationCheckFrequency = EvaluationParameters.getParameter< uint_t >("evaluationCheckFrequency");
const uint_t rampUpTime = EvaluationParameters.getParameter< uint_t >("rampUpTime");
const bool evaluationLogToStream = EvaluationParameters.getParameter< bool >("logToStream");
const bool evaluationLogToFile = EvaluationParameters.getParameter< bool >("logToFile");
const std::string evaluationFilename = EvaluationParameters.getParameter< std::string >("filename");
std::function< void() > getMacroFields = [&]() {
for (auto& block : *blocks)
{
sweepCollection.calculateMacroscopicParameters(&block);
}
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< VelocityField_T, gpu::GPUField< real_t > >(blocks, ids.velocityField, ids.velocityFieldGPU);
gpu::fieldCpy< ScalarField_T, gpu::GPUField< real_t > >(blocks, ids.densityField, ids.densityFieldGPU);
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif
};
shared_ptr< Evaluation > evaluation(new Evaluation(blocks, evaluationCheckFrequency, rampUpTime, getMacroFields, boundaryCollection, ids,
fluidFlagUID, obstacleFlagUID, setup, evaluationLogToStream, evaluationLogToFile, evaluationFilename));
// create time loop
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
std::shared_ptr< lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollection_T, BoundaryCollection_T > >
LBMRefinement;
LBMRefinement = std::make_shared<
lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollection_T, BoundaryCollection_T > >(
blocks, ids.pdfFieldGPU, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
LBMRefinement->addPostBoundaryHandlingBlockFunction(evaluation->forceCalculationFunctor());
LBMRefinement->addRefinementToTimeLoop(timeloop);
#else
SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
std::shared_ptr< lbm_generated::BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T > >
LBMRefinement;
LBMRefinement =
std::make_shared< lbm_generated::BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T > >(
blocks, ids.pdfField, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
LBMRefinement->addPostBoundaryHandlingBlockFunction(evaluation->forceCalculationFunctor());
LBMRefinement->addRefinementToTimeLoop(timeloop);
#endif
//////////////////
/// VTK OUTPUT ///
//////////////////
WALBERLA_LOG_INFO_ON_ROOT("SWEEPS AND TIMELOOP done")
auto VTKWriter = config->getOneBlock("VTKWriter");
const uint_t vtkWriteFrequency = VTKWriter.getParameter< uint_t >("vtkWriteFrequency", 0);
const bool writeVelocity = VTKWriter.getParameter< bool >("velocity");
const bool writeDensity = VTKWriter.getParameter< bool >("density");
const bool writeAverageFields = VTKWriter.getParameter< bool >("averageFields", false);
const bool writeFlag = VTKWriter.getParameter< bool >("flag");
const bool writeOnlySlice = VTKWriter.getParameter< bool >("writeOnlySlice", true);
const bool amrFileFormat = VTKWriter.getParameter< bool >("amrFileFormat", false);
const bool oneFilePerProcess = VTKWriter.getParameter< bool >("oneFilePerProcess", false);
auto finalDomain = blocks->getDomain();
if (vtkWriteFrequency > 0)
{
auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false,
"vtk_FlowAroundCylinder", "simulation_step", false, true, true,
false, 0, amrFileFormat, oneFilePerProcess);
vtkOutput->addBeforeFunction([&]() {
for (auto& block : *blocks)
{
sweepCollection.calculateMacroscopicParameters(&block);
}
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< VelocityField_T, gpu::GPUField< real_t > >(blocks, ids.velocityField, ids.velocityFieldGPU);
gpu::fieldCpy< ScalarField_T, gpu::GPUField< real_t > >(blocks, ids.densityField, ids.densityFieldGPU);
#endif
});
if (writeOnlySlice)
{
const AABB sliceAABB(finalDomain.xMin(), finalDomain.yMin(), finalDomain.center()[2] - coarseMeshSize,
finalDomain.xMax(), finalDomain.yMax(), finalDomain.center()[2] + coarseMeshSize);
vtkOutput->addCellInclusionFilter(vtk::AABBCellFilter(sliceAABB));
}
if(!amrFileFormat)
{
field::FlagFieldCellFilter<FlagField_T> fluidFilter( ids.flagField );
fluidFilter.addFlag(obstacleFlagUID);
vtkOutput->addCellExclusionFilter(fluidFilter);
}
if (writeVelocity)
{
auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(ids.velocityField, "velocity");
vtkOutput->addCellDataWriter(velWriter);
}
if (writeDensity)
{
auto densityWriter = make_shared< field::VTKWriter< ScalarField_T, float32 > >(ids.densityField, "density");
vtkOutput->addCellDataWriter(densityWriter);
}
if (writeAverageFields)
{
auto avgVelWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(ids.avgVelField, "avgVelocity");
vtkOutput->addCellDataWriter(avgVelWriter);
auto avgVelSqrWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(ids.avgVelSqrField, "avgVelocitySqr");
vtkOutput->addCellDataWriter(avgVelSqrWriter);
auto avgPressureWriter = make_shared< field::VTKWriter< ScalarField_T, float32 > >(ids.avgPressureField, "avgPressure");
vtkOutput->addCellDataWriter(avgPressureWriter);
}
if (writeFlag)
{
auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(ids.flagField, "flag");
vtkOutput->addCellDataWriter(flagWriter);
}
timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
}
// log remaining time
const real_t remainingTimeLoggerFrequency =
loggingParameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0); // in seconds
if (uint_c(remainingTimeLoggerFrequency) > 0)
{
timeloop.addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
"remaining time logger");
}
// LBM stability check
auto CheckerParameters = config->getOneBlock("StabilityChecker");
const uint_t checkFrequency = CheckerParameters.getParameter< uint_t >("checkFrequency", uint_t(0));
if (checkFrequency > 0)
{
auto checkFunction = [](PdfField_T::value_type value) { return value < math::abs(PdfField_T::value_type(10)); };
timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
config, blocks, ids.pdfField, ids.flagField, fluidFlagUID, checkFunction)),
"Stability check");
}
timeloop.addFuncBeforeTimeStep(SharedFunctor< Evaluation >(evaluation), "evaluation");
timeloop.addFuncBeforeTimeStep(evaluation->resetForceFunctor(), "evaluation: reset force");
// WALBERLA_LOG_INFO_ON_ROOT("Execute single timestep to fully complete the preprocessing")
// Do a single timestep to make sure all setups are completed before benchmarking
// timeloop.singleStep();
// timeloop.setCurrentTimeStepToZero();
WALBERLA_MPI_BARRIER()
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif
// WALBERLA_LOG_INFO_ON_ROOT("Execute single timestep to fully complete the preprocessing done")
//////////////////////
/// RUN SIMULATION ///
//////////////////////
const lbm_generated::PerformanceEvaluation< FlagField_T > performance(blocks, ids.flagField, fluidFlagUID);
field::CellCounter< FlagField_T > fluidCells(blocks, ids.flagField, fluidFlagUID);
fluidCells();
WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << blocks->getNumberOfBlocks())
for (uint_t level = 0; level <= refinementLevels; level++)
{
WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << blocks->getNumberOfBlocks(level))
}
WALBERLA_LOG_INFO_ON_ROOT(
"Benchmark run data:"
"\n- simulation parameters:"
"\n + collision model: "
<< infoCollisionOperator << "\n + stencil: " << infoStencil << "\n + streaming: "
<< infoStreamingPattern << "\n + compressible: " << (StorageSpecification_T::compressible ? "yes" : "no")
<< "\n + mesh levels: " << refinementLevels + uint_c(1) << "\n + resolution: " << coarseMeshSize
<< " - on the coarsest grid"
<< "\n + resolution: " << fineMeshSize << " - on the finest grid"
<< "\n- simulation properties:"
"\n + fluid cells: "
<< fluidCells.numberOfCells() << " (in total on all levels)"
<< "\n + H: " << setup.H << " [m]"
<< "\n + L: " << setup.L << " [m]"
<< "\n + cylinder pos.(x): " << setup.cylinderXPosition << " [m]"
<< "\n + cylinder pos.(y): " << setup.cylinderYPosition << " [m]"
<< "\n + cylinder radius: " << setup.cylinderRadius << " [m]"
<< "\n + circular profile: " << (setup.circularCrossSection ? "yes" : "no (= box)")
<< "\n + kin. viscosity: " << setup.kinViscosity << " [m^2/s] (on the coarsest grid)"
<< "\n + omega: " << omega << " (on the coarsest grid)"
<< "\n + rho: " << setup.rho << " [kg/m^3]"
<< "\n + inflow velocity: " << setup.inflowVelocity << " [m/s]"
<< "\n + lattice velocity: " << maxLatticeVelocity
<< "\n + Reynolds number: " << reynoldsNumber
<< "\n + dx (coarsest grid): " << setup.dxC << " [m]"
<< "\n + dx (finest grid): " << setup.dxF << " [m]"
<< "\n + dt (coarsest grid): " << setup.dt << " [s]"
<< "\n + #time steps: " << timeloop.getNrOfTimeSteps() << " on the coarsest grid, "
<< (real_t(1) / setup.dt) << " for 1s of real time)"
<< "\n + simulation time: " << (real_c(timeloop.getNrOfTimeSteps()) * setup.dt) << " [s]")
WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation")
WcTimingPool timeloopTiming;
WcTimer simTimer;
WALBERLA_MPI_BARRIER()
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif
simTimer.start();
timeloop.run(timeloopTiming);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif
WALBERLA_MPI_BARRIER()
simTimer.end();
evaluation->writeDrag();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
real_t time = simTimer.max();
WALBERLA_MPI_SECTION() { walberla::mpi::reduceInplace(time, walberla::mpi::MAX); }
performance.logResultOnRoot(timesteps, time);
const auto reducedTimeloopTiming = timeloopTiming.getReduced();
WALBERLA_LOG_RESULT_ON_ROOT("Time loop timing:\n" << *reducedTimeloopTiming)
printResidentMemoryStatistics();
return EXIT_SUCCESS;
}
\ No newline at end of file
import sympy as sp
import numpy as np
from pystencils import TypedSymbol, Target
from pystencils.field import fields
from pystencils.simp.subexpression_insertion import insert_constants, insert_aliases
from lbmpy import Stencil, LBStencil, Method, LBMConfig, LBMOptimisation
from lbmpy.boundaries.boundaryconditions import ExtrapolationOutflow, UBB, QuadraticBounceBack, FreeSlip, NoSlip
from lbmpy.creationfunctions import create_lb_collision_rule
from pystencils_walberla import CodeGeneration, generate_info_header
from lbmpy_walberla import generate_lbm_package, lbm_boundary_generator
info_header = """
#pragma once
const char * infoStencil = "{stencil}";
const char * infoStreamingPattern = "{streaming_pattern}";
const char * infoCollisionOperator = "{collision_operator}";
"""
omega = sp.symbols("omega")
with CodeGeneration() as ctx:
dtype = 'float64'
pdf_dtype = 'float64'
stencil = LBStencil(Stencil.D3Q27)
q = stencil.Q
dim = stencil.D
streaming_pattern = 'aa'
pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {pdf_dtype}[3D]", layout='fzyx')
velocity_field, density_field = fields(f"velocity({dim}), density(1) : {dtype}[{dim}D]", layout='fzyx')
macroscopic_fields = {'density': density_field, 'velocity': velocity_field}
method_enum = Method.CUMULANT
lbm_config = LBMConfig(
method=method_enum,
stencil=stencil,
relaxation_rate=omega,
compressible=True,
galilean_correction=False,
fourth_order_correction=1e-4 if method_enum == Method.CUMULANT and stencil.Q == 27 else False,
field_name='pdfs',
streaming_pattern=streaming_pattern,
)
lbm_opt = LBMOptimisation(cse_global=False, cse_pdfs=False)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
collision_rule = insert_constants(collision_rule)
collision_rule = insert_aliases(collision_rule)
lb_method = collision_rule.method
if ctx.gpu:
target = Target.GPU
openmp = False
cpu_vec = None
vp = [('int64_t', 'cudaBlockSize0'),
('int64_t', 'cudaBlockSize1'),
('int64_t', 'cudaBlockSize2')]
sweep_block_size = (TypedSymbol("cudaBlockSize0", np.int64),
TypedSymbol("cudaBlockSize1", np.int64),
TypedSymbol("cudaBlockSize2", np.int64))
sweep_params = {'block_size': sweep_block_size}
else:
if ctx.optimize_for_localhost:
cpu_vec = {"instruction_set": None}
else:
cpu_vec = None
openmp = True if ctx.openmp else False
target = Target.CPU
sweep_params = {}
vp = ()
freeslip = lbm_boundary_generator("FreeSlip", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil))
no_slip = lbm_boundary_generator(class_name='NoSlip', flag_uid='NoSlip',
boundary_object=NoSlip(), field_data_type=pdf_dtype)
quadratic_bounce_back = QuadraticBounceBack(omega, calculate_force_on_boundary=True)
no_slip_interpolated = lbm_boundary_generator(class_name='Obstacle', flag_uid='Obstacle',
boundary_object=quadratic_bounce_back, field_data_type=pdf_dtype)
ubb = lbm_boundary_generator(class_name='UBB', flag_uid='UBB',
boundary_object=UBB(lambda *args: None, dim=stencil.D, data_type=dtype),
field_data_type=pdf_dtype)
outflow = lbm_boundary_generator(class_name='Outflow', flag_uid='Outflow',
boundary_object=ExtrapolationOutflow(stencil[4], lb_method),
field_data_type=pdf_dtype)
generate_lbm_package(ctx, name="FlowAroundCylinder", collision_rule=collision_rule,
lbm_config=lbm_config, lbm_optimisation=lbm_opt,
nonuniform=True, boundaries=[freeslip, no_slip, no_slip_interpolated, ubb, outflow],
macroscopic_fields=macroscopic_fields, set_pre_collision_pdfs=True,
gpu_indexing_params=sweep_params,
target=target, data_type=dtype, pdfs_data_type=pdf_dtype,
cpu_vectorize_info=cpu_vec)
field_typedefs = {'VelocityField_T': velocity_field,
'ScalarField_T': density_field}
# Info header containing correct template definitions for stencil and field
generate_info_header(ctx, 'FlowAroundCylinderInfoHeader',
field_typedefs=field_typedefs)
infoHeaderParams = {
'stencil': stencil.name.lower(),
'streaming_pattern': streaming_pattern,
'collision_operator': lbm_config.method.name.lower(),
}
ctx.write_file("FlowAroundCylinderStaticDefines.h", info_header.format(**infoHeaderParams))
Parameters
{
kinViscosity 0.001; // [m^2/s]
rho 1; // [kg/m^3]
inflowVelocity 0.45; // [m/s]
maxLatticeVelocity 0.01;
timesteps 20001;
diameterCylinder 0.1; // On the finest mesh
cylinderXPosition 0.5;
cylinderyPosition 0.2;
coarseMeshSize 0.01; // 9.765625e-4;
circularCrossSection true;
inflowProfile 1; // 1: Parabolic (SchaeferTurek), 2: uniform (Jacob et al)
evaluatePressure true;
pAlpha < 0.45, 0.2, 0.205 >; // points for evaluating
pOmega < 0.55, 0.2, 0.205 >; // the pressure difference
evaluateStrouhal false;
pStrouhal < 1, 0.325, 0.205 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
cylinderRefinementBuffer 0;
processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>;
// GPU related Parameters, only used if GPU is enabled
sendDirectlyFromGPU true;
gpuBlockSize <128, 1, 1>;
}
DomainSetup
{
cellsPerBlock < 250, 41, 41 >;
domainSize < 2.5, 0.41, 0.41 >;
periodic < false, false, false>;
refinementLevels 0;
numberProcesses 1; // This is for load balancing, overwritten if more than one proc is used
}
Boundaries
{
Border { direction N; walldistance -1; flag NoSlip; }
Border { direction S; walldistance -1; flag NoSlip; }
Border { direction T; walldistance -1; flag NoSlip; }
Border { direction B; walldistance -1; flag NoSlip; }
Border { direction W; walldistance -1; flag UBB; }
Border { direction E; walldistance -1; flag Outflow; }
}
StabilityChecker
{
checkFrequency 0;
streamOutput false;
vtkOutput true;
}
VTKWriter
{
vtkWriteFrequency 0;
velocity true;
density false;
averageFields true;
flag false;
writeOnlySlice true;
amrFileFormat false;
oneFilePerProcess false;
}
Logging
{
logLevel info; // info progress detail tracing
writeSetupForestAndReturn false;
readSetupFromFile false;
remainingTimeLoggerFrequency 60; // in seconds
}
Evaluation
{
evaluationCheckFrequency 2000;
rampUpTime 10000;
logToStream true;
logToFile true;
filename SchaeferTurek3D-1Z.txt;
}
AABBRefinementSelection
{
// coordinates of Regions are always in [0,1] -> the entire simulation space is mapped to [ <0,0,0>, <1,1,1> ]
Region
/*{
level 1;
region [ <0.2, 0.2, 0>, <0.8, 0.8, 1> ];
}
Region
{
level 2;
region [ <0.26, 0.41, 0>, <0.7, 0.59, 1> ];
}
Region
{
level 3;
region [ <0.26, 0.41, 0>, <0.5, 0.59, 1> ];
}
Region
{
level 4;
region [ <0.26, 0.41, 0>, <0.4, 0.59, 1> ];
}
Region
{
level 5;
region [ <0.26, 0.41, 0>, <0.33, 0.59, 1> ];
}*/
}
\ No newline at end of file
Parameters
{
kinViscosity 0.001; // [m^2/s]
rho 1; // [kg/m^3]
inflowVelocity 2.25; // [m/s]
maxLatticeVelocity 0.01;
timesteps 500001;
diameterCylinder 0.1; // On the finest mesh
cylinderXPosition 0.5;
cylinderyPosition 0.2;
coarseMeshSize 9.765625e-4;
circularCrossSection true;
inflowProfile 1; // 1: Parabolic (SchaeferTurek), 2: uniform (Jacob et al)
evaluatePressure false;
pAlpha < 0.45, 0.2, 0.205 >; // points for evaluating
pOmega < 0.55, 0.2, 0.205 >; // the pressure difference
evaluateStrouhal true;
pStrouhal < 1, 0.325, 0.205 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
cylinderRefinementBuffer 0;
processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>;
// GPU related Parameters, only used if GPU is enabled
sendDirectlyFromGPU true;
gpuBlockSize <128, 1, 1>;
}
DomainSetup
{
cellsPerBlock < 320, 420, 420 >;
domainSize < 2.5, 0.41, 0.41 >;
periodic < false, false, false>;
refinementLevels 0;
numberProcesses 1; // This is for load balancing, overwritten if more than one proc is used
}
Boundaries
{
Border { direction N; walldistance -1; flag NoSlip; }
Border { direction S; walldistance -1; flag NoSlip; }
Border { direction T; walldistance -1; flag NoSlip; }
Border { direction B; walldistance -1; flag NoSlip; }
Border { direction W; walldistance -1; flag UBB; }
Border { direction E; walldistance -1; flag Outflow; }
}
StabilityChecker
{
checkFrequency 0;
streamOutput false;
vtkOutput true;
}
VTKWriter
{
vtkWriteFrequency 50000;
velocity true;
density false;
flag false;
writeOnlySlice true;
amrFileFormat false;
oneFilePerProcess false;
}
Logging
{
logLevel info; // info progress detail tracing
writeSetupForestAndReturn true;
readSetupFromFile false;
remainingTimeLoggerFrequency 60; // in seconds
}
Evaluation
{
evaluationCheckFrequency 5000;
rampUpTime 300000;
logToStream true;
logToFile true;
filename SchaeferTurek3D-2Z.txt;
}
AABBRefinementSelection
{
// coordinates of Regions are always in [0,1] -> the entire simulation space is mapped to [ <0,0,0>, <1,1,1> ]
Region
/*{
level 1;
region [ <0.2, 0.2, 0>, <0.8, 0.8, 1> ];
}
Region
{
level 2;
region [ <0.26, 0.41, 0>, <0.7, 0.59, 1> ];
}
Region
{
level 3;
region [ <0.26, 0.41, 0>, <0.5, 0.59, 1> ];
}
Region
{
level 4;
region [ <0.26, 0.41, 0>, <0.4, 0.59, 1> ];
}
Region
{
level 5;
region [ <0.26, 0.41, 0>, <0.33, 0.59, 1> ];
}*/
}
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file Setup.h
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#pragma once
#include "core/DataTypes.h"
#include "core/math/Vector3.h"
namespace walberla
{
struct Setup
{
uint_t xBlocks;
uint_t yBlocks;
uint_t zBlocks;
uint_t xCells;
uint_t yCells;
uint_t zCells;
Vector3< uint_t > cellsPerBlock;
Vector3< real_t > domainSize;
Vector3< bool > periodic;
uint_t numGhostLayers;
real_t H;
real_t L;
real_t cylinderXPosition;
real_t cylinderYPosition;
real_t cylinderRadius;
bool circularCrossSection;
bool evaluateForceComponents;
uint_t nbrOfEvaluationPointsForCoefficientExtremas;
bool evaluatePressure;
Vector3< real_t > pAlpha;
Vector3< real_t > pOmega;
bool evaluateStrouhal;
Vector3< real_t > pStrouhal;
real_t kinViscosity;
real_t rho;
real_t inflowVelocity;
real_t uMean;
real_t dxC;
real_t dxF;
real_t dt;
};
}
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file Types.h
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
# pragma once
#include "domain_decomposition/BlockDataID.h"
using namespace walberla;
struct IDs
{
BlockDataID pdfField;
BlockDataID velocityField;
BlockDataID densityField;
BlockDataID avgVelField;
BlockDataID avgVelSqrField;
BlockDataID avgPressureField;
BlockDataID flagField;
BlockDataID pdfFieldGPU;
BlockDataID velocityFieldGPU;
BlockDataID densityFieldGPU;
};
waLBerla_link_files_to_builddir( *.py )
waLBerla_link_files_to_builddir( *.prm )
waLBerla_link_files_to_builddir( *.png )
waLBerla_link_files_to_builddir( *.obj )
waLBerla_link_files_to_builddir( *.stl )
waLBerla_link_files_to_builddir( *.mtl )
waLBerla_generate_target_from_python(NAME FlowAroundSphereCodeGen
FILE FlowAroundSphere.py
OUT_FILES FlowAroundSphereSweepCollection.h FlowAroundSphereSweepCollection.${CODEGEN_FILE_SUFFIX}
FlowAroundSphereStorageSpecification.h FlowAroundSphereStorageSpecification.${CODEGEN_FILE_SUFFIX}
FreeSlip.h FreeSlip.${CODEGEN_FILE_SUFFIX}
NoSlip.h NoSlip.${CODEGEN_FILE_SUFFIX}
Obstacle.h Obstacle.${CODEGEN_FILE_SUFFIX}
UBB.h UBB.${CODEGEN_FILE_SUFFIX}
Outflow.h Outflow.${CODEGEN_FILE_SUFFIX}
FixedDensity.h FixedDensity.${CODEGEN_FILE_SUFFIX}
FlowAroundSphereBoundaryCollection.h
FlowAroundSphereInfoHeader.h
FlowAroundSphereStaticDefines.h)
if (WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
waLBerla_add_executable ( NAME FlowAroundSphere
FILES FlowAroundSphere.cpp Sphere.cpp Evaluation.cpp GridGeneration.h Types.h
DEPENDS blockforest boundary core field gpu lbm_generated stencil timeloop vtk FlowAroundSphereCodeGen )
else()
waLBerla_add_executable ( NAME FlowAroundSphere
FILES FlowAroundSphere.cpp Sphere.cpp Evaluation.cpp GridGeneration.h Types.h
DEPENDS blockforest boundary core field lbm_generated stencil timeloop vtk FlowAroundSphereCodeGen )
endif(WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
//======================================================================================================================
//
// 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/>.
//
//! \file Evaluation.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "Evaluation.h"
namespace walberla
{
void Evaluation::operator()()
{
if (checkFrequency_ == uint_t(0)) return;
++coarseExecutionCounter_;
if (rampUpTime_ > coarseExecutionCounter_) return;
real_t pressureDifference_L(real_t(0));
real_t pressureDifference(real_t(0));
evaluate(pressureDifference_L, pressureDifference);
if ((coarseExecutionCounter_ - uint_c(1)) % checkFrequency_ != 0) return;
WALBERLA_CHECK_NOT_NULLPTR(blocks_)
// Strouhal number (needs vortex shedding frequency)
real_t vortexVelocity(real_t(0));
if (setup_.evaluateStrouhal){
auto block = blocks_->getBlock(setup_.pStrouhal);
if (block != nullptr){
const VelocityField_T* const velocityField = block->template getData< VelocityField_T >(ids_.velocityField);
const auto cell = blocks_->getBlockLocalCell(*block, setup_.pStrouhal);
WALBERLA_ASSERT(velocityField->xyzSize().contains(cell))
vortexVelocity += velocityField->get(cell.x(), cell.y(), cell.z(), cell_idx_c(0));
}
mpi::reduceInplace(vortexVelocity, mpi::SUM);
}
WALBERLA_ROOT_SECTION()
{
for(auto it = dragResults.begin(); it != dragResults.end(); ++it){
coefficients_[0].push_back(it->cDRealArea);
coefficients_[1].push_back(it->cLRealArea);
coefficients_[2].push_back(it->cDDiscreteArea);
coefficients_[3].push_back(it->cLDiscreteArea);
}
if (coefficients_[0].size() > setup_.nbrOfEvaluationPointsForCoefficientExtremas){
for (uint_t i = uint_t(0); i < uint_t(4); ++i)
coefficients_[i].pop_front();
}
for (uint_t i = uint_t(0); i < uint_t(4); ++i){
coefficientExtremas_[i] = std::make_pair(*(coefficients_[i].begin()), *(coefficients_[i].begin()));
for (auto v = coefficients_[i].begin(); v != coefficients_[i].end(); ++v){
coefficientExtremas_[i].first = std::min(coefficientExtremas_[i].first, *v);
coefficientExtremas_[i].second = std::max(coefficientExtremas_[i].second, *v);
}
}
std::ostringstream oss;
if (logToStream_ and !dragResults.empty()){
WALBERLA_LOG_RESULT_ON_ROOT(
"force acting on sphere (in dimensionless lattice units of the coarsest grid - evaluated in time step "
<< coarseExecutionCounter_ - uint_c(1) << "):\n " << force_ << oss.str()
<< "\ndrag and lift coefficients (including extrema of last " << (coefficients_[0].size() * checkFrequency_)
<< " time steps):"
"\n \"real\" area:"
"\n c_D: "
<< dragResults.back().cDRealArea << " (min = " << coefficientExtremas_[0].first << ", max = " << coefficientExtremas_[0].second
<< ")"
<< "\n c_L: " << dragResults.back().cLRealArea << " (min = " << coefficientExtremas_[1].first
<< ", max = " << coefficientExtremas_[1].second << ")"
<< "\n discrete area:"
"\n c_D: "
<< dragResults.back().cDDiscreteArea << " (min = " << coefficientExtremas_[2].first
<< ", max = " << coefficientExtremas_[2].second << ")"
<< "\n c_L: " << dragResults.back().cLDiscreteArea << " (min = " << coefficientExtremas_[3].first
<< ", max = " << coefficientExtremas_[3].second << ")")
}
if (setup_.evaluatePressure && logToStream_){
WALBERLA_LOG_RESULT_ON_ROOT("pressure:\n difference: " << pressureDifference << " Pa ("<< pressureDifference_L << ")")
}
if (setup_.evaluateStrouhal)
{
// We evaluate the derivative (-> strouhalRising_) in order to find the local minima and maxima of the velocity
// over time. If we know the time between a local minimum and a local maximum, we can calculate the frequency.
// -> "Smooth noise-robust differentiators"
// (http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/)
if (strouhalVelocities_.size() < uint_t(11)) { strouhalVelocities_.push_back(vortexVelocity); }
else{
for (uint_t i = uint_t(0); i < uint_t(10); ++i)
strouhalVelocities_[i] = strouhalVelocities_[i + 1];
strouhalVelocities_[10] = vortexVelocity;
const real_t f1 = strouhalVelocities_[6] - strouhalVelocities_[4];
const real_t f2 = strouhalVelocities_[7] - strouhalVelocities_[3];
const real_t f3 = strouhalVelocities_[8] - strouhalVelocities_[2];
const real_t f4 = strouhalVelocities_[9] - strouhalVelocities_[1];
const real_t f5 = strouhalVelocities_[10] - strouhalVelocities_[0];
const real_t diff =
(real_c(322) * f1 + real_c(256) * f2 + real_c(39) * f3 - real_c(32) * f4 - real_c(11) * f5) /
real_c(1536);
if ((diff > real_t(0)) != strouhalRising_){
strouhalRising_ = (diff > real_t(0));
if (strouhalTimeStep_.size() < uint_t(3)) { strouhalTimeStep_.push_back(coarseExecutionCounter_); }
else{
strouhalTimeStep_[0] = strouhalTimeStep_[1];
strouhalTimeStep_[1] = strouhalTimeStep_[2];
strouhalTimeStep_[2] = coarseExecutionCounter_;
}
}
}
if (strouhalTimeStep_.size() == uint_t(3)){
strouhalNumberRealD_ = diameterSphere / (meanVelocity * real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]));
strouhalEvaluationExecutionCount_ = coarseExecutionCounter_ - uint_t(1);
if (logToStream_){
WALBERLA_LOG_RESULT_ON_ROOT(
"Strouhal number (evaluated in time step "
<< strouhalEvaluationExecutionCount_
<< "):"
"\n D/U (in lattice units): "
<< (diameterSphere / meanVelocity) << " (\"real\" D), "
<< "\n T: "
<< (real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]) * setup_.dt) << " s ("
<< real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0])
<< ")"
"\n f: "
<< (real_t(1) / (real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]) * setup_.dt)) << " Hz ("
<< (real_t(1) / real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]))
<< ")"
"\n St (\"real\" D): "
<< strouhalNumberRealD_ << "\n")
}
}
}
std::ofstream outfile( dragFilename_.c_str(), std::ios_base::app );
for(auto it = dragResults.begin(); it != dragResults.end(); ++it){
outfile << it->timestep << ",";
outfile << it->Fx << "," << it->Fy << "," << it->Fz << ",";
outfile << it->cDRealArea << ",";
outfile << it->cLRealArea << ",";
outfile << it->cDDiscreteArea << ",";
outfile << it->cLDiscreteArea;
outfile << "\n";
}
outfile.close();
if (logToFile_ and !dragResults.empty()){
std::ofstream file(filename_.c_str(), std::ios_base::app);
file << coarseExecutionCounter_ - uint_t(1) << " " << force_[0] << " " << force_[1] << " " << force_[2] << " "
<< dragResults.back().cDRealArea << " " << dragResults.back().cLRealArea << " " << dragResults.back().cDDiscreteArea << " " << dragResults.back().cLDiscreteArea << " "
<< pressureDifference_L << " " << pressureDifference << " " << vortexVelocity << " "
<< strouhalNumberRealD_ << '\n';
file.close();
}
dragResults.clear();
}
// WALBERLA_MPI_WORLD_BARRIER();
}
void Evaluation::resetForce()
{
if (!initialized_) refresh();
}
void Evaluation::forceCalculation(const uint_t level)
{
if (rampUpTime_ > coarseExecutionCounter_) return;
if(level == maxLevel_){
for (auto b : finestBlocks_){
force_ += Vector3<double>(boundaryCollection_.ObstacleObject->getForce(b));
}
mpi::reduceInplace(force_, mpi::SUM);
WALBERLA_ROOT_SECTION(){
const double meanU2 = double_c(meanVelocity) * double_c(meanVelocity);
const double cDRealArea = (double_c(8.0) * force_[0]) / (meanU2 * double_c(surfaceAreaSphere));
const double cLRealArea = (double_c(8.0) * force_[1]) / (meanU2 * double_c(surfaceAreaSphere));
const double cDDiscreteArea = (double_c(8.0) * force_[0]) / (meanU2 * double_c(AD_));
const double cLDiscreteArea = (double_c(8.0) * force_[1]) / (meanU2 * double_c(AL_));
DragCoefficient DC(fineExecutionCounter_, force_, cDRealArea, cLRealArea, cDDiscreteArea, cLDiscreteArea);
dragResults.push_back(DC);
fineExecutionCounter_++;
}
}
force_[0] = double_c(0.0);
force_[1] = double_c(0.0);
force_[2] = double_c(0.0);
}
void Evaluation::refresh()
{
WALBERLA_CHECK_NOT_NULLPTR(blocks_)
const uint_t finestLevel = blocks_->getDepth();
real_t AD(real_t(0));
real_t AL(real_t(0));
for (auto block = blocks_->begin(); block != blocks_->end(); ++block){
const uint_t blockLevel = blocks_->getLevel(*block);
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
auto fluid = flagField->getFlag(setup_.fluidUID);
auto obstacle = flagField->getFlag(setup_.obstacleUID);
const real_t area = real_c(4.0);
auto xyzSize = flagField->xyzSize();
for (auto z = xyzSize.zMin(); z <= xyzSize.zMax(); ++z){
for (auto y = xyzSize.yMin(); y <= xyzSize.yMax(); ++y){
for (auto x = xyzSize.xMin(); x <= xyzSize.xMax(); ++x){
if (flagField->isFlagSet(x, y, z, fluid)){
for (auto it = Stencil_T::beginNoCenter(); it != Stencil_T::end(); ++it){
auto nx = x + cell_idx_c(it.cx());
auto ny = y + cell_idx_c(it.cy());
auto nz = z + cell_idx_c(it.cz());
if (flagField->isFlagSet(nx, ny, nz, obstacle)){
WALBERLA_CHECK(blockLevel == finestLevel, "The sphere must be completely located on the finest level")
if (it.cx() == 1 && it.cy() == 0 && it.cz() == 0) { AD += area; }
else if (it.cx() == 0 && it.cz() == 0){
if (it.cy() == 1){
AL += area;
}
}
}
}
}
}
}
}
}
mpi::reduceInplace(AD, mpi::SUM);
mpi::reduceInplace(AL, mpi::SUM);
WALBERLA_ROOT_SECTION(){
AD_ = AD;
AL_ = AL;
}
// Check if points alpha and omega (required for evaluating the pressure difference) are located in fluid cells
if (setup_.evaluatePressure){
auto block = blocks_->getBlock(setup_.pAlpha);
if (block != nullptr){
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks_->getBlockLocalCell(*block, setup_.pAlpha);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(setup_.fluidUID);
if (!flagField->isFlagSet(cell, fluid)){
WALBERLA_ABORT("Cell for evaluating pressure difference at point alpha " << setup_.pAlpha << " is not a fluid cell!")
}
}
block = blocks_->getBlock(setup_.pOmega);
if (block != nullptr){
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks_->getBlockLocalCell(*block, setup_.pOmega);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(setup_.fluidUID);
if (!flagField->isFlagSet(cell, fluid)){
WALBERLA_ABORT("Cell for evaluating pressure difference at point omega " << setup_.pOmega << " is not a fluid cell!")
}
}
}
// Check if point for evaluating the Strouhal number is located inside a fluid cell
if (setup_.evaluateStrouhal){
auto block = blocks_->getBlock(setup_.pStrouhal);
if (block != nullptr){
const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
const auto cell = blocks_->getBlockLocalCell(*block, setup_.pStrouhal);
WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
const auto fluid = flagField->getFlag(setup_.fluidUID);
if (!flagField->isFlagSet(cell, fluid))
WALBERLA_ABORT("Cell for evaluating the Strouhal number at point " << setup_.pStrouhal << " is not a fluid cell!")
}
}
initialized_ = true;
}
void Evaluation::evaluate(real_t& pressureDifference_L, real_t& pressureDifference)
{
if (!initialized_) refresh();
real_t pAlpha(real_t(0));
real_t pOmega(real_t(0));
WALBERLA_CHECK_NOT_NULLPTR(blocks_)
if (setup_.evaluatePressure){
auto block = blocks_->getBlock(setup_.pAlpha);
if (block != nullptr){
const ScalarField_T* const densityField = block->template getData< ScalarField_T >(ids_.densityField);
const auto cell = blocks_->getBlockLocalCell(*block, setup_.pAlpha);
WALBERLA_ASSERT(densityField->xyzSize().contains(cell))
pAlpha += densityField->get(cell) / real_c(3);
}
block = blocks_->getBlock(setup_.pOmega);
if (block != nullptr){
const ScalarField_T* const densityField = block->template getData< ScalarField_T >(ids_.densityField);
const auto cell = blocks_->getBlockLocalCell(*block, setup_.pOmega);
WALBERLA_ASSERT(densityField->xyzSize().contains(cell))
pOmega += densityField->get(cell) / real_c(3);
}
mpi::reduceInplace(pAlpha, mpi::SUM);
mpi::reduceInplace(pOmega, mpi::SUM);
}
WALBERLA_ROOT_SECTION(){
pressureDifference_L = pAlpha - pOmega;
pressureDifference = (pressureDifference_L * setup_.rho * setup_.dxC * setup_.dxC) / (setup_.dt * setup_.dt);
}
}
}
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file Evaluation.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
# pragma once
#include "core/config/Config.h"
#include "core/math/Sample.h"
#include "core/math/Constants.h"
#include "lbm_generated/field/PdfField.h"
#include "sqlite/SQLite.h"
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <functional>
#include <iostream>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
#include "Types.h"
#include "Setup.h"
#include "FlowAroundSphereInfoHeader.h"
using namespace walberla;
using StorageSpecification_T = lbm::FlowAroundSphereStorageSpecification;
using Stencil_T = StorageSpecification_T::Stencil;
using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
using FlagField_T = FlagField< uint8_t >;
using BoundaryCollection_T = lbm::FlowAroundSphereBoundaryCollection< FlagField_T >;
using VoidFunction = std::function<void ()>;
namespace walberla
{
struct DragCoefficient {
uint_t timestep;
double Fx;
double Fy;
double Fz;
double cDRealArea;
double cLRealArea;
double cDDiscreteArea;
double cLDiscreteArea;
DragCoefficient(uint_t t, Vector3<double> f, double cdR, double clR, double cdD, double clD) : timestep(t), Fx(f[0]), Fy(f[1]), Fz(f[2]), cDRealArea(cdR), cLRealArea(clR), cDDiscreteArea(cdD), cLDiscreteArea(clD) {}
};
class Evaluation
{
public:
Evaluation(std::shared_ptr< StructuredBlockForest >& blocks, const uint_t checkFrequency, const uint_t rampUpTime,
BoundaryCollection_T & boundaryCollection,
const IDs& ids, const Setup& setup,
const bool logToStream = true, const bool logToFile = true,
const std::string& filename = std::string("FlowAroundSphere.txt"))
: blocks_(blocks), checkFrequency_(checkFrequency), rampUpTime_(rampUpTime),
boundaryCollection_(boundaryCollection), ids_(ids), setup_(setup),
logToStream_(logToStream), logToFile_(logToFile), filename_(filename)
{
WALBERLA_ASSERT_NOT_NULLPTR(blocks)
maxLevel_ = blocks->getDepth();
blocks->getBlocks(finestBlocks_, maxLevel_);
coefficients_.resize(uint_c(4));
coefficientExtremas_.resize(uint_c(4));
const real_t factor = setup_.dxC / setup_.dxF;
diameterSphere = real_c(2.0) * (setup_.sphereRadius * factor);
surfaceAreaSphere = math::pi * diameterSphere * diameterSphere;
meanVelocity = setup_.inflowVelocity; // (real_c(4.0) * setup_.inflowVelocity) / real_c(9.0);
dragFilename_ = std::string("dragSphereRe_") + std::to_string(uint_t(setup_.reynoldsNumber)) + std::string("_meshLevels_") + std::to_string(uint_t(setup_.refinementLevels + 1)) + std::string(".csv");
WALBERLA_ROOT_SECTION()
{
filesystem::path dataFile( dragFilename_.c_str() );
if( filesystem::exists( dataFile ) )
std::remove( dataFile.string().c_str() );
std::ofstream outfile( dragFilename_.c_str() );
outfile << "SEP=," << "\n";
setup_.writeParameterHeader(outfile);
outfile << "timestep," << "Fx," << "Fy," << "Fz," << "cDRealArea," << "cLRealArea," << "cDDiscreteArea," << "cLDiscreteArea" << "\n";
outfile.close();
if (logToFile_)
{
std::ofstream file(filename_.c_str());
file << "# time step [1], force (x) [2], force (y) [3], force (z) [4], "
"cD [5], cL [6], "
"pressure difference (in lattice units) [7], pressure difference (in Pa) [8], vortex velocity (in "
"lattice units) [9], "
"Strouhal number (real D) [10]"
<< '\n';
if (!setup_.evaluatePressure)
file << "# ATTENTION: pressure was not evaluated, pressure difference is set to zero!" << '\n';
if (!setup_.evaluateStrouhal)
file << "# ATTENTION: vortex velocities were not evaluated, Strouhal number is set to zero!"
<< '\n';
file.close();
}
}
refresh();
WALBERLA_LOG_INFO_ON_ROOT(
"Evaluation initialised:"
"\n + Sphere real area: " << surfaceAreaSphere
<< "\n + Sphere real diameter: " << diameterSphere
<< "\n + Sphere discrete area drag coefficient: " << AD_
<< "\n + Sphere discrete area lift coefficient: " << AL_
)
}
void operator()();
void forceCalculation(const uint_t level); // for calculating the force
void resetForce();
std::function<void (const uint_t)> forceCalculationFunctor()
{
return [this](uint_t level) { forceCalculation(level); };
}
std::function<void()> resetForceFunctor()
{
return [this]() { resetForce(); };
}
void refresh();
protected:
void evaluate(real_t& pressureDifference_L, real_t& pressureDifference);
bool initialized_{false};
std::shared_ptr< StructuredBlockForest > blocks_;
uint_t maxLevel_;
std::vector<Block *> finestBlocks_;
uint_t coarseExecutionCounter_{ uint_c(0) };
uint_t fineExecutionCounter_{ uint_c(0) };
uint_t checkFrequency_;
uint_t rampUpTime_;
BoundaryCollection_T & boundaryCollection_;
IDs ids_;
Setup setup_;
real_t diameterSphere;
real_t surfaceAreaSphere;
real_t meanVelocity;
Vector3< double > force_;
real_t AD_;
real_t AL_;
std::vector<DragCoefficient> dragResults;
std::vector< std::deque< real_t > > coefficients_;
std::vector< std::pair< real_t, real_t > > coefficientExtremas_;
std::vector< real_t > strouhalVelocities_;
std::vector< uint_t > strouhalTimeStep_;
bool strouhalRising_{false};
real_t strouhalNumberRealD_ { real_c(0.0) };
uint_t strouhalEvaluationExecutionCount_ { uint_c(0) };
bool logToStream_;
bool logToFile_;
std::string filename_;
std::string dragFilename_;
std::map< std::string, double > sqlResults_;
}; // class Evaluation
}
\ No newline at end of file
//======================================================================================================================
//
// 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/>.
//
//! \file FlowAroundSphere.cpp
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "blockforest/StructuredBlockForest.h"
#include "blockforest/communication/NonUniformBufferedScheme.h"
#include "core/Abort.h"
#include "core/DataTypes.h"
#include "core/SharedFunctor.h"
#include "core/debug/CheckFunctions.h"
#include "core/logging/Logging.h"
#include "core/math/Vector3.h"
#include "core/Environment.h"
#include "core/mpi/MPIManager.h"
#include "core/MemoryUsage.h"
#include "core/mpi/Reduce.h"
#include "core/timing/RemainingTimeLogger.h"
#include "core/timing/TimingPool.h"
#include "field/AddToStorage.h"
#include "field/CellCounter.h"
#include "field/FlagField.h"
#include "field/StabilityChecker.h"
#include "field/adaptors/AdaptorCreators.h"
#include "field/iterators/FieldIterator.h"
#include "field/vtk/VTKWriter.h"
#include "field/vtk/FlagFieldCellFilter.h"
#include "geometry/InitBoundaryHandling.h"
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
# include "gpu/AddGPUFieldToStorage.h"
# include "gpu/DeviceSelectMPI.h"
# include "gpu/ErrorChecking.h"
# include "gpu/FieldCopy.h"
# include "gpu/HostFieldAllocator.h"
# include "gpu/ParallelStreams.h"
# include "gpu/communication/BlockNonUniformGPUScheme.h"
# include "gpu/communication/UniformGPUScheme.h"
#endif
#include "lbm_generated/communication/NonuniformGeneratedPdfPackInfo.h"
#include "lbm_generated/communication/UniformGeneratedPdfPackInfo.h"
#include "lbm_generated/evaluation/PerformanceEvaluation.h"
#include "lbm_generated/field/AddToStorage.h"
#include "lbm_generated/field/PdfField.h"
#include "lbm_generated/refinement/BasicRecursiveTimeStep.h"
#include "lbm_generated/refinement/RefinementScaling.h"
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
# include "lbm_generated/gpu/AddToStorage.h"
# include "lbm_generated/gpu/BasicRecursiveTimeStepGPU.h"
# include "lbm_generated/gpu/GPUPdfField.h"
# include "lbm_generated/gpu/NonuniformGeneratedGPUPdfPackInfo.h"
# include "lbm_generated/gpu/UniformGeneratedGPUPdfPackInfo.h"
#endif
#include "timeloop/SweepTimeloop.h"
#include "vtk/VTKOutput.h"
#include <cstdlib>
#include <iostream>
#include <memory>
#include "Evaluation.h"
#include "GridGeneration.h"
#include "Setup.h"
#include "Sphere.h"
#include "Types.h"
#include "FlowAroundSphereStaticDefines.h"
using namespace walberla;
using BoundaryCollection_T = lbm::FlowAroundSphereBoundaryCollection< FlagField_T >;
using SweepCollection_T = lbm::FlowAroundSphereSweepCollection;
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
using GPUPdfField_T = lbm_generated::GPUPdfField< StorageSpecification_T >;
using gpu::communication::BlockNonUniformGPUScheme;
using gpu::communication::UniformGPUScheme;
using lbm_generated::NonuniformGeneratedGPUPdfPackInfo;
using lbm_generated::UniformGeneratedGPUPdfPackInfo;
#else
using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
using blockforest::communication::NonUniformBufferedScheme;
using blockforest::communication::UniformBufferedScheme;
using lbm_generated::NonuniformGeneratedPdfPackInfo;
using lbm_generated::UniformGeneratedPdfPackInfo;
#endif
int main(int argc, char** argv)
{
Environment env( argc, argv );
mpi::MPIManager::instance()->useWorldComm();
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::selectDeviceBasedOnMpiRank();
#endif
auto config = env.config();
///////////////////////
/// PARAMETER INPUT ///
///////////////////////
// read general simulation parameters
auto parameters = config->getOneBlock("Parameters");
auto domainParameters = config->getOneBlock("DomainSetup");
auto boundariesConfig = config->getBlock("Boundaries");
auto loggingParameters= config->getOneBlock("Logging");
Setup setup(parameters, domainParameters, loggingParameters, boundariesConfig, infoMap);
bool writeSetupForestAndReturn = loggingParameters.getParameter< bool >("writeSetupForestAndReturn", false);
if (uint_c(MPIManager::instance()->numProcesses()) > 1) writeSetupForestAndReturn = false;
WALBERLA_LOG_INFO_ON_ROOT("Diameter of the Sphere is resolved with " << setup.resolutionSphere << " lattice cells.")
Sphere Sphere( setup );
///////////////////////////
/// CREATE BLOCK FOREST ///
///////////////////////////
shared_ptr< BlockForest > bfs;
createBlockForest(bfs, setup);
if (writeSetupForestAndReturn && mpi::MPIManager::instance()->numProcesses() == 1)
{
WALBERLA_LOG_INFO_ON_ROOT("BlockForest has been created and writen to file. Returning program")
const uint_t nBlocks = bfs->getNumberOfBlocks();
const uint_t nCells = nBlocks * (setup.cellsPerBlock[0] * setup.cellsPerBlock[1] * setup.cellsPerBlock[2]);
const memory_t totalMemory = memory_t(nCells) * setup.memoryPerCell;
WALBERLA_LOG_INFO_ON_ROOT( "Benchmark run data:"
"\n- simulation parameters:" << std::setprecision(16) <<
"\n + collision model: " << infoMap["collisionOperator"] <<
"\n + stencil: " << infoMap["stencil"] <<
"\n + streaming: " << infoMap["streamingPattern"] <<
"\n + compressible: " << ( StorageSpecification_T::compressible ? "yes" : "no" ) <<
"\n + mesh levels: " << setup.refinementLevels + uint_c(1) <<
"\n + resolution: " << setup.dxC << " - on the coarsest grid" <<
"\n + resolution: " << setup.dxF << " - on the finest grid" <<
"\n- domain properties: "
"\n + # blocks: " << nBlocks <<
"\n + # cells: " << uint_c(real_c(nCells) / real_c(1e6)) << " [1e6]"
"\n + # cells sphere D: " << setup.resolutionSphere <<
"\n + total memory: " << totalMemory / 1e9 << " [GB]" <<
"\n- simulation properties: "
"\n + sphere pos.(x): " << setup.sphereXPosition * setup.dxC << " [m]" <<
"\n + sphere pos.(y): " << setup.sphereYPosition * setup.dxC << " [m]" <<
"\n + sphere pos.(z): " << setup.sphereZPosition * setup.dxC << " [m]" <<
"\n + sphere radius: " << setup.sphereRadius * setup.dxC << " [m]" <<
"\n + kin. viscosity: " << setup.viscosity * setup.dxC * setup.dxC / setup.dt << " [m^2/s] (on the coarsest grid)" <<
"\n + viscosity LB: " << setup.viscosity << " [dx*dx/dt] (on the coarsest grid)" <<
"\n + omega: " << setup.omega << " (on the coarsest grid)" <<
"\n + rho: " << setup.rho << " [kg/m^3]" <<
"\n + inflow velocity: " << setup.referenceVelocity << " [m/s]" <<
"\n + lattice velocity: " << setup.inflowVelocity << " [dx/dt]" <<
"\n + Reynolds number: " << setup.reynoldsNumber <<
"\n + Mach number: " << setup.machNumber <<
"\n + dx (coarsest grid): " << setup.dxC << " [m]" <<
"\n + dt (coarsest grid): " << setup.dt << " [s]"
"\n + #time steps: " << setup.timesteps << " (on the coarsest grid, " << ( real_c(1.0) / setup.dt ) << " for 1s of real time)" <<
"\n + simulation time: " << ( real_c( setup.timesteps ) * setup.dt ) << " [s]" )
logging::Logging::printFooterOnStream();
return EXIT_SUCCESS;
}
auto blocks = std::make_shared< StructuredBlockForest >(bfs, setup.cellsPerBlock[0], setup.cellsPerBlock[1], setup.cellsPerBlock[2]);
blocks->createCellBoundingBoxes();
////////////////////////////////////
/// CREATE AND INITIALIZE FIELDS ///
////////////////////////////////////
// create fields
const StorageSpecification_T StorageSpec = StorageSpecification_T();
IDs ids;
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
auto PDFAllocator = make_shared< gpu::HostFieldAllocator< PdfField_T::value_type > >();
auto allocator = make_shared< gpu::HostFieldAllocator< real_t > >();
ids.pdfField = lbm_generated::addPdfFieldToStorage(blocks, "pdfs", StorageSpec, uint_c(2), field::fzyx, PDFAllocator);
ids.velocityField = field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx, uint_c(2), allocator);
ids.densityField = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx, uint_c(2), allocator);
ids.omegaField = field::addToStorage< ScalarField_T >(blocks, "omega", real_c(0.0), field::fzyx, uint_c(2), allocator);
ids.flagField = field::addFlagFieldToStorage< FlagField_T >(blocks, "Boundary Flag Field", uint_c(3));
ids.pdfFieldGPU = lbm_generated::addGPUPdfFieldToStorage< PdfField_T >(blocks, ids.pdfField, StorageSpec, "pdfs on GPU", true);
ids.velocityFieldGPU = gpu::addGPUFieldToStorage< VelocityField_T >(blocks, ids.velocityField, "velocity on GPU", true);
ids.densityFieldGPU = gpu::addGPUFieldToStorage< ScalarField_T >(blocks, ids.densityField, "density on GPU", true);
ids.omegaFieldGPU = gpu::addGPUFieldToStorage< ScalarField_T >(blocks, ids.omegaField, "omega on GPU", true);
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#else
ids.pdfField = lbm_generated::addPdfFieldToStorage(blocks, "pdfs", StorageSpec, uint_c(2), field::fzyx);
ids.velocityField = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx, uint_c(2));
ids.densityField = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx, uint_c(2));
ids.omegaField = field::addToStorage< ScalarField_T >(blocks, "omega", real_c(0.0), field::fzyx, uint_c(2));
ids.flagField = field::addFlagFieldToStorage< FlagField_T >(blocks, "Boundary Flag Field", uint_c(3));
#endif
auto spongeLayer = config->getOneBlock("SpongeLayer");
const bool deactivateSpongeLayer = spongeLayer.getParameter< bool >("deactivateSpongeLayer");
const real_t targetOmega = spongeLayer.getParameter< real_t >("targetOmega");
const real_t spongeStart = spongeLayer.getParameter< real_t >("spongeStart");
const real_t omegaDiff = setup.omega - targetOmega;
const real_t spongeWidth = real_c(blocks->getDomain().xMax()) - spongeStart;
const real_t spongeMid = spongeStart + (spongeWidth / real_c(2.0));
if(!deactivateSpongeLayer)
WALBERLA_LOG_WARNING_ON_ROOT("Using a sponge layer at the Outlet. The sponge layer starts at " << spongeStart << " [m] and targets a relaxation rate of " << targetOmega << " at the outlet" )
for (auto& block : *blocks)
{
Block& b = dynamic_cast< Block& >(block);
const uint_t level = blocks->getLevel(b);
auto * omegaField = b.getData< ScalarField_T > ( ids.omegaField );
for( auto it = omegaField->beginWithGhostLayer(0); it != omegaField->end(); ++it ){
if(deactivateSpongeLayer){
omegaField->get(it.cell()) = lbm_generated::relaxationRateScaling(setup.omega, level);
}
else{
Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, it.cell());
Vector3<real_t> cellCenter = blocks->getCellCenter(globalCell, level);
const real_t factor = real_c(0.5) + real_c(0.5) * std::tanh(real_c(2.0) * (cellCenter[0] - spongeMid) / spongeWidth);
omegaField->get(it.cell()) = lbm_generated::relaxationRateScaling(setup.omega - (factor * omegaDiff), level);
}
}
}
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< gpu::GPUField< real_t >, ScalarField_T>(blocks, ids.omegaFieldGPU, ids.omegaField);
#endif
WALBERLA_MPI_BARRIER()
const Cell innerOuterSplit =
Cell(parameters.getParameter< Vector3< cell_idx_t > >("innerOuterSplit", Vector3< cell_idx_t >(1, 1, 1)));
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
SweepCollection_T sweepCollection(blocks, ids.pdfFieldGPU, ids.omegaFieldGPU, ids.densityFieldGPU, ids.velocityFieldGPU, innerOuterSplit);
#else
SweepCollection_T sweepCollection(blocks, ids.pdfField, ids.omegaField, ids.densityField, ids.velocityField, innerOuterSplit);
#endif
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
WALBERLA_LOG_INFO_ON_ROOT("Setting up communication...")
const bool gpuEnabledMPI = parameters.getParameter< bool >("gpuEnabledMPI", false);
const bool asyncCommunication = parameters.getParameter< bool >("asyncCommunication");
auto nonUniformCommunication = std::make_shared< BlockNonUniformGPUScheme< CommunicationStencil_T > >(blocks, gpuEnabledMPI);
auto nonUniformPackInfo = lbm_generated::setupNonuniformGPUPdfCommunication< GPUPdfField_T >(blocks, ids.pdfFieldGPU);
nonUniformCommunication->addPackInfo(nonUniformPackInfo);
if(asyncCommunication) { nonUniformCommunication->activateGPUStreams(); }
#else
WALBERLA_LOG_INFO_ON_ROOT("Setting up communication...")
auto nonUniformCommunication = std::make_shared< NonUniformBufferedScheme< CommunicationStencil_T > >(blocks);
auto nonUniformPackInfo = lbm_generated::setupNonuniformPdfCommunication< PdfField_T >(blocks, ids.pdfField);
nonUniformCommunication->addPackInfo(nonUniformPackInfo);
#endif
WALBERLA_MPI_BARRIER()
WALBERLA_LOG_INFO_ON_ROOT("Setting up communication done")
/////////////////////////
/// BOUNDARY HANDLING ///
/////////////////////////
WALBERLA_LOG_INFO_ON_ROOT("Start BOUNDARY HANDLING")
// create and initialize boundary handling
Sphere.setupBoundary(blocks, ids.flagField);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, setup.fluidUID, cell_idx_c(0));
if(parameters.getParameter< bool >("initialiseWithInletVelocity"))
{
for (auto& block : *blocks)
{
auto * flagField = block.getData< FlagField_T > ( ids.flagField );
auto * velField = block.getData< VelocityField_T > ( ids.velocityField );
// auto domainFlag = flagField->getFlag(fluidFlagUID);
for( auto it = flagField->beginWithGhostLayer(2); it != flagField->end(); ++it )
{
// if (!isFlagSet(it, domainFlag))
// continue;
velField->get(it.cell(), 0) = setup.inflowVelocity;
}
}
}
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< gpu::GPUField< real_t >, VelocityField_T>(blocks, ids.velocityFieldGPU, ids.velocityField);
#endif
for (auto& block : *blocks)
{
sweepCollection.initialise(&block, cell_idx_c(1));
}
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< PdfField_T, gpu::GPUField< PdfField_T::value_type > >(blocks, ids.pdfField, ids.pdfFieldGPU);
#endif
sweepCollection.initialiseBlockPointer();
std::function< real_t(const Cell&, const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
wallDistanceFunctor = wallDistance(Sphere);
const real_t omegaFinestLevel = lbm_generated::relaxationRateScaling(setup.omega, setup.refinementLevels);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, setup.fluidUID, omegaFinestLevel, setup.inflowVelocity, wallDistanceFunctor, ids.pdfField);
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#else
BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfField, setup.fluidUID, omegaFinestLevel, setup.inflowVelocity, wallDistanceFunctor);
#endif
WALBERLA_MPI_BARRIER()
WALBERLA_LOG_INFO_ON_ROOT("BOUNDARY HANDLING done")
//////////////////////////////////
/// SET UP SWEEPS AND TIMELOOP ///
//////////////////////////////////
WALBERLA_LOG_INFO_ON_ROOT("Start SWEEPS AND TIMELOOP")
// flow evaluation
auto EvaluationParameters = config->getOneBlock("Evaluation");
const uint_t evaluationCheckFrequency = EvaluationParameters.getParameter< uint_t >("evaluationCheckFrequency");
const uint_t rampUpTime = EvaluationParameters.getParameter< uint_t >("rampUpTime");
const bool evaluationLogToStream = EvaluationParameters.getParameter< bool >("logToStream");
const bool evaluationLogToFile = EvaluationParameters.getParameter< bool >("logToFile");
const std::string evaluationFilename = EvaluationParameters.getParameter< std::string >("filename");
shared_ptr< Evaluation > evaluation( new Evaluation( blocks, evaluationCheckFrequency, rampUpTime, boundaryCollection,
ids, setup, evaluationLogToStream, evaluationLogToFile, evaluationFilename));
// create time loop
SweepTimeloop timeloop(blocks->getBlockStorage(), setup.timesteps);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
int streamHighPriority = 0;
int streamLowPriority = 0;
WALBERLA_GPU_CHECK(gpuDeviceGetStreamPriorityRange(&streamLowPriority, &streamHighPriority))
sweepCollection.setOuterPriority(streamHighPriority);
auto defaultStream = gpu::StreamRAII::newPriorityStream(streamLowPriority);
lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollection_T, BoundaryCollection_T >
LBMMeshRefinement(blocks, ids.pdfFieldGPU, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
LBMMeshRefinement.addRefinementToTimeLoop(timeloop, uint_c(0));
LBMMeshRefinement.addPostBoundaryHandlingBlockFunction(evaluation->forceCalculationFunctor());
#else
std::shared_ptr< lbm_generated::BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T > >
LBMRefinement;
LBMRefinement = std::make_shared<
lbm_generated::BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T > >(
blocks, ids.pdfField, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
LBMRefinement->addPostBoundaryHandlingBlockFunction(evaluation->forceCalculationFunctor());
LBMRefinement->addRefinementToTimeLoop(timeloop);
#endif
//////////////////
/// VTK OUTPUT ///
//////////////////
WALBERLA_LOG_INFO_ON_ROOT("SWEEPS AND TIMELOOP done")
auto VTKWriter = config->getOneBlock("VTKWriter");
const uint_t vtkWriteFrequency = VTKWriter.getParameter< uint_t >("vtkWriteFrequency", 0);
const uint_t vtkGhostLayers = VTKWriter.getParameter< uint_t >("vtkGhostLayers", 0);
const bool amrFileFormat = VTKWriter.getParameter< bool >("amrFileFormat", false);
const bool oneFilePerProcess = VTKWriter.getParameter< bool >("oneFilePerProcess", false);
const real_t samplingResolution = VTKWriter.getParameter< real_t >("samplingResolution", real_c(-1.0));
const uint_t initialWriteCallsToSkip = VTKWriter.getParameter< uint_t >("initialWriteCallsToSkip", uint_c(0.0));
const real_t sliceThickness = VTKWriter.getParameter< real_t >("sliceThickness", real_c(1.0));
auto finalDomain = blocks->getDomain();
if (vtkWriteFrequency > 0)
{
const std::string vtkFolder = "VTKSphereRE_" + std::to_string(uint64_c(setup.reynoldsNumber));
auto vtkOutput =
vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, vtkGhostLayers, false, vtkFolder,
"simulation_step", false, true, true, false, 0, amrFileFormat, oneFilePerProcess);
vtkOutput->setInitialWriteCallsToSkip(initialWriteCallsToSkip);
vtkOutput->addBeforeFunction([&]() {
for (auto& block : *blocks)
{
sweepCollection.calculateMacroscopicParameters(&block);
}
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< VelocityField_T, gpu::GPUField< real_t > >(blocks, ids.velocityField, ids.velocityFieldGPU);
gpu::fieldCpy< ScalarField_T, gpu::GPUField< real_t > >(blocks, ids.densityField, ids.densityFieldGPU);
#endif
});
vtkOutput->setSamplingResolution(samplingResolution);
field::FlagFieldCellFilter<FlagField_T> fluidFilter( ids.flagField );
fluidFilter.addFlag( setup.obstacleUID );
vtkOutput->addCellExclusionFilter(fluidFilter);
if (VTKWriter.getParameter< bool >("writeOnlySlice", true)){
const AABB sliceXY(finalDomain.xMin(), finalDomain.yMin(), finalDomain.center()[2] - sliceThickness * setup.dxF,
finalDomain.xMax(), finalDomain.yMax(), finalDomain.center()[2] + sliceThickness * setup.dxF);
vtkOutput->addCellInclusionFilter(vtk::AABBCellFilter(sliceXY));
if (VTKWriter.getParameter< bool >("writeXZSlice", false)){
const AABB sliceXZ(finalDomain.xMin(), finalDomain.center()[1] - setup.dxF, finalDomain.zMin(),
finalDomain.xMax(), finalDomain.center()[1] + setup.dxF, finalDomain.zMax());
vtkOutput->addCellInclusionFilter(vtk::AABBCellFilter(sliceXZ));
}
}
if (VTKWriter.getParameter< bool >("velocity"))
{
auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(ids.velocityField, "velocity");
vtkOutput->addCellDataWriter(velWriter);
}
if (VTKWriter.getParameter< bool >("density"))
{
auto densityWriter = make_shared< field::VTKWriter< ScalarField_T, float32 > >(ids.densityField, "density");
vtkOutput->addCellDataWriter(densityWriter);
}
if (VTKWriter.getParameter< bool >("omega"))
{
auto omegaWriter = make_shared< field::VTKWriter< ScalarField_T, float32 > >(ids.omegaField, "omega");
vtkOutput->addCellDataWriter(omegaWriter);
}
if (VTKWriter.getParameter< bool >("flag"))
{
auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(ids.flagField, "flag");
vtkOutput->addCellDataWriter(flagWriter);
}
timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
}
// log remaining time
const real_t remainingTimeLoggerFrequency =
loggingParameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0); // in seconds
if (uint_c(remainingTimeLoggerFrequency) > 0)
{
timeloop.addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
"remaining time logger");
}
// LBM stability check
auto CheckerParameters = config->getOneBlock("StabilityChecker");
const uint_t checkFrequency = CheckerParameters.getParameter< uint_t >("checkFrequency", uint_t(0));
if (checkFrequency > 0)
{
auto checkFunction = [](PdfField_T::value_type value) { return value < math::abs(PdfField_T::value_type(10)); };
timeloop.addFuncAfterTimeStep(
makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
config, blocks, ids.pdfField, ids.flagField, setup.fluidUID, checkFunction)),
"Stability check");
}
timeloop.addFuncBeforeTimeStep( SharedFunctor< Evaluation >(evaluation), "evaluation" );
WALBERLA_MPI_BARRIER()
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif
// WALBERLA_LOG_INFO_ON_ROOT("Execute single timestep to fully complete the preprocessing done")
//////////////////////
/// RUN SIMULATION ///
//////////////////////
const lbm_generated::PerformanceEvaluation< FlagField_T > performance(blocks, ids.flagField, setup.fluidUID);
field::CellCounter< FlagField_T > fluidCells(blocks, ids.flagField, setup.fluidUID);
fluidCells();
WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << blocks->getNumberOfBlocks())
for (uint_t level = 0; level <= setup.refinementLevels; level++)
{
WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << blocks->getNumberOfBlocks(level))
}
WALBERLA_LOG_INFO_ON_ROOT( "Benchmark run data:"
"\n- simulation parameters:" << std::setprecision(16) <<
"\n + collision model: " << infoMap["collisionOperator"] <<
"\n + stencil: " << infoMap["stencil"] <<
"\n + streaming: " << infoMap["streamingPattern"] <<
"\n + compressible: " << ( StorageSpecification_T::compressible ? "yes" : "no" ) <<
"\n + mesh levels: " << setup.refinementLevels + uint_c(1) <<
"\n + resolution: " << setup.dxC << " - on the coarsest grid" <<
"\n + resolution: " << setup.dxF << " - on the finest grid" <<
"\n- simulation properties: "
"\n + sphere pos.(x): " << setup.sphereXPosition * setup.dxC << " [m]" <<
"\n + sphere pos.(y): " << setup.sphereYPosition * setup.dxC << " [m]" <<
"\n + sphere pos.(z): " << setup.sphereZPosition * setup.dxC << " [m]" <<
"\n + sphere radius: " << setup.sphereRadius * setup.dxC << " [m]" <<
"\n + kin. viscosity: " << setup.viscosity * setup.dxC * setup.dxC / setup.dt << " [m^2/s] (on the coarsest grid)" <<
"\n + viscosity LB: " << setup.viscosity << " [dx*dx/dt] (on the coarsest grid)" <<
"\n + omega: " << setup.omega << " (on the coarsest grid)" <<
"\n + rho: " << setup.rho << " [kg/m^3]" <<
"\n + inflow velocity: " << setup.referenceVelocity << " [m/s]" <<
"\n + lattice velocity: " << setup.inflowVelocity << " [dx/dt]" <<
"\n + Reynolds number: " << setup.reynoldsNumber <<
"\n + Mach number: " << setup.machNumber <<
"\n + dt (coarsest grid): " << setup.dt << " [s]"
"\n + #time steps: " << timeloop.getNrOfTimeSteps() << " (on the coarsest grid, " << ( real_t(1) / setup.dt ) << " for 1s of real time)" <<
"\n + simulation time: " << ( real_c( timeloop.getNrOfTimeSteps() ) * setup.dt ) << " [s]" )
WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation")
WcTimingPool timeloopTiming;
WcTimer simTimer;
WALBERLA_MPI_BARRIER()
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif
simTimer.start();
timeloop.run(timeloopTiming);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
LBMMeshRefinement.waitAllCommunications();
WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif
WALBERLA_MPI_BARRIER()
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
real_t time = simTimer.max();
WALBERLA_MPI_SECTION() { walberla::mpi::reduceInplace(time, walberla::mpi::MAX); }
performance.logResultOnRoot(setup.timesteps, time);
const auto reducedTimeloopTiming = timeloopTiming.getReduced();
WALBERLA_LOG_RESULT_ON_ROOT("Time loop timing:\n" << *reducedTimeloopTiming)
printResidentMemoryStatistics();
return EXIT_SUCCESS;
}
\ No newline at end of file
Parameters
{
reynoldsNumber 1000000;
diameterSphere 1.0;
sphereXPosition 12;
referenceVelocity 1.0;
latticeVelocity 0.05;
initialiseWithInletVelocity true;
// coarseMeshSize 0.0078125;
// coarseMeshSize 0.015625;
// coarseMeshSize 0.03125;
coarseMeshSize 0.0625;
evaluatePressure false;
pAlpha < 0.45, 0.2, 0.205 >; // points for evaluating
pOmega < 0.55, 0.2, 0.205 >; // the pressure difference
evaluateStrouhal false;
pStrouhal < 1, 0.325, 0.205 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
timesteps 60001;
// timesteps 101;
processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>;
// GPU related Parameters, only used if GPU is enabled
gpuEnabledMPI true;
asyncCommunication false;
}
//! [domainSetup]
DomainSetup
{
// cellsPerBlock < 130, 240, 240 >;
// domainSize < 32.5, 15, 15 >;
cellsPerBlock < 64, 64, 64 >;
domainSize < 40, 20, 20 >;
periodic < false, false, false >;
refinementLevels 5;
numberProcesses 1; // This is for load balancing, overwritten if more than one proc is used
blockForestFilestem flowAroundSphereBlockForest;
}
//! [domainSetup]
Boundaries
{
sphere Obstacle;
inflow UBB;
outflow FixedDensity;
walls FreeSlip;
}
SpongeLayer
{
deactivateSpongeLayer false;
targetOmega 1.9;
spongeStart 36;
}
StabilityChecker
{
checkFrequency 0;
streamOutput false;
vtkOutput true;
}
VTKWriter
{
vtkWriteFrequency 5000;
vtkGhostLayers 0;
velocity true;
density true;
flag false;
omega false;
writeOnlySlice true;
sliceThickness 1;
writeXZSlice false;
amrFileFormat true;
oneFilePerProcess false;
samplingResolution -1;
initialWriteCallsToSkip 0;
}
Logging
{
logLevel info; // info progress detail tracing
writeSetupForestAndReturn true;
readSetupFromFile false;
remainingTimeLoggerFrequency 600; // in seconds
}
Evaluation
{
evaluationCheckFrequency 500;
rampUpTime 0;
logToStream true;
logToFile false;
filename FlowAroundSphere.txt;
}
import sympy as sp
from pystencils import Target
from pystencils.field import fields
from pystencils.simp import insert_aliases, insert_constants
from lbmpy import LBStencil, LBMConfig, LBMOptimisation
from lbmpy.boundaries.boundaryconditions import (ExtrapolationOutflow, UBB, QuadraticBounceBack,
FreeSlip, NoSlip, FixedDensity)
from lbmpy.creationfunctions import create_lb_collision_rule
from lbmpy.enums import Method, Stencil, SubgridScaleModel
from pystencils_walberla import CodeGeneration, generate_info_header
from lbmpy_walberla import generate_lbm_package, lbm_boundary_generator
info_header = """
#pragma once
#include <map>
std::map<std::string, std::string> infoMap{{{{"stencil", "{stencil}"}},
{{"streamingPattern", "{streaming_pattern}"}},
{{"collisionOperator", "{collision_operator}"}}}};
"""
omega = sp.symbols("omega")
inlet_velocity = sp.symbols("u_x")
max_threads = 256
with CodeGeneration() as ctx:
target = Target.GPU if ctx.gpu else Target.CPU
sweep_params = {'block_size': (128, 1, 1)} if ctx.gpu else {}
dtype = 'float64'
pdf_dtype = 'float64'
stencil = LBStencil(Stencil.D3Q27)
q = stencil.Q
dim = stencil.D
streaming_pattern = 'esopull'
pdfs = fields(f"pdfs_src({stencil.Q}): {pdf_dtype}[3D]", layout='fzyx')
velocity_field, density_field = fields(f"velocity({dim}), density(1) : {dtype}[{dim}D]", layout='fzyx')
omega_field = fields(f"rr(1) : {dtype}[{dim}D]", layout='fzyx')
macroscopic_fields = {'density': density_field, 'velocity': velocity_field}
method_enum = Method.CUMULANT
lbm_config = LBMConfig(
method=method_enum,
stencil=stencil,
relaxation_rate=omega_field.center,
compressible=True,
# subgrid_scale_model=SubgridScaleModel.QR,
fourth_order_correction=0.01 if method_enum == Method.CUMULANT and stencil.Q == 27 else False,
field_name='pdfs',
streaming_pattern=streaming_pattern,
)
lbm_opt = LBMOptimisation(cse_global=False, cse_pdfs=False, field_layout="fzyx",
symbolic_field=pdfs)
collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
if lbm_config.method == Method.CUMULANT:
collision_rule = insert_constants(collision_rule)
collision_rule = insert_aliases(collision_rule)
lb_method = collision_rule.method
freeslip = lbm_boundary_generator("FreeSlip", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil),
field_data_type=pdf_dtype)
no_slip = lbm_boundary_generator(class_name='NoSlip', flag_uid='NoSlip',
boundary_object=NoSlip(), field_data_type=pdf_dtype)
quadratic_bounce_back = QuadraticBounceBack(omega, calculate_force_on_boundary=True)
no_slip_interpolated = lbm_boundary_generator(class_name='Obstacle', flag_uid='Obstacle',
boundary_object=quadratic_bounce_back,
field_data_type=pdf_dtype)
ubb = lbm_boundary_generator(class_name='UBB', flag_uid='UBB',
boundary_object=UBB((inlet_velocity, 0.0, 0.0), density=1.0, data_type=dtype),
field_data_type=pdf_dtype)
outflow_boundary = ExtrapolationOutflow(stencil[4], lb_method)
outflow = lbm_boundary_generator(class_name='Outflow', flag_uid='Outflow',
boundary_object=ExtrapolationOutflow(stencil[4], lb_method),
field_data_type=pdf_dtype)
fixed_density = lbm_boundary_generator(class_name='FixedDensity', flag_uid='FixedDensity',
boundary_object=FixedDensity(1.0),
field_data_type=pdf_dtype)
generate_lbm_package(ctx, name="FlowAroundSphere", collision_rule=collision_rule,
lbm_config=lbm_config, lbm_optimisation=lbm_opt,
nonuniform=True, boundaries=[freeslip, no_slip, no_slip_interpolated,
ubb, outflow, fixed_density],
macroscopic_fields=macroscopic_fields, gpu_indexing_params=sweep_params,
target=target, data_type=dtype, pdfs_data_type=pdf_dtype,
max_threads=max_threads)
field_typedefs = {'VelocityField_T': velocity_field,
'ScalarField_T': density_field}
# Info header containing correct template definitions for stencil and field
generate_info_header(ctx, 'FlowAroundSphereInfoHeader',
field_typedefs=field_typedefs)
infoHeaderParams = {
'stencil': stencil.name.lower(),
'streaming_pattern': streaming_pattern,
'collision_operator': lbm_config.method.name.lower(),
}
ctx.write_file("FlowAroundSphereStaticDefines.h", info_header.format(**infoHeaderParams))