Skip to content
Snippets Groups Projects
Commit bf963f4d authored by Markus Holzer's avatar Markus Holzer
Browse files

[skip ci] update

parent 76ce5c85
Branches
No related tags found
No related merge requests found
Pipeline #72104 skipped
...@@ -18,11 +18,11 @@ waLBerla_generate_target_from_python(NAME PorousMediaGPUCodeGen ...@@ -18,11 +18,11 @@ waLBerla_generate_target_from_python(NAME PorousMediaGPUCodeGen
if (WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP) if (WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
waLBerla_add_executable ( NAME PorousMediaGPU waLBerla_add_executable ( NAME PorousMediaGPU
FILES PorousMediaGPU.cpp Types.h InitializerFunctions.cpp FILES PorousMediaGPU.cpp GridGeneration.h Setup.h Types.h InitializerFunctions.cpp
DEPENDS blockforest boundary core field gpu lbm_generated stencil timeloop vtk PorousMediaGPUCodeGen ) DEPENDS blockforest boundary core field gpu lbm_generated postprocessing stencil timeloop vtk PorousMediaGPUCodeGen )
else() else()
waLBerla_add_executable ( NAME PorousMediaGPU waLBerla_add_executable ( NAME PorousMediaGPU
FILES PorousMediaGPU.cpp Types.h InitializerFunctions.cpp FILES PorousMediaGPU.cpp GridGeneration.h Setup.h Types.h InitializerFunctions.cpp
DEPENDS blockforest boundary core field lbm_generated stencil timeloop vtk PorousMediaGPUCodeGen ) DEPENDS blockforest boundary core field lbm_generated postprocessing stencil timeloop vtk PorousMediaGPUCodeGen )
endif(WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP) 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 GridGeneration.h
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
#pragma once
#include "blockforest/Initialization.h"
#include "blockforest/SetupBlock.h"
#include "blockforest/SetupBlockForest.h"
#include "blockforest/loadbalancing/StaticCurve.h"
#include "core/Environment.h"
#include "core/logging/Initialization.h"
#include "core/timing/RemainingTimeLogger.h"
#include "core/timing/TimingPool.h"
#include <string>
#include <fstream>
#include "Types.h"
#include "Setup.h"
using namespace walberla;
uint_t blockCells(const Setup& setup, const bool withGhostLayers){
uint_t cells;
if(!withGhostLayers){
cells = setup.cellsPerBlock[0] * setup.cellsPerBlock[1] * setup.cellsPerBlock[2];
}
else{
cells = (setup.cellsPerBlock[0] + 2 * setup.numGhostLayers) *
(setup.cellsPerBlock[1] + 2 * setup.numGhostLayers) *
(setup.cellsPerBlock[2] + 2 * setup.numGhostLayers);
}
return cells;
}
void createSetupBlockForest(SetupBlockForest& setupBfs, const Setup& setup, const bool useMPIManager=false)
{
WALBERLA_LOG_INFO_ON_ROOT("Generating SetupBlockForest...")
uint_t numProcesses = setup.numProcesses;
const std::string blockForestFilestem = setup.blockForestFilestem;
if(useMPIManager) {numProcesses = uint_c(mpi::MPIManager::instance()->numProcesses());}
// Sphere Sphere( setup );
// SphereRefinementSelection SphereRefinementSelection( Sphere, setup.refinementLevels );
// SphereBlockExclusion SphereBlockExclusion( Sphere );
// setupBfs.addRefinementSelectionFunction(std::function<void(SetupBlockForest &)>(SphereRefinementSelection));
// setupBfs.addBlockExclusionFunction(SphereBlockExclusion);
const AABB domain(real_t(0.0), real_t(0.0), real_t(0.0), setup.domainSize[0], setup.domainSize[1], setup.domainSize[2]);
setupBfs.addWorkloadMemorySUIDAssignmentFunction(blockforest::uniformWorkloadAndMemoryAssignment);
setupBfs.init(domain, setup.rootBlocks[0], setup.rootBlocks[1], setup.rootBlocks[2],
setup.periodic[0], setup.periodic[1], setup.periodic[2]);
setupBfs.balanceLoad(blockforest::StaticLevelwiseCurveBalanceWeighted(), numProcesses);
WALBERLA_LOG_INFO_ON_ROOT("=========================== BLOCK FOREST STATISTICS ============================");
WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << setupBfs.getNumberOfBlocks())
for (uint_t level = 0; level <= setup.refinementLevels; level++){
const uint_t numberOfBlocks = setupBfs.getNumberOfBlocks(level);
WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << numberOfBlocks)
}
const real_t avgBlocksPerProc = real_c(setupBfs.getNumberOfBlocks()) / real_c(setupBfs.getNumberOfProcesses());
WALBERLA_LOG_INFO_ON_ROOT("Average blocks per process: " << avgBlocksPerProc);
const uint_t cellsPerBlock = blockCells(setup, false);
const uint_t totalNumberCells = setupBfs.getNumberOfBlocks() * cellsPerBlock;
const real_t averageCellsPerGPU = avgBlocksPerProc * real_c(cellsPerBlock);
const uint_t cellsPerBlockGL = blockCells(setup, true);
const uint_t totalNumberCellsGL = setupBfs.getNumberOfBlocks() * cellsPerBlockGL;
const real_t averageCellsPerGPUGL = avgBlocksPerProc * real_c(cellsPerBlockGL);
const uint_t PDFsPerCell = StorageSpecification_T::inplace ? Stencil_T::Q : 2 * Stencil_T::Q;
const uint_t valuesPerCell = (PDFsPerCell + VelocityField_T::F_SIZE + ScalarField_T::F_SIZE);
const uint_t sizePerValue = sizeof(StorageSpecification_T::value_type);
const double expectedMemory = double_c(totalNumberCells * valuesPerCell * sizePerValue) * 1e-9;
const double expectedMemoryPerGPU = double_c(averageCellsPerGPU * valuesPerCell * sizePerValue) * 1e-9;
const double expectedMemoryGL = double_c(totalNumberCellsGL * valuesPerCell * sizePerValue) * 1e-9;
const double expectedMemoryPerGPUGL = double_c(averageCellsPerGPUGL * valuesPerCell * sizePerValue) * 1e-9;
WALBERLA_LOG_INFO_ON_ROOT( "Total number of cells will be " << totalNumberCells << " fluid cells (in total on all levels)")
WALBERLA_LOG_INFO_ON_ROOT( "Expected total memory demand will be " << expectedMemory << " GB")
WALBERLA_LOG_INFO_ON_ROOT( "Average memory demand per GPU will be " << expectedMemoryPerGPU << " GB")
WALBERLA_LOG_INFO_ON_ROOT( "Expected total memory demand (with GL) will be " << expectedMemoryGL << " GB")
WALBERLA_LOG_INFO_ON_ROOT( "Average memory demand per GPU (with GL) will be " << expectedMemoryPerGPUGL << " GB")
WALBERLA_LOG_INFO_ON_ROOT("=================================================================================")
if(mpi::MPIManager::instance()->numProcesses() > 1)
return;
if(setup.writeSetupForestAndReturn){
std::ostringstream oss;
oss << blockForestFilestem << ".bfs";
setupBfs.saveToFile(oss.str().c_str());
setupBfs.writeVTKOutput(blockForestFilestem);
}
}
void createBlockForest(shared_ptr< BlockForest >& bfs, const Setup& setup){
if (mpi::MPIManager::instance()->numProcesses() > 1){
// Load structured block forest from file
std::ostringstream oss;
oss << setup.blockForestFilestem << ".bfs";
const std::string setupBlockForestFilepath = oss.str();
std::ifstream infile(setupBlockForestFilepath.c_str());
if(!infile.good()){
WALBERLA_LOG_WARNING_ON_ROOT("Blockforest was not created beforehand and thus needs to be created on the fly. For large simulation runs this can be a severe problem!")
SetupBlockForest setupBfs;
createSetupBlockForest(setupBfs, setup, true);
bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs);
}
else{
bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()),
setupBlockForestFilepath.c_str(), false);
}
}
else{
SetupBlockForest setupBfs;
createSetupBlockForest(setupBfs, setup);
bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs);
}
}
\ No newline at end of file
...@@ -46,6 +46,7 @@ ...@@ -46,6 +46,7 @@
#include "field/vtk/FlagFieldCellFilter.h" #include "field/vtk/FlagFieldCellFilter.h"
#include "geometry/InitBoundaryHandling.h" #include "geometry/InitBoundaryHandling.h"
#include "geometry/mesh/TriangleMeshIO.h"
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT) #if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
# include "gpu/AddGPUFieldToStorage.h" # include "gpu/AddGPUFieldToStorage.h"
...@@ -68,6 +69,8 @@ ...@@ -68,6 +69,8 @@
# include "lbm_generated/gpu/UniformGeneratedGPUPdfPackInfo.h" # include "lbm_generated/gpu/UniformGeneratedGPUPdfPackInfo.h"
#endif #endif
#include "postprocessing/FieldToSurfaceMesh.h"
#include "timeloop/SweepTimeloop.h" #include "timeloop/SweepTimeloop.h"
#include "vtk/VTKOutput.h" #include "vtk/VTKOutput.h"
...@@ -76,6 +79,8 @@ ...@@ -76,6 +79,8 @@
#include <memory> #include <memory>
#include "InitializerFunctions.h" #include "InitializerFunctions.h"
#include "GridGeneration.h"
#include "Setup.h"
#include "Types.h" #include "Types.h"
#include "PorousMediaGPUStaticDefines.h" #include "PorousMediaGPUStaticDefines.h"
...@@ -98,8 +103,7 @@ using blockforest::communication::UniformBufferedScheme; ...@@ -98,8 +103,7 @@ using blockforest::communication::UniformBufferedScheme;
using lbm_generated::UniformGeneratedPdfPackInfo; using lbm_generated::UniformGeneratedPdfPackInfo;
#endif #endif
int main(int argc, char** argv) int main(int argc, char** argv){
{
Environment env( argc, argv ); Environment env( argc, argv );
mpi::MPIManager::instance()->useWorldComm(); mpi::MPIManager::instance()->useWorldComm();
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT) #if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
...@@ -112,54 +116,48 @@ int main(int argc, char** argv) ...@@ -112,54 +116,48 @@ int main(int argc, char** argv)
/// SETUP AND CONFIGURATION /// /// SETUP AND CONFIGURATION ///
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
auto blocks = blockforest::createUniformBlockGridFromConfig(config); auto parameters = config->getOneBlock("Parameters");
auto domainParameters = config->getOneBlock("DomainSetup");
auto domainParameters = config->getOneBlock("DomainSetup"); auto loggingParameters = config->getOneBlock("Logging");
const real_t coarseMeshSize = domainParameters.getParameter< real_t >("dx"); auto boundariesConfig = config->getBlock("Boundaries");
auto parameters = config->getOneBlock("Parameters");
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(50));
const bool gpuEnabledMPI = parameters.getParameter< bool >("gpuEnabledMPI", false);
const real_t reynoldsNumber = parameters.getParameter< real_t >("reynoldsNumber");
const real_t referenceVelocity = parameters.getParameter< real_t >("referenceVelocity");
const real_t referenceLength = parameters.getParameter< real_t >("referenceLength");
const real_t latticeVelocity = parameters.getParameter< real_t >("latticeVelocity");
const real_t waveNumber = parameters.getParameter< real_t >("waveNumber");
const real_t scaling = parameters.getParameter< real_t >("scaling");
const real_t speedOfSound = real_c(real_c(1.0) / std::sqrt( real_c(3.0) ));
const real_t machNumber = latticeVelocity / speedOfSound;
// TODO define RE
const real_t viscosity = real_c((latticeVelocity * referenceLength) / reynoldsNumber);
const real_t omega = real_c(real_c(1.0) / (real_c(3.0) * viscosity + real_c(0.5)));
const real_t dt = latticeVelocity / referenceVelocity * coarseMeshSize;
const StorageSpecification_T StorageSpec = StorageSpecification_T();
WALBERLA_LOG_INFO_ON_ROOT( "Benchmark run data:" Setup setup(parameters, domainParameters, loggingParameters, infoMap);
WALBERLA_LOG_INFO_ON_ROOT( "Simulation run data:"
"\n- simulation parameters:" << std::setprecision(16) << "\n- simulation parameters:" << std::setprecision(16) <<
"\n + collision model: " << infoMap["collisionOperator"] << "\n + collision model: " << setup.collisionModel <<
"\n + stencil: " << infoMap["stencil"] << "\n + stencil: " << setup.stencil <<
"\n + streaming: " << infoMap["streamingPattern"] << "\n + streaming: " << setup.streamingPattern <<
"\n + compressible: " << ( StorageSpecification_T::compressible ? "yes" : "no" ) << "\n + compressible: " << ( StorageSpecification_T::compressible ? "yes" : "no" ) <<
"\n + resolution: " << coarseMeshSize << " - on the coarsest grid" << "\n + resolution: " << setup.dxC << " - on the coarsest grid" <<
"\n- simulation properties: " "\n- simulation properties: "
"\n + kin. viscosity: " << viscosity * coarseMeshSize * coarseMeshSize / dt << " [m^2/s] (on the coarsest grid)" << "\n + kin. viscosity: " << setup.kinViscosity << " [m^2/s] (on the coarsest grid)" <<
"\n + viscosity LB: " << viscosity << " [dx*dx/dt] (on the coarsest grid)" << "\n + viscosity LB: " << setup.viscosity << " [dx*dx/dt] (on the coarsest grid)" <<
"\n + omega: " << omega << " (on the coarsest grid)" << "\n + omega: " << setup.omega << " (on the coarsest grid)" <<
"\n + inflow velocity: " << referenceVelocity << " [m/s]" << "\n + inflow velocity: " << setup.referenceVelocity << " [m/s]" <<
"\n + lattice velocity: " << latticeVelocity << " [dx/dt]" << "\n + lattice velocity: " << setup.latticeVelocity << " [dx/dt]" <<
"\n + Reynolds number: " << reynoldsNumber << "\n + Reynolds number: " << setup.reynoldsNumber <<
"\n + Mach number: " << machNumber << "\n + Mach number: " << setup.machNumber <<
"\n + dx (coarsest grid): " << coarseMeshSize << " [m]" << "\n + dx (coarsest grid): " << setup.dxC << " [m]" <<
"\n + dt (coarsest grid): " << dt << " [s]" "\n + dt (coarsest grid): " << setup.dt << " [s]" <<
"\n + #time steps: " << timesteps << " (on the coarsest grid, " << ( real_c(1.0) / dt ) << " for 1s of real time)" << "\n + conversion Faktor: " << setup.conversionFactor <<
"\n + simulation time: " << ( real_c( timesteps ) * 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]" )
bool writeSetupForestAndReturn = loggingParameters.getParameter< bool >("writeSetupForestAndReturn", false);
if (uint_c(MPIManager::instance()->numProcesses()) > 1) writeSetupForestAndReturn = false;
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")
logging::Logging::printFooterOnStream();
return EXIT_SUCCESS;
}
auto blocks = std::make_shared< StructuredBlockForest >(bfs, setup.cellsPerBlock[0], setup.cellsPerBlock[1], setup.cellsPerBlock[2]);
blocks->createCellBoundingBoxes();
const StorageSpecification_T StorageSpec = StorageSpecification_T();
// Creating fields // Creating fields
IDs ids; IDs ids;
ids.pdfField = lbm_generated::addPdfFieldToStorage(blocks, "pdfs", StorageSpec, uint_c(1), field::fzyx); ids.pdfField = lbm_generated::addPdfFieldToStorage(blocks, "pdfs", StorageSpec, uint_c(1), field::fzyx);
...@@ -174,10 +172,8 @@ int main(int argc, char** argv) ...@@ -174,10 +172,8 @@ int main(int argc, char** argv)
ids.densityFieldGPU = gpu::addGPUFieldToStorage< ScalarField_T >(blocks, ids.densityField, "velocity on GPU", true); ids.densityFieldGPU = gpu::addGPUFieldToStorage< ScalarField_T >(blocks, ids.densityField, "velocity on GPU", true);
const Cell innerOuterSplit = Cell(parameters.getParameter< Vector3<cell_idx_t> >("innerOuterSplit", Vector3<cell_idx_t>(1, 1, 1))); const Cell innerOuterSplit = Cell(parameters.getParameter< Vector3<cell_idx_t> >("innerOuterSplit", Vector3<cell_idx_t>(1, 1, 1)));
SweepCollection_T sweepCollection(blocks, ids.pdfFieldGPU, ids.densityFieldGPU, ids.velocityFieldGPU, omega, innerOuterSplit); SweepCollection_T sweepCollection(blocks, ids.pdfFieldGPU, ids.densityFieldGPU, ids.velocityFieldGPU, setup.omega, innerOuterSplit);
for (auto& block : *blocks){
for (auto& block : *blocks)
{
sweepCollection.initialise(&block); sweepCollection.initialise(&block);
} }
...@@ -186,25 +182,28 @@ int main(int argc, char** argv) ...@@ -186,25 +182,28 @@ int main(int argc, char** argv)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Boundaries // Boundaries
const FlagUID fluidFlagUID("Fluid"); if (boundariesConfig){
const FlagUID wallFlagUID("NoSlip");
auto boundariesConfig = config->getBlock("Boundaries");
if (boundariesConfig)
{
WALBERLA_LOG_INFO_ON_ROOT("Setting boundary conditions") WALBERLA_LOG_INFO_ON_ROOT("Setting boundary conditions")
geometry::initBoundaryHandling< FlagField_T >(*blocks, ids.flagField, boundariesConfig); geometry::initBoundaryHandling< FlagField_T >(*blocks, ids.flagField, boundariesConfig);
init_TPMS(blocks, ids.flagField, wallFlagUID, waveNumber, scaling); init_TPMS(blocks, ids.flagField, setup.obstacleUID, setup.waveNumber, setup.scaling);
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, fluidFlagUID); geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, setup.fluidUID);
} }
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, fluidFlagUID); geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, setup.fluidUID);
// BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, fluidFlagUID, latticeVelocity); // BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, fluidFlagUID, latticeVelocity);
BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, fluidFlagUID, latticeVelocity, ids.pdfField); BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, setup.fluidUID, setup.latticeVelocity, ids.pdfField);
if( loggingParameters.getParameter< bool >("writeSurfaceMeshOfTPMS", false) ){
WALBERLA_LOG_INFO_ON_ROOT("Writing TPMS as surface file")
auto mesh = postprocessing::flagFieldToSurfaceMesh< FlagField_T >(blocks, ids.flagField, 4, true);
WALBERLA_EXCLUSIVE_WORLD_SECTION(0){
geometry::writeMesh("TPMS.obj", *mesh);
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// COMMUNICATION SCHEME /// /// COMMUNICATION SCHEME ///
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
const bool gpuEnabledMPI = parameters.getParameter< bool >("gpuEnabledMPI", false);
UniformGPUScheme< Stencil_T > communication(blocks, gpuEnabledMPI); UniformGPUScheme< Stencil_T > communication(blocks, gpuEnabledMPI);
auto packInfo = std::make_shared<lbm_generated::UniformGeneratedGPUPdfPackInfo< GPUPdfField_T >>(ids.pdfFieldGPU); auto packInfo = std::make_shared<lbm_generated::UniformGeneratedGPUPdfPackInfo< GPUPdfField_T >>(ids.pdfFieldGPU);
communication.addPackInfo(packInfo); communication.addPackInfo(packInfo);
...@@ -215,7 +214,7 @@ int main(int argc, char** argv) ...@@ -215,7 +214,7 @@ int main(int argc, char** argv)
const int streamLowPriority = 0; const int streamLowPriority = 0;
auto defaultStream = gpu::StreamRAII::newPriorityStream(streamLowPriority); auto defaultStream = gpu::StreamRAII::newPriorityStream(streamLowPriority);
SweepTimeloop timeLoop(blocks->getBlockStorage(), timesteps); SweepTimeloop timeLoop(blocks->getBlockStorage(), setup.timesteps);
timeLoop.add() << BeforeFunction(communication.getCommunicateFunctor(), "communication") timeLoop.add() << BeforeFunction(communication.getCommunicateFunctor(), "communication")
<< Sweep(boundaryCollection.getSweep(BoundaryCollection_T::ALL, defaultStream), "Boundary Conditions"); << Sweep(boundaryCollection.getSweep(BoundaryCollection_T::ALL, defaultStream), "Boundary Conditions");
timeLoop.add() << Sweep(sweepCollection.streamCollide(SweepCollection_T::ALL, defaultStream), "LBM StreamCollide"); timeLoop.add() << Sweep(sweepCollection.streamCollide(SweepCollection_T::ALL, defaultStream), "LBM StreamCollide");
...@@ -225,38 +224,81 @@ int main(int argc, char** argv) ...@@ -225,38 +224,81 @@ int main(int argc, char** argv)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// VTK // VTK
auto vtkParameters = config->getOneBlock("VTKWriter"); auto VTKWriter = config->getOneBlock("VTKWriter");
const uint_t vtkWriteFrequency = vtkParameters.getParameter< uint_t >("vtkWriteFrequency", 0); const uint_t vtkWriteFrequency = VTKWriter.getParameter< uint_t >("vtkWriteFrequency", 0);
if (vtkWriteFrequency > 0) const uint_t vtkGhostLayers = VTKWriter.getParameter< uint_t >("vtkGhostLayers", 0);
{ const bool amrFileFormat = VTKWriter.getParameter< bool >("amrFileFormat", false);
auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, 0, false, "vtk_out", const bool oneFilePerProcess = VTKWriter.getParameter< bool >("oneFilePerProcess", false);
"simulation_step", false, true, true, false, 0); const real_t samplingResolution = VTKWriter.getParameter< real_t >("samplingResolution", real_c(-1.0));
auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(ids.velocityField, "vel"); const uint_t initialWriteCallsToSkip = VTKWriter.getParameter< uint_t >("initialWriteCallsToSkip", uint_c(0.0));
vtkOutput->addCellDataWriter(velWriter); const real_t sliceThickness = VTKWriter.getParameter< real_t >("sliceThickness", real_c(1.0));
auto finalDomain = blocks->getDomain();
if (vtkWriteFrequency > 0){
const std::string vtkFolder = "VTKPorousMediaRE_" + 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([&]() { vtkOutput->addBeforeFunction([&]() {
for (auto& block : *blocks) for (auto& block : *blocks){
sweepCollection.calculateMacroscopicParameters(&block); sweepCollection.calculateMacroscopicParameters(&block);
}
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
gpu::fieldCpy< VelocityField_T, gpu::GPUField< VelocityField_T::value_type > >(blocks, ids.velocityField, ids.velocityFieldGPU); gpu::fieldCpy< VelocityField_T, gpu::GPUField< VelocityField_T::value_type > >(blocks, ids.velocityField, ids.velocityFieldGPU);
gpu::fieldCpy< ScalarField_T, gpu::GPUField< ScalarField_T::value_type > >(blocks, ids.densityField, ids.densityFieldGPU);
#endif
}); });
if (vtkParameters.getParameter< bool >("flag")) 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")){
if (VTKWriter.getParameter< bool >("writeVelocityAsMagnitude")){
auto velWriter = make_shared< field::VTKMagnitudeWriter< VelocityField_T, float32 > >(ids.velocityField, "velocityMagnitude");
vtkOutput->addCellDataWriter(velWriter);
}
else{
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 >("flag")){
auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(ids.flagField, "flag"); auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(ids.flagField, "flag");
vtkOutput->addCellDataWriter(flagWriter); vtkOutput->addCellDataWriter(flagWriter);
} }
timeLoop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output"); timeLoop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
} }
// log remaining time // log remaining time
auto loggingParameters= config->getOneBlock("Logging");
const real_t remainingTimeLoggerFrequency = if (uint_c(setup.remainingTimeLoggerFrequency) > 0){
loggingParameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0); // in seconds
if (uint_c(remainingTimeLoggerFrequency) > 0)
{
timeLoop.addFuncAfterTimeStep( timeLoop.addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeLoop.getNrOfTimeSteps(), remainingTimeLoggerFrequency), timing::RemainingTimeLogger(timeLoop.getNrOfTimeSteps(), setup.remainingTimeLoggerFrequency),
"remaining time logger"); "remaining time logger");
} }
...@@ -266,8 +308,8 @@ int main(int argc, char** argv) ...@@ -266,8 +308,8 @@ int main(int argc, char** argv)
WALBERLA_GPU_CHECK(gpuPeekAtLastError()) WALBERLA_GPU_CHECK(gpuPeekAtLastError())
#endif #endif
const lbm_generated::PerformanceEvaluation< FlagField_T > performance(blocks, ids.flagField, fluidFlagUID); const lbm_generated::PerformanceEvaluation< FlagField_T > performance(blocks, ids.flagField, setup.fluidUID);
field::CellCounter< FlagField_T > fluidCells(blocks, ids.flagField, fluidFlagUID); field::CellCounter< FlagField_T > fluidCells(blocks, ids.flagField, setup.fluidUID);
fluidCells(); fluidCells();
WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation") WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation")
...@@ -291,7 +333,7 @@ int main(int argc, char** argv) ...@@ -291,7 +333,7 @@ int main(int argc, char** argv)
WALBERLA_LOG_INFO_ON_ROOT("Simulation finished") WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
double time = double_c(simTimer.max()); double time = double_c(simTimer.max());
WALBERLA_MPI_SECTION() { walberla::mpi::reduceInplace(time, walberla::mpi::MAX); } WALBERLA_MPI_SECTION() { walberla::mpi::reduceInplace(time, walberla::mpi::MAX); }
performance.logResultOnRoot(timesteps, time); performance.logResultOnRoot(setup.timesteps, time);
const auto reducedTimeloopTiming = timeloopTiming.getReduced(); const auto reducedTimeloopTiming = timeloopTiming.getReduced();
WALBERLA_LOG_RESULT_ON_ROOT("Time loop timing:\n" << *reducedTimeloopTiming) WALBERLA_LOG_RESULT_ON_ROOT("Time loop timing:\n" << *reducedTimeloopTiming)
......
...@@ -10,7 +10,7 @@ Parameters ...@@ -10,7 +10,7 @@ Parameters
waveNumber 6.2831853072; waveNumber 6.2831853072;
scaling 10.0; scaling 10.0;
timesteps 50001; timesteps 10;
processMemoryLimit 512; // MiB processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>; innerOuterSplit <1, 1, 1>;
...@@ -22,10 +22,13 @@ Parameters ...@@ -22,10 +22,13 @@ Parameters
//! [domainSetup] //! [domainSetup]
DomainSetup DomainSetup
{ {
cellsPerBlock < 1536, 192, 192 >; domainSize < 8, 1, 1 >;
cellsPerBlock < 384, 48, 48 >;
blocks < 1, 1, 1 >; blocks < 1, 1, 1 >;
periodic < false, true, true >; periodic < false, true, true >;
dx 0.00520833333333; refinementLevels 0;
numberProcesses 1; // This is for load balancing, overwritten if more than one proc is used
blockForestFilestem porousMediaBlockForest;
} }
//! [domainSetup] //! [domainSetup]
...@@ -45,16 +48,16 @@ StabilityChecker ...@@ -45,16 +48,16 @@ StabilityChecker
VTKWriter VTKWriter
{ {
vtkWriteFrequency 5000; vtkWriteFrequency 5;
vtkGhostLayers 0; vtkGhostLayers 0;
velocity true; velocity true;
density true; writeVelocityAsMagnitude true;
flag true; density false;
omega false; flag false;
writeOnlySlice true; writeOnlySlice true;
sliceThickness 1; sliceThickness 1;
writeXZSlice false; writeXZSlice false;
amrFileFormat true; amrFileFormat false;
oneFilePerProcess false; oneFilePerProcess false;
samplingResolution -1; samplingResolution -1;
initialWriteCallsToSkip 0; initialWriteCallsToSkip 0;
...@@ -63,9 +66,10 @@ VTKWriter ...@@ -63,9 +66,10 @@ VTKWriter
Logging Logging
{ {
logLevel info; // info progress detail tracing logLevel info; // info progress detail tracing
writeSetupForestAndReturn true; // When only one process is used the decomposition is writen and the program terminates writeSetupForestAndReturn false; // When only one process is used the decomposition is writen and the program terminates
readSetupFromFile false; readSetupFromFile false;
remainingTimeLoggerFrequency 20; // in seconds remainingTimeLoggerFrequency 20; // in seconds
writeSurfaceMeshOfTPMS false;
} }
Evaluation Evaluation
......
//======================================================================================================================
//
// 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/config/Config.h"
#include "core/DataTypes.h"
#include "core/math/Vector3.h"
#include <string>
#include "Types.h"
namespace walberla
{
struct Setup
{
Setup(const Config::BlockHandle & parameters, const Config::BlockHandle & domainParameters,
const Config::BlockHandle & logging, std::map<std::string, std::string>& infoMap)
{
blockForestFilestem = domainParameters.getParameter< std::string >("blockForestFilestem", "blockforest");
numProcesses = domainParameters.getParameter< uint_t >("numberProcesses");
domainSize = domainParameters.getParameter< Vector3< real_t > >("domainSize");
cellsPerBlock = domainParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
rootBlocks = domainParameters.getParameter< Vector3< uint_t > >("blocks");
cells = Vector3<uint_t>(rootBlocks[0] * cellsPerBlock[0], rootBlocks[1] * cellsPerBlock[1], rootBlocks[2] * cellsPerBlock[2]);
periodic = domainParameters.getParameter< Vector3< bool > >("periodic");
refinementLevels = domainParameters.getParameter< uint_t >("refinementLevels");
const real_t coarseMeshSize = domainSize[0] / real_c(rootBlocks[0] * cellsPerBlock[0]);
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(coarseMeshSize, domainSize[1] / real_c(rootBlocks[1] * cellsPerBlock[1]), real_c(1e-12), "Resulting resolution is not equal in x and y direxction. This is not supported");
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(coarseMeshSize, domainSize[2] / real_c(rootBlocks[2] * cellsPerBlock[2]), real_c(1e-12), "Resulting resolution is not equal in x and z direxction. This is not supported");
referenceLength = domainSize[1];
reynoldsNumber = parameters.getParameter< real_t >("reynoldsNumber");
referenceVelocity = parameters.getParameter< real_t >("referenceVelocity");
latticeVelocity = parameters.getParameter< real_t >("latticeVelocity");
waveNumber = parameters.getParameter< real_t >("waveNumber");
scaling = parameters.getParameter< real_t >("scaling");
const real_t speedOfSound = real_c(real_c(1.0) / std::sqrt( real_c(3.0) ));
machNumber = latticeVelocity / speedOfSound;
viscosity = real_c((latticeVelocity * referenceLength) / reynoldsNumber);
omega = real_c(real_c(1.0) / (real_c(3.0) * viscosity + real_c(0.5)));
rho = real_c(1.0);
dxC = coarseMeshSize;
dxF = real_c(coarseMeshSize) / real_c( 1 << refinementLevels );
dt = latticeVelocity / referenceVelocity * coarseMeshSize;
conversionFactor = dxC / dt;
kinViscosity = viscosity * coarseMeshSize * coarseMeshSize / dt;
timesteps = parameters.getParameter< uint_t >("timesteps");
numGhostLayers = refinementLevels == 0 ? uint_c(1) : uint_c(2);
valuesPerCell = (Stencil_T::Q + VelocityField_T::F_SIZE + uint_c(1) * ScalarField_T::F_SIZE);
memoryPerCell = memory_t(valuesPerCell * sizeof(PdfField_T::value_type) + uint_c(1));
processMemoryLimit = parameters.getParameter< memory_t >( "processMemoryLimit", memory_t( 512 ) ) * memory_t( 1024 * 1024 );
stencil = infoMap["stencil"];
streamingPattern = infoMap["streamingPattern"];
collisionModel = infoMap["collisionOperator"];
fluidUID = FlagUID("Fluid");
obstacleUID = FlagUID("NoSlip");
inflowUID = FlagUID("UBB");
outflowUID = FlagUID("Outflow");
writeSetupForestAndReturn = logging.getParameter< bool >("writeSetupForestAndReturn", false);
remainingTimeLoggerFrequency = logging.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
}
std::string blockForestFilestem;
uint_t numProcesses;
Vector3<uint_t> rootBlocks;
Vector3<uint_t> cells;
Vector3<real_t> domainSize;
Vector3< bool > periodic;
uint_t refinementLevels;
Vector3< uint_t > cellsPerBlock;
uint_t numGhostLayers;
real_t referenceVelocity;
real_t referenceLength;
real_t latticeVelocity;
real_t waveNumber;
real_t scaling;
real_t machNumber;
real_t reynoldsNumber;
real_t viscosity;
real_t kinViscosity;
real_t omega;
real_t rho;
real_t dxC;
real_t dxF;
real_t dt;
real_t conversionFactor;
uint_t timesteps;
uint_t valuesPerCell;
memory_t memoryPerCell;
memory_t processMemoryLimit;
std::string stencil;
std::string streamingPattern;
std::string collisionModel;
FlagUID fluidUID;
FlagUID obstacleUID;
FlagUID inflowUID;
FlagUID outflowUID;
bool writeSetupForestAndReturn;
real_t remainingTimeLoggerFrequency;
};
}
\ No newline at end of file
...@@ -105,7 +105,7 @@ class VTKWriter : public vtk::BlockCellDataWriter< OutputType, ...@@ -105,7 +105,7 @@ class VTKWriter : public vtk::BlockCellDataWriter< OutputType,
public: public:
using OutputTrait = VectorTrait<typename Field_T::value_type>; using OutputTrait = VectorTrait<typename Field_T::value_type>;
using base_t = vtk::BlockCellDataWriter<OutputType, OutputTrait::F_SIZE * Field_T::F_SIZE>; using base_t = vtk::BlockCellDataWriter<OutputType, OutputTrait::F_SIZE * Field_T::F_SIZE>;
VTKWriter( const ConstBlockDataID bdid, const std::string& id ) : VTKWriter( const ConstBlockDataID bdid, const std::string& id ) :
base_t( id ), bdid_( bdid ), field_( nullptr ) {} base_t( id ), bdid_( bdid ), field_( nullptr ) {}
...@@ -140,6 +140,46 @@ protected: ...@@ -140,6 +140,46 @@ protected:
}; };
template< typename Field_T, typename OutputType = float >
class VTKMagnitudeWriter : public vtk::BlockCellDataWriter< OutputType, 1 >
{
public:
using base_t = vtk::BlockCellDataWriter<OutputType, 1>;
using type_t = typename Field_T::value_type;
VTKMagnitudeWriter( const ConstBlockDataID bdid, const std::string& id ) :
base_t( id ), bdid_( bdid ), field_( nullptr ) {}
protected:
void configure() override {
WALBERLA_ASSERT_NOT_NULLPTR( this->block_ );
field_ = this->block_->template getData< Field_T >( bdid_ );
}
OutputType evaluate( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const cell_idx_t /*f*/ ) override
{
WALBERLA_ASSERT_NOT_NULLPTR( field_ );
if ( Field_T::F_SIZE == 1 ){
return numeric_cast<OutputType>( field_->get(x,y,z,0) );
}
else if (Field_T::F_SIZE == 2){
return numeric_cast<OutputType>( std::sqrt( field_->get(x,y,z,0)*field_->get(x,y,z,0) + field_->get(x,y,z,1)*field_->get(x,y,z,1) ));
}
else if (Field_T::F_SIZE == 3){
return numeric_cast<OutputType>( std::sqrt( field_->get(x,y,z,0)*field_->get(x,y,z,0) + field_->get(x,y,z,1)*field_->get(x,y,z,1) + field_->get(x,y,z,2)*field_->get(x,y,z,2) ));
}
else{
WALBERLA_ABORT("Only fields with 1, 2 or 3 entries per cell are supported");
}
}
const ConstBlockDataID bdid_;
const Field_T * field_;
}; // class VTKMagnitudeWriter
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment