diff --git a/apps/gpu/FlowAroundSphere/FlowAroundSphere.cpp b/apps/gpu/FlowAroundSphere/FlowAroundSphere.cpp index e1fea619653f017968cbd1284dc22cbadf73567c..bbc69fa186805c95f6c106ed0ddbb92e00d54731 100644 --- a/apps/gpu/FlowAroundSphere/FlowAroundSphere.cpp +++ b/apps/gpu/FlowAroundSphere/FlowAroundSphere.cpp @@ -68,26 +68,6 @@ ////////////////////////////////////////////////////////////////////////// -#include "GksMeshAdapter/GksMeshAdapter.h" - -#include "GksGpu/DataBase/DataBase.h" -#include "GksGpu/Initializer/Initializer.h" -#include "GksGpu/Parameters/Parameters.h" - -#include "GksGpu/FlowStateData/FlowStateDataConversion.cuh" - -#include "GksGpu/BoundaryConditions/BoundaryCondition.h" -#include "GksGpu/BoundaryConditions/IsothermalWall.h" - -#include "GksGpu/TimeStepping/NestedTimeStep.h" - -#include "GksGpu/Analyzer/ConvergenceAnalyzer.h" -#include "GksGpu/Analyzer/CupsAnalyzer.h" - -#include "GksGpu/CudaUtility/CudaUtility.h" - -#include "GksGpu/Output/VtkWriter.h" - int main(int argc, char *argv[]) { try { @@ -97,19 +77,15 @@ int main(int argc, char *argv[]) std::string outputPath("./output/Sphere"); std::string simulationName("Sphere"); - const real L = 1.0; - const real Re = 1000.0; + const real L = 1.0; + const real Re = 1000.0; const real velocity = 1.0; - const real dt = (real)0.5e-3; - const uint nx = 64; + const real dt = (real)0.5e-3; + const uint nx = 64; const uint timeStepOut = 10000; const uint timeStepEnd = 250000; - // switch between LBM and GKS solver here - // LbmOrGks lbmOrGks = GKS; - LbmOrGks lbmOrGks = LBM; - ////////////////////////////////////////////////////////////////////////// // setup logger ////////////////////////////////////////////////////////////////////////// @@ -134,217 +110,89 @@ int main(int argc, char *argv[]) gridBuilder->addCoarseGrid(-1.0 * L, -1.0 * L, -1.0 * L, 1.0 * L, 1.0 * L, 1.0 * L, dx); // use primitive - Object* sphere = new Sphere(0.0, 0.0, 0.0, 0.2); + Object *sphere = new Sphere(0.0, 0.0, 0.0, 0.2); // use stl // Object* sphere = TriangularMesh::make("stl/sphere.stl"); gridBuilder->addGeometry(sphere); gridBuilder->setPeriodicBoundaryCondition(false, false, false); - gridBuilder->buildGrids(lbmOrGks, false); + gridBuilder->buildGrids(LBM, false); ////////////////////////////////////////////////////////////////////////// - // branch between LBM and GKS + // setup simulation parameters ////////////////////////////////////////////////////////////////////////// - if (lbmOrGks == LBM) { - SPtr<Parameter> para = Parameter::make(); - - ////////////////////////////////////////////////////////////////////////// - // compute parameters in lattice units - ////////////////////////////////////////////////////////////////////////// - - const real velocityLB = velocity * dt / dx; // LB units - - const real vx = velocityLB / sqrt(2.0); // LB units - const real vy = velocityLB / sqrt(2.0); // LB units - - const real viscosityLB = nx * velocityLB / Re; // LB units - - *logging::out << logging::Logger::INFO_HIGH << "velocity [dx/dt] = " << velocityLB << " \n"; - *logging::out << logging::Logger::INFO_HIGH << "viscosity [dx^2/dt] = " << viscosityLB << "\n"; - - ////////////////////////////////////////////////////////////////////////// - // set parameters - ////////////////////////////////////////////////////////////////////////// - - para->setOutputPath(outputPath); - para->setOutputPrefix(simulationName); - - para->setPathAndFilename(para->getOutputPath() + "/" + para->getOutputPrefix()); - - para->setPrintFiles(true); - - para->setVelocityLB(velocityLB); - para->setViscosityLB(viscosityLB); - - para->setVelocityRatio(velocity / velocityLB); - - para->setTimestepOut(timeStepOut); - para->setTimestepEnd(timeStepEnd); - - ////////////////////////////////////////////////////////////////////////// - // set boundary conditions - ////////////////////////////////////////////////////////////////////////// - - gridBuilder->setNoSlipBoundaryCondition(SideType::PX); - gridBuilder->setNoSlipBoundaryCondition(SideType::MX); - gridBuilder->setNoSlipBoundaryCondition(SideType::PY); - gridBuilder->setNoSlipBoundaryCondition(SideType::MY); - gridBuilder->setVelocityBoundaryCondition(SideType::PZ, vx, vy, 0.0); - gridBuilder->setNoSlipBoundaryCondition(SideType::MZ); - gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, 0.0, 0.0); - - ////////////////////////////////////////////////////////////////////////// - // set copy mesh to simulation - ////////////////////////////////////////////////////////////////////////// - - SPtr<CudaMemoryManager> cudaMemoryManager = CudaMemoryManager::make(para); - - SPtr<GridProvider> gridGenerator = GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager); - - ////////////////////////////////////////////////////////////////////////// - // run simulation - ////////////////////////////////////////////////////////////////////////// - - Simulation sim; - SPtr<FileWriter> fileWriter = SPtr<FileWriter>(new FileWriter()); - sim.init(para, gridGenerator, fileWriter, cudaMemoryManager); - sim.run(); - sim.free(); - } else { - CudaUtility::setCudaDevice(0); - - Parameters parameters; - - ////////////////////////////////////////////////////////////////////////// - // compute remaining parameters - ////////////////////////////////////////////////////////////////////////// - - const real vx = velocity / sqrt(2.0); - const real vy = velocity / sqrt(2.0); - - parameters.K = 2.0; - parameters.Pr = 1.0; + SPtr<Parameter> para = Parameter::make(); - const real Ma = (real)0.1; - - real rho = 1.0; - - real cs = velocity / Ma; - real lambda = c1o2 * ((parameters.K + 5.0) / (parameters.K + 3.0)) / (cs * cs); - - const real mu = velocity * L * rho / Re; - - *logging::out << logging::Logger::INFO_HIGH << "mu = " << mu << " m^2/s\n"; - - *logging::out << logging::Logger::INFO_HIGH << "CFL = " << dt * (velocity + cs) / dx << "\n"; - - ////////////////////////////////////////////////////////////////////////// - // set parameters - ////////////////////////////////////////////////////////////////////////// - - parameters.mu = mu; - - parameters.dt = dt; - parameters.dx = dx; - - parameters.lambdaRef = lambda; - - ////////////////////////////////////////////////////////////////////////// - // set copy mesh to simulation - ////////////////////////////////////////////////////////////////////////// - - GksMeshAdapter meshAdapter(gridBuilder); - - meshAdapter.inputGrid(); - - auto dataBase = std::make_shared<DataBase>("GPU"); - - ////////////////////////////////////////////////////////////////////////// - // set boundary conditions - ////////////////////////////////////////////////////////////////////////// - - SPtr<BoundaryCondition> bcLid = - std::make_shared<IsothermalWall>(dataBase, Vec3(vx, vy, 0.0), lambda, false); - SPtr<BoundaryCondition> bcWall = - std::make_shared<IsothermalWall>(dataBase, Vec3(0.0, 0.0, 0.0), lambda, false); - - bcLid->findBoundaryCells(meshAdapter, false, [&](Vec3 center) { - return center.z > 0.5 && center.x > -0.5 && center.x < 0.5 && center.y > -0.5 && center.y < 0.5; - }); - - bcWall->findBoundaryCells(meshAdapter, true, [&](Vec3 center) { - return center.x < -0.5 || center.x > 0.5 || center.y < -0.5 || center.y > 0.5 || center.z < -0.5; - }); - - dataBase->boundaryConditions.push_back(bcLid); - dataBase->boundaryConditions.push_back(bcWall); - - ////////////////////////////////////////////////////////////////////////// - // set initial condition and upload mesh and initial condition to GPGPU - ////////////////////////////////////////////////////////////////////////// - - dataBase->setMesh(meshAdapter); + ////////////////////////////////////////////////////////////////////////// + // compute parameters in lattice units + ////////////////////////////////////////////////////////////////////////// - Initializer::interpret(dataBase, [&](Vec3 cellCenter) -> ConservedVariables { - return toConservedVariables(PrimitiveVariables(rho, 0.0, 0.0, 0.0, lambda), parameters.K); - }); + const real velocityLB = velocity * dt / dx; // LB units - dataBase->copyDataHostToDevice(); + const real vx = velocityLB / sqrt(2.0); // LB units + const real vy = velocityLB / sqrt(2.0); // LB units - Initializer::initializeDataUpdate(dataBase); + const real viscosityLB = nx * velocityLB / Re; // LB units - VtkWriter::write(dataBase, parameters, outputPath + "/" + simulationName + "_0"); + *logging::out << logging::Logger::INFO_HIGH << "velocity [dx/dt] = " << velocityLB << " \n"; + *logging::out << logging::Logger::INFO_HIGH << "viscosity [dx^2/dt] = " << viscosityLB << "\n"; - ////////////////////////////////////////////////////////////////////////// - // set analyzers - ////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////// + // set parameters + ////////////////////////////////////////////////////////////////////////// - CupsAnalyzer cupsAnalyzer(dataBase, false, 60.0, true, 10000); + para->setOutputPath(outputPath); + para->setOutputPrefix(simulationName); - ConvergenceAnalyzer convergenceAnalyzer(dataBase, 10000); + para->setPathAndFilename(para->getOutputPath() + "/" + para->getOutputPrefix()); - cupsAnalyzer.start(); + para->setPrintFiles(true); - ////////////////////////////////////////////////////////////////////////// - // run simulation - ////////////////////////////////////////////////////////////////////////// + para->setVelocityLB(velocityLB); + para->setViscosityLB(viscosityLB); - for (uint iter = 1; iter <= timeStepEnd; iter++) { - TimeStepping::nestedTimeStep(dataBase, parameters, 0); + para->setVelocityRatio(velocity / velocityLB); - if (iter % timeStepOut == 0) { - dataBase->copyDataDeviceToHost(); + para->setTimestepOut(timeStepOut); + para->setTimestepEnd(timeStepEnd); - VtkWriter::write(dataBase, parameters, outputPath + "/" + simulationName + "_" + std::to_string(iter)); - } + ////////////////////////////////////////////////////////////////////////// + // set boundary conditions + ////////////////////////////////////////////////////////////////////////// - int crashCellIndex = dataBase->getCrashCellIndex(); - if (crashCellIndex >= 0) { - *logging::out << logging::Logger::LOGGER_ERROR - << "Simulation crashed at CellIndex = " << crashCellIndex << "\n"; - dataBase->copyDataDeviceToHost(); - VtkWriter::write(dataBase, parameters, outputPath + "/" + simulationName + "_" + std::to_string(iter)); + gridBuilder->setNoSlipBoundaryCondition(SideType::PX); + gridBuilder->setNoSlipBoundaryCondition(SideType::MX); + gridBuilder->setNoSlipBoundaryCondition(SideType::PY); + gridBuilder->setNoSlipBoundaryCondition(SideType::MY); + gridBuilder->setVelocityBoundaryCondition(SideType::PZ, vx, vy, 0.0); + gridBuilder->setNoSlipBoundaryCondition(SideType::MZ); + gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, 0.0, 0.0); - break; - } + ////////////////////////////////////////////////////////////////////////// + // set copy mesh to simulation + ////////////////////////////////////////////////////////////////////////// - dataBase->getCrashCellIndex(); + SPtr<CudaMemoryManager> cudaMemoryManager = CudaMemoryManager::make(para); - cupsAnalyzer.run(iter, parameters.dt); + SPtr<GridProvider> gridGenerator = GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager); - convergenceAnalyzer.run(iter); - } - } - } catch (const std::bad_alloc& e) { + ////////////////////////////////////////////////////////////////////////// + // run simulation + ////////////////////////////////////////////////////////////////////////// + Simulation sim; + SPtr<FileWriter> fileWriter = SPtr<FileWriter>(new FileWriter()); + sim.init(para, gridGenerator, fileWriter, cudaMemoryManager); + sim.run(); + sim.free(); + + } catch (const std::bad_alloc &e) { *logging::out << logging::Logger::LOGGER_ERROR << "Bad Alloc:" << e.what() << "\n"; - } catch (const std::exception& e) { - + } catch (const std::exception &e) { *logging::out << logging::Logger::LOGGER_ERROR << e.what() << "\n"; } catch (std::string &s) { - *logging::out << logging::Logger::LOGGER_ERROR << s << "\n"; } catch (...) { *logging::out << logging::Logger::LOGGER_ERROR << "Unknown exception!\n";