Flow Around A Sphere (Uniform Grid Version)

View on GitLab

This example application simulates the flow behavior around a spherical obstacle in a channel.

Code Generation

The numerical kernels involved in the application are generated by the script FlowAroundSphereExample.py:

# 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/>.
 
from argparse import ArgumentParser
from dataclasses import replace
 
from lbmpy import (
    LBStencil,
    Stencil,
    LBMConfig,
    Method,
)
from lbmpy.boundaries import NoSlip, UBB, SimpleExtrapolationOutflow, QuadraticBounceBack
 
import sympy as sp
import pystencils as ps
 
from pystencilssfg import SourceFileGenerator
from sweepgen import Sweep, get_build_config
from sweepgen.boundaries import SparseBoundary
from sweepgen.prefabs import LbmBulk
 
from sweepgen.build_config import DEBUG
 
DEBUG.use_cpu_default()
 
with SourceFileGenerator(keep_unknown_argv=True) as sfg:
    Sweep.use_v8core_fields()
    sfg.namespace("FlowAroundSphereExample::gen")
 
    stencil = LBStencil(Stencil.D3Q19)
    omega = sp.symbols("omega")
    inflow_vel = sp.symbols("inflow_vel")
 
    base_lbm_config = LBMConfig(
        stencil=stencil,
        method=Method.CUMULANT,
        compressible=True,
        relaxation_rate=omega,
    )
 
    lbm_bulk = LbmBulk(sfg, "LBM", base_lbm_config)
    sfg.generate(lbm_bulk)
 
    rho, u = lbm_bulk.rho, lbm_bulk.u
 
    initial_state_assignments = [
        ps.Assignment(rho(), 1),
        ps.Assignment(u(0), inflow_vel)
    ]
 
    init_fields = Sweep("IntializeMacroFields", initial_state_assignments)
    sfg.generate(init_fields)
 
    noSlip = SparseBoundary(NoSlip(name="NoSlip"), lbm_bulk.lb_method, lbm_bulk.pdfs)
    sfg.generate(noSlip)
 
    qbb = SparseBoundary(QuadraticBounceBack(omega, name="QBB"), lbm_bulk.lb_method, lbm_bulk.pdfs)
    sfg.generate(qbb)
 
    inflow_velocity = (inflow_vel, 0, 0)
    ubb = SparseBoundary(UBB(inflow_velocity, name="UBB"), lbm_bulk.lb_method, lbm_bulk.pdfs)
    sfg.generate(ubb)
 
    outflow = SparseBoundary(SimpleExtrapolationOutflow((1, 0, 0), stencil, name="Outflow"), lbm_bulk.lb_method, lbm_bulk.pdfs)
    sfg.generate(outflow)

Application Frame

The simulation app itself is implemented in FlowAroundSphereExample.cpp:

//======================================================================================================================
//
//  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/>.
//
//
//======================================================================================================================
 
#include "core/all.h"
#include "blockforest/all.h"
#include "geometry/all.h"
#include "lbm/lattice_model/CollisionModel.h"
#include "stencil/all.h"
#include "field/all.h"
#include "timeloop/all.h"
 
#include "blockforest/communication/UniformBufferedScheme.h"
#include "field/communication/StencilRestrictedPackInfo.h"
#include "domain_decomposition/SharedSweep.h"
 
#include "walberla/v8/HaloExchange.hpp"
#include "walberla/v8/Memory.hpp"
#include "walberla/v8/Sweep.hpp"
 
#include "gen/FlowAroundSphereExample.hpp"
 
namespace FlowAroundSphere
{
    using namespace walberla;
    using namespace walberla::v8;
 
    using MemoryTag = memtag::automatic;
 
    using ScalarField_T = memory::Field<real_t, 1, MemoryTag>;
    using VectorField_T = memory::Field<real_t, 3, MemoryTag>;
 
    using LbStencil = stencil::D3Q19;
    using PdfField_T = memory::Field<real_t, LbStencil::Q, MemoryTag>;
 
 
    FlowAroundSphereExample::gen::QBBData sqSignedDistanceToSphere(geometry::Sphere sphere, Vector3<real_t> point) {
       real_t distance = sqrt(pow(sphere.midpoint()[0] - point[0], 2) + pow(sphere.midpoint()[1] - point[1], 2) + pow(sphere.midpoint()[2] - point[2], 2)) - sphere.radius();
       real_t sign = 1.0 ? distance >= 0 : -1.0;
       return {sign * distance * distance};
    }
 
    void run(int argc, char **argv)
    {
        Environment env{argc, argv};
        auto config = env.config();
 
        auto blocks = blockforest::createUniformBlockGridFromConfig(config);
        auto dx = blocks->dx();
        auto domainAABB = blocks->getDomain();
        WALBERLA_LOG_INFO_ON_ROOT("dx " << dx << " domain " << domainAABB)
        PdfField_T pdfs{ *blocks };
        ScalarField_T rho{ *blocks };
        VectorField_T u{ *blocks };
        
        Config::BlockHandle simParams = config->getBlock("Parameters");
        const real_t reynoldsNumber{ simParams.getParameter< real_t >("re") };
        const real_t vel_inflow = simParams.getParameter< real_t >("u_max");
        const real_t refLen = domainAABB.xSize()*0.1;
        const real_t viscosity = vel_inflow * refLen / reynoldsNumber;
        real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
 
        // Initialize macroscopic fields and PDF field
        FlowAroundSphereExample::gen::IntializeMacroFields initFields{rho, u, vel_inflow};
        FlowAroundSphereExample::gen::LBM::InitPdfs lbInit{pdfs, rho, u};
 
        for (auto &b : *blocks)
        {
           initFields(&b);
           lbInit(&b);
        }
 
        //  Set up LB stream/collide sweep
        auto streamCollide = std::make_shared< FlowAroundSphereExample::gen::LBM::StreamCollide >(pdfs, rho, u, omega);
 
        //  Set up ghost layer communication
        auto haloExchange = HaloExchange::create< LbStencil, MemoryTag >(blocks)
                               .sync(halo_exchange::streamPullSync< LbStencil >(pdfs)).makeShared();
 
        //  Set up boundary conditions
        geometry::Sphere sphere(Vector3<real_t> (domainAABB.xSize()*0.35,domainAABB.ySize()*0.5,domainAABB.zSize()*0.5), refLen*0.5);
 
        auto inflowLinks = [&](auto link) -> bool {
           blocks->transformBlockLocalToGlobalCell(link.wallCell, link.block);
           return link.wallCell.x() < blocks->getDomainCellBB().xMin();
        };
 
        auto outflowLinks = [&](auto link) -> bool {
           blocks->transformBlockLocalToGlobalCell(link.wallCell, link.block);
           return link.wallCell.x() > blocks->getDomainCellBB().xMax();
        };
 
        auto sideWallLinks = [&](auto link) -> bool {
           blocks->transformBlockLocalToGlobalCell(link.wallCell, link.block);
           return link.wallCell.y() < blocks->getDomainCellBB().yMin() || link.wallCell.y() > blocks->getDomainCellBB().yMax() ||
                  link.wallCell.z() < blocks->getDomainCellBB().zMin() || link.wallCell.z() > blocks->getDomainCellBB().zMax();
        };
 
        auto sphereLinks = [&](auto link) -> std::optional< FlowAroundSphereExample::gen::QBBData >  {
           blocks->transformBlockLocalToGlobalCell(link.wallCell, link.block);
           blocks->transformBlockLocalToGlobalCell(link.fluidCell, link.block);
           auto cellCenterSolid = blocks->getCellCenter(link.wallCell);
           auto cellCenterFluid = blocks->getCellCenter(link.fluidCell);
           if (contains(sphere, cellCenterSolid) && !contains(sphere, cellCenterFluid))
              return sqSignedDistanceToSphere(sphere, cellCenterFluid);
           else
              return std::nullopt;
        };
 
 
        auto qbb = FlowAroundSphereExample::gen::QBBFactory{ blocks, pdfs, omega }.fromLinks(sphereLinks);
        auto noSlip = FlowAroundSphereExample::gen::NoSlipFactory{ blocks, pdfs }.fromLinks(sideWallLinks);
        auto inflow = FlowAroundSphereExample::gen::UBBFactory{ blocks, pdfs, vel_inflow }.fromLinks(inflowLinks);
        auto outflow = FlowAroundSphereExample::gen::OutflowFactory{ blocks, pdfs }.fromLinks(outflowLinks);
        
        //  Timeloop
        const uint_t numTimesteps{simParams.getParameter<uint_t>("timesteps")};
        SweepTimeloop loop{blocks->getBlockStorage(), numTimesteps};
 
        loop.add() << Sweep(makeSharedSweep(streamCollide)) << AfterFunction(SharedFunctor(haloExchange));
        loop.add() << Sweep(noSlip);
        loop.add() << Sweep(inflow);
        loop.add() << Sweep(outflow);
        loop.add() << Sweep(qbb);
 
 
        RemainingTimeLogger logger{numTimesteps};
        loop.addFuncAfterTimeStep(logger);
 
        //  VTK Output
        Config::BlockHandle outputParams = config->getBlock("Output");
 
        const uint_t vtkWriteFrequency = outputParams.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);
 
            auto densityWriter = make_shared<memory::FieldVtkWriter<ScalarField_T, float32>>(rho, "density");
            vtkOutput->addCellDataWriter(densityWriter);
 
            auto velWriter = make_shared<memory::FieldVtkWriter<VectorField_T, float32>>(u, "velocity");
            vtkOutput->addCellDataWriter(velWriter);
 
            loop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
        }
 
        //  Run the Simulation
        WALBERLA_LOG_INFO_ON_ROOT("Commencing simulation with " << numTimesteps << " timesteps")
        loop.run();
    }
}
 
int main(int argc, char **argv)
{
    FlowAroundSphere::run(argc, argv);
    return EXIT_SUCCESS;
}