Dear CS10-Gitlab-users, on Thursday, Feb 3 there will be maintenance. That will lead to a downtime of the CS10-Gitlab-service including Subversion and Mattermost chat from 09:30. This might take the whole day since we don't know how long it is going to take. We are sorry for the inconvenience! Best regards, CS10-Admin-Team

Commit 139d60a7 authored by Markus Holzer's avatar Markus Holzer
Browse files

Added Flow around sphere test case

parent fd9ad0cd
Pipeline #30145 passed with stages
in 244 minutes and 45 seconds
......@@ -19,8 +19,9 @@ if ( WALBERLA_BUILD_WITH_PYTHON )
add_subdirectory( FieldCommunication )
if ( WALBERLA_BUILD_WITH_CODEGEN )
add_subdirectory( UniformGridGenerated )
add_subdirectory( FlowAroundSphereCodeGen )
add_subdirectory( PhaseFieldAllenCahn )
add_subdirectory( UniformGridGenerated )
endif()
if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_CUDA )
......
waLBerla_link_files_to_builddir( "*.py" )
if (WALBERLA_BUILD_WITH_CUDA)
waLBerla_generate_target_from_python(NAME FlowAroundSphereGenerated
FILE FlowAroundSphereCodeGen.py
OUT_FILES FlowAroundSphereCodeGen_EvenSweep.cu FlowAroundSphereCodeGen_EvenSweep.h
FlowAroundSphereCodeGen_OddSweep.cu FlowAroundSphereCodeGen_OddSweep.h
FlowAroundSphereCodeGen_MacroSetter.cu FlowAroundSphereCodeGen_MacroSetter.h
FlowAroundSphereCodeGen_UBB.cu FlowAroundSphereCodeGen_UBB.h
FlowAroundSphereCodeGen_NoSlip.cu FlowAroundSphereCodeGen_NoSlip.h
FlowAroundSphereCodeGen_Outflow.cu FlowAroundSphereCodeGen_Outflow.h
FlowAroundSphereCodeGen_PackInfoEven.cu FlowAroundSphereCodeGen_PackInfoEven.h
FlowAroundSphereCodeGen_PackInfoOdd.cu FlowAroundSphereCodeGen_PackInfoOdd.h
FlowAroundSphereCodeGen_InfoHeader.h)
waLBerla_add_executable( NAME FlowAroundSphereCodeGen FILE FlowAroundSphereCodeGen.cpp
DEPENDS blockforest boundary core domain_decomposition field geometry timeloop vtk FlowAroundSphereGenerated)
set_target_properties( FlowAroundSphereCodeGen PROPERTIES CXX_VISIBILITY_PRESET hidden)
else ()
waLBerla_generate_target_from_python(NAME FlowAroundSphereGenerated
FILE FlowAroundSphereCodeGen.py
OUT_FILES FlowAroundSphereCodeGen_EvenSweep.cpp FlowAroundSphereCodeGen_EvenSweep.h
FlowAroundSphereCodeGen_OddSweep.cpp FlowAroundSphereCodeGen_OddSweep.h
FlowAroundSphereCodeGen_MacroSetter.cpp FlowAroundSphereCodeGen_MacroSetter.h
FlowAroundSphereCodeGen_UBB.cpp FlowAroundSphereCodeGen_UBB.h
FlowAroundSphereCodeGen_NoSlip.cpp FlowAroundSphereCodeGen_NoSlip.h
FlowAroundSphereCodeGen_Outflow.cpp FlowAroundSphereCodeGen_Outflow.h
FlowAroundSphereCodeGen_PackInfoEven.cpp FlowAroundSphereCodeGen_PackInfoEven.h
FlowAroundSphereCodeGen_PackInfoOdd.cpp FlowAroundSphereCodeGen_PackInfoOdd.h
FlowAroundSphereCodeGen_InfoHeader.h)
waLBerla_add_executable( NAME FlowAroundSphereCodeGen FILE FlowAroundSphereCodeGen.cpp
DEPENDS blockforest boundary core domain_decomposition field geometry timeloop vtk FlowAroundSphereGenerated)
set_target_properties( FlowAroundSphereCodeGen PROPERTIES CXX_VISIBILITY_PRESET hidden)
endif()
\ 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 FlowAroundSphereCodeGen.cpp
//! \author Frederik Hennig <frederik.hennig@fau.de>
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#include "blockforest/all.h"
#include "core/all.h"
#include "domain_decomposition/all.h"
#include "field/all.h"
#include "geometry/all.h"
#include "lbm/vtk/QCriterion.h"
#include "python_coupling/CreateConfig.h"
#include "python_coupling/PythonCallback.h"
#include "timeloop/all.h"
#if defined(WALBERLA_BUILD_WITH_CUDA)
# include "cuda/AddGPUFieldToStorage.h"
# include "cuda/DeviceSelectMPI.h"
# include "cuda/HostFieldAllocator.h"
# include "cuda/NVTX.h"
# include "cuda/ParallelStreams.h"
# include "cuda/communication/GPUPackInfo.h"
# include "cuda/communication/UniformGPUScheme.h"
#endif
// CodeGen includes
#include "FlowAroundSphereCodeGen_InfoHeader.h"
#include "FlowAroundSphereCodeGen_EvenSweep.h"
#include "FlowAroundSphereCodeGen_MacroSetter.h"
#include "FlowAroundSphereCodeGen_NoSlip.h"
#include "FlowAroundSphereCodeGen_OddSweep.h"
#include "FlowAroundSphereCodeGen_Outflow.h"
#include "FlowAroundSphereCodeGen_PackInfoEven.h"
#include "FlowAroundSphereCodeGen_PackInfoOdd.h"
#include "FlowAroundSphereCodeGen_UBB.h"
typedef lbm::FlowAroundSphereCodeGen_PackInfoEven PackInfoEven_T;
typedef lbm::FlowAroundSphereCodeGen_PackInfoOdd PackInfoOdd_T;
typedef walberla::uint8_t flag_t;
typedef FlagField< flag_t > FlagField_T;
#if defined(WALBERLA_BUILD_WITH_CUDA)
typedef cuda::GPUField< real_t > GPUField;
#endif
using namespace std::placeholders;
auto pdfFieldAdder = [](IBlock* const block, StructuredBlockStorage* const storage) {
return new PdfField_T(storage->getNumberOfXCells(*block), storage->getNumberOfYCells(*block),
storage->getNumberOfZCells(*block), uint_t(1), field::fzyx,
make_shared< field::AllocateAligned< real_t, 64 > >());
};
auto VelocityCallback = [](const Cell& pos, const shared_ptr< StructuredBlockForest >& SbF, IBlock& block,
real_t inflow_velocity) {
Cell globalCell;
CellInterval domain = SbF->getDomainCellBB();
real_t h_y = domain.yMax() - domain.yMin();
real_t h_z = domain.zMax() - domain.zMin();
SbF->transformBlockLocalToGlobalCell(globalCell, block, pos);
real_t y1 = globalCell[1] - (h_y / 2.0 + 0.5);
real_t z1 = globalCell[2] - (h_z / 2.0 + 0.5);
real_t u = (inflow_velocity * 16) / (h_y * h_y * h_z * h_z) * (h_y / 2.0 - y1) * (h_y / 2 + y1) * (h_z / 2 - z1) *
(h_z / 2 + z1);
Vector3< real_t > result(u, 0.0, 0.0);
return result;
};
class TimestepModulusTracker
{
private:
uint_t modulus_;
public:
TimestepModulusTracker(uint_t initialTimestep) : modulus_(initialTimestep & 1){};
void setTimestep(uint_t timestep) { modulus_ = timestep & 1; }
std::function< void() > advancementFunction()
{
return [this]() { this->modulus_ = (this->modulus_ + 1) & 1; };
}
uint_t modulus() const { return modulus_; }
bool evenStep() const { return static_cast< bool >(modulus_ ^ 1); }
bool oddStep() const { return static_cast< bool >(modulus_ & 1); }
};
class AlternatingSweep
{
public:
typedef std::function< void(IBlock*) > SweepFunction;
AlternatingSweep(SweepFunction evenSweep, SweepFunction oddSweep, std::shared_ptr< TimestepModulusTracker > tracker)
: tracker_(tracker), sweeps_{ evenSweep, oddSweep } {};
void operator()(IBlock* block) { sweeps_[tracker_->modulus()](block); }
private:
std::shared_ptr< TimestepModulusTracker > tracker_;
std::vector< SweepFunction > sweeps_;
};
class AlternatingBeforeFunction
{
public:
typedef std::function< void() > BeforeFunction;
AlternatingBeforeFunction(BeforeFunction evenFunc, BeforeFunction oddFunc,
std::shared_ptr< TimestepModulusTracker >& tracker)
: tracker_(tracker), funcs_{ evenFunc, oddFunc } {};
void operator()() { funcs_[tracker_->modulus()](); }
private:
std::shared_ptr< TimestepModulusTracker > tracker_;
std::vector< BeforeFunction > funcs_;
};
class Filter
{
public:
explicit Filter(Vector3< uint_t > numberOfCells) : numberOfCells_(numberOfCells) {}
void operator()(const IBlock& /*block*/) {}
bool operator()(const cell_idx_t x, const cell_idx_t y, const cell_idx_t z) const
{
return x >= -1 && x <= cell_idx_t(numberOfCells_[0]) && y >= -1 && y <= cell_idx_t(numberOfCells_[1]) &&
z >= -1 && z <= cell_idx_t(numberOfCells_[2]);
}
private:
Vector3< uint_t > numberOfCells_;
};
using FluidFilter_T = Filter;
int main(int argc, char** argv)
{
walberla::Environment walberlaEnv(argc, argv);
#if defined(WALBERLA_BUILD_WITH_CUDA)
cuda::selectDeviceBasedOnMpiRank();
#endif
for (auto cfg = python_coupling::configBegin(argc, argv); cfg != python_coupling::configEnd(); ++cfg)
{
WALBERLA_MPI_WORLD_BARRIER()
auto config = *cfg;
logging::configureLogging(config);
auto blocks = blockforest::createUniformBlockGridFromConfig(config);
// read parameters
Vector3< uint_t > cellsPerBlock =
config->getBlock("DomainSetup").getParameter< Vector3< uint_t > >("cellsPerBlock");
auto parameters = config->getOneBlock("Parameters");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(10));
const real_t omega = parameters.getParameter< real_t >("omega", real_t(1.9));
const real_t u_max = parameters.getParameter< real_t >("u_max", real_t(0.05));
const real_t reynolds_number = parameters.getParameter< real_t >("reynolds_number", real_t(1000));
const uint_t diameter_sphere = parameters.getParameter< uint_t >("diameter_sphere", uint_t(5));
const double remainingTimeLoggerFrequency =
parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
// create fields
BlockDataID pdfFieldID = blocks->addStructuredBlockData< PdfField_T >(pdfFieldAdder, "PDFs");
BlockDataID velFieldID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_t(0), field::fzyx);
BlockDataID densityFieldID = field::addToStorage< ScalarField_T >(blocks, "density", real_t(0), field::fzyx);
#if defined(WALBERLA_BUILD_WITH_CUDA)
BlockDataID pdfFieldIDGPU = cuda::addGPUFieldToStorage< PdfField_T >(blocks, pdfFieldID, "PDFs on GPU", true);
BlockDataID velFieldIDGPU =
cuda::addGPUFieldToStorage< VelocityField_T >(blocks, velFieldID, "velocity on GPU", true);
BlockDataID densityFieldIDGPU =
cuda::addGPUFieldToStorage< ScalarField_T >(blocks, densityFieldID, "density on GPU", true);
#endif
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
// initialise all PDFs
#if defined(WALBERLA_BUILD_WITH_CUDA)
pystencils::FlowAroundSphereCodeGen_MacroSetter setterSweep(pdfFieldIDGPU, velFieldIDGPU);
for (auto& block : *blocks)
setterSweep(&block);
cuda::fieldCpy< PdfField_T , GPUField >(blocks, pdfFieldID, pdfFieldIDGPU);
#else
pystencils::FlowAroundSphereCodeGen_MacroSetter setterSweep(pdfFieldID, velFieldID);
for (auto& block : *blocks)
setterSweep(&block);
#endif
// Create communication
#if defined(WALBERLA_BUILD_WITH_CUDA)
cuda::communication::UniformGPUScheme< Stencil_T > comEven(blocks, false);
comEven.addPackInfo(make_shared< PackInfoEven_T >(pdfFieldIDGPU));
auto evenComm = std::function< void() >([&]() { comEven.communicate(nullptr); });
cuda::communication::UniformGPUScheme< Stencil_T > comODD(blocks, false);
comODD.addPackInfo(make_shared< PackInfoOdd_T >(pdfFieldIDGPU));
auto oddComm = std::function< void() >([&]() { comODD.communicate(nullptr); });
#else
blockforest::communication::UniformBufferedScheme< Stencil_T > evenComm(blocks);
evenComm.addPackInfo(make_shared< PackInfoEven_T >(pdfFieldID));
blockforest::communication::UniformBufferedScheme< Stencil_T > oddComm(blocks);
oddComm.addPackInfo(make_shared< PackInfoOdd_T >(pdfFieldID));
#endif
// create and initialize boundary handling
const FlagUID fluidFlagUID("Fluid");
auto boundariesConfig = config->getOneBlock("Boundaries");
std::function< Vector3< real_t >(const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
velocity_initialisation = std::bind(VelocityCallback, _1, _2, _3, u_max);
#if defined(WALBERLA_BUILD_WITH_CUDA)
lbm::FlowAroundSphereCodeGen_UBB ubb(blocks, pdfFieldIDGPU, velocity_initialisation);
lbm::FlowAroundSphereCodeGen_NoSlip noSlip(blocks, pdfFieldIDGPU);
lbm::FlowAroundSphereCodeGen_Outflow outflow(blocks, pdfFieldIDGPU, pdfFieldID);
#else
lbm::FlowAroundSphereCodeGen_UBB ubb(blocks, pdfFieldID, velocity_initialisation);
lbm::FlowAroundSphereCodeGen_NoSlip noSlip(blocks, pdfFieldID);
lbm::FlowAroundSphereCodeGen_Outflow outflow(blocks, pdfFieldID);
#endif
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldId, boundariesConfig);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldId, fluidFlagUID);
ubb.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("UBB"), fluidFlagUID);
noSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("NoSlip"), fluidFlagUID);
outflow.fillFromFlagField< FlagField_T >(blocks, flagFieldId, FlagUID("Outflow"), fluidFlagUID);
// create time loop
SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
#if defined(WALBERLA_BUILD_WITH_CUDA)
pystencils::FlowAroundSphereCodeGen_EvenSweep LBEvenSweep(densityFieldIDGPU, pdfFieldIDGPU, velFieldIDGPU, omega);
pystencils::FlowAroundSphereCodeGen_OddSweep LBOddSweep(densityFieldIDGPU, pdfFieldIDGPU, velFieldIDGPU, omega);
#else
pystencils::FlowAroundSphereCodeGen_EvenSweep LBEvenSweep(densityFieldID, pdfFieldID, velFieldID, omega);
pystencils::FlowAroundSphereCodeGen_OddSweep LBOddSweep(densityFieldID, pdfFieldID, velFieldID, omega);
#endif
// All the sweeps
auto tracker = make_shared< TimestepModulusTracker >(0);
AlternatingSweep LBSweep(LBEvenSweep, LBOddSweep, tracker);
AlternatingSweep noSlipSweep(noSlip.getEvenSweep(), noSlip.getOddSweep(), tracker);
AlternatingSweep outflowSweep(outflow.getEvenSweep(), outflow.getOddSweep(), tracker);
AlternatingSweep ubbSweep(ubb.getEvenSweep(), ubb.getOddSweep(), tracker);
AlternatingBeforeFunction communication(evenComm, oddComm, tracker);
// add LBM sweep and communication to time loop
timeloop.add() << BeforeFunction(communication, "communication")
<< Sweep(noSlipSweep, "noSlip boundary");
timeloop.add() << Sweep(outflowSweep, "outflow boundary");
timeloop.add() << Sweep(ubbSweep, "ubb boundary");
timeloop.add() << BeforeFunction(tracker->advancementFunction(), "Timestep Advancement")
<< Sweep(LBSweep, "LB update rule");
// LBM stability check
timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
config, blocks, pdfFieldID, flagFieldId, fluidFlagUID)),
"LBM stability check");
// log remaining time
timeloop.addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
"remaining time logger");
// add VTK output to time loop
uint_t vtkWriteFrequency = parameters.getParameter< uint_t >("vtkWriteFrequency", 0);
if (vtkWriteFrequency > 0)
{
auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out",
"simulation_step", false, true, true, false, 0);
#if defined(WALBERLA_BUILD_WITH_CUDA)
vtkOutput->addBeforeFunction([&]() {
cuda::fieldCpy< VelocityField_T, GPUField >(blocks, velFieldID, velFieldIDGPU);
cuda::fieldCpy< ScalarField_T, GPUField >(blocks, densityFieldID, densityFieldIDGPU);
});
#endif
auto velWriter = make_shared< field::VTKWriter< VelocityField_T > >(velFieldID, "velocity");
auto densityWriter = make_shared< field::VTKWriter< ScalarField_T > >(densityFieldID, "density");
FluidFilter_T filter(cellsPerBlock);
auto QCriterionWriter = make_shared<lbm::QCriterionVTKWriter<VelocityField_T, FluidFilter_T>>(blocks,
filter, velFieldID, "Q-Criterion");
vtkOutput->addCellDataWriter(velWriter);
vtkOutput->addCellDataWriter(densityWriter);
vtkOutput->addCellDataWriter(QCriterionWriter);
timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
}
WcTimer simTimer;
WALBERLA_LOG_INFO_ON_ROOT(
"Simulating flow around sphere:"
"\n timesteps: " << timesteps <<
"\n reynolds number: " << reynolds_number <<
"\n relaxation rate: " << omega <<
"\n maximum inflow velocity: " << u_max <<
"\n diameter_sphere: " << diameter_sphere)
simTimer.start();
timeloop.run();
simTimer.end();
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
auto time = simTimer.last();
auto nrOfCells = real_c(cellsPerBlock[0] * cellsPerBlock[1] * cellsPerBlock[2]);
auto mlupsPerProcess = nrOfCells * real_c(timesteps) / time * 1e-6;
WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process " << mlupsPerProcess)
WALBERLA_LOG_RESULT_ON_ROOT("Time per time step " << time / real_c(timesteps))
}
return EXIT_SUCCESS;
}
from lbmpy.advanced_streaming.utility import get_timesteps, Timestep
from pystencils.field import fields
from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
from lbmpy.stencils import get_stencil
from lbmpy.creationfunctions import create_lb_collision_rule, create_lb_method, create_lb_update_rule
from lbmpy.boundaries import NoSlip, UBB, ExtrapolationOutflow
from lbmpy_walberla.additional_data_handler import UBBAdditionalDataHandler, OutflowAdditionalDataHandler
from pystencils_walberla import CodeGeneration, generate_sweep
from lbmpy_walberla import RefinementScaling, generate_boundary, generate_lb_pack_info
import sympy as sp
stencil = get_stencil("D3Q27")
q = len(stencil)
dim = len(stencil[0])
streaming_pattern = 'esotwist'
timesteps = get_timesteps(streaming_pattern)
pdfs, velocity_field, density_field = fields(f"pdfs({q}), velocity({dim}), density(1) : double[{dim}D]", layout='fzyx')
omega = sp.Symbol("omega")
u_max = sp.Symbol("u_max")
output = {
'density': density_field,
'velocity': velocity_field
}
options = {'method': 'cumulant',
'stencil': stencil,
'relaxation_rate': omega,
'galilean_correction': True,
'field_name': 'pdfs',
'streaming_pattern': streaming_pattern,
'output': output,
'optimization': {'symbolic_field': pdfs,
'cse_global': False,
'cse_pdfs': False}}
method = create_lb_method(**options)
# getter & setter
setter_assignments = macroscopic_values_setter(method, velocity=velocity_field.center_vector,
pdfs=pdfs, density=1,
streaming_pattern=streaming_pattern,
previous_timestep=timesteps[0])
# opt = {'instruction_set': 'sse', 'assume_aligned': True, 'nontemporal': False, 'assume_inner_stride_one': True}
collision_rule = create_lb_collision_rule(lb_method=method, **options)
update_rule_even = create_lb_update_rule(collision_rule=collision_rule, timestep=Timestep.EVEN, **options)
update_rule_odd = create_lb_update_rule(collision_rule=collision_rule, timestep=Timestep.ODD, **options)
info_header = f"""
using namespace walberla;
#include "stencil/D{dim}Q{q}.h"
using Stencil_T = walberla::stencil::D{dim}Q{q};
using PdfField_T = GhostLayerField<real_t, {q}>;
using VelocityField_T = GhostLayerField<real_t, {dim}>;
using ScalarField_T = GhostLayerField<real_t, 1>;
"""
stencil = method.stencil
with CodeGeneration() as ctx:
if ctx.cuda:
target = 'gpu'
else:
target = 'cpu'
# sweeps
generate_sweep(ctx, 'FlowAroundSphereCodeGen_EvenSweep', update_rule_even, target=target)
generate_sweep(ctx, 'FlowAroundSphereCodeGen_OddSweep', update_rule_odd, target=target)
generate_sweep(ctx, 'FlowAroundSphereCodeGen_MacroSetter', setter_assignments, target=target)
# boundaries
ubb = UBB(lambda *args: None, dim=dim)
ubb_data_handler = UBBAdditionalDataHandler(stencil, ubb)
outflow = ExtrapolationOutflow(stencil[4], method, streaming_pattern=streaming_pattern)
outflow_data_handler = OutflowAdditionalDataHandler(stencil, outflow, target=target)
generate_boundary(ctx, 'FlowAroundSphereCodeGen_UBB', ubb, method,
target=target, streaming_pattern=streaming_pattern, always_generate_separate_classes=True,
additional_data_handler=ubb_data_handler)
generate_boundary(ctx, 'FlowAroundSphereCodeGen_NoSlip', NoSlip(), method, target=target,
streaming_pattern=streaming_pattern, always_generate_separate_classes=True)
generate_boundary(ctx, 'FlowAroundSphereCodeGen_Outflow', outflow, method, target=target,
streaming_pattern=streaming_pattern, always_generate_separate_classes=True,
additional_data_handler=outflow_data_handler)
# communication
generate_lb_pack_info(ctx, 'FlowAroundSphereCodeGen_PackInfo', stencil, pdfs,
streaming_pattern=streaming_pattern, always_generate_separate_classes=True, target=target)
# Info header containing correct template definitions for stencil and field
ctx.write_file("FlowAroundSphereCodeGen_InfoHeader.h", info_header)
import waLBerla as wlb
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
class Scenario:
def __init__(self):
self.timesteps = 1501
self.vtkWriteFrequency = 1500
self.cells = (384, 128, 128)
self.blocks = (1, 1, 1)
self.periodic = (0, 0, 0)
self.diameter_sphere = min(self.cells) // 2
self.u_max = 0.05
self.reynolds_number = 1000000
kinematic_viscosity = (self.diameter_sphere * self.u_max) / self.reynolds_number
self.omega = relaxation_rate_from_lattice_viscosity(kinematic_viscosity)
self.total_cells = (self.cells[0] * self.blocks[0],
self.cells[1] * self.blocks[1],
self.cells[2] * self.blocks[2])
@wlb.member_callback
def config(self):
return {
'DomainSetup': {
'blocks': self.blocks,
'cellsPerBlock': self.cells,
'periodic': self.periodic,
'oneBlockPerProcess': True
},
'Parameters': {
'timesteps': self.timesteps,
'vtkWriteFrequency': self.vtkWriteFrequency,
'omega': self.omega,
'u_max': self.u_max,
'reynolds_number': self.reynolds_number,
'diameter_sphere': self.diameter_sphere
},
'Boundaries': {
'Border': [
{'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'},
{'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'},
{'direction': 'W', 'walldistance': -1, 'flag': 'UBB'},
{'direction': 'E', 'walldistance': -1, 'flag': 'Outflow'},
{'direction': 'T', 'walldistance': -1, 'flag': 'NoSlip'},
{'direction': 'B', 'walldistance': -1, 'flag': 'NoSlip'},
],
'Body': [
{'shape': "sphere",
'midpoint': (int(0.40 * self.total_cells[0]), self.total_cells[1] // 2, self.total_cells[2] // 2),
'radius': self.diameter_sphere // 2,
'flag': 'NoSlip'}
]
},
}
scenarios = wlb.ScenarioManager()
scenarios.add(Scenario())
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment