From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 24 Jun 2019 14:43:36 +0200 Subject: Initialize at openlb-1-3 --- examples/particles/bifurcation3d/bifurcation3d.stl | Bin 0 -> 79184 bytes .../particles/bifurcation3d/eulerEuler/Makefile | 105 ++++ .../bifurcation3d/eulerEuler/bifurcation3d.cpp | 420 ++++++++++++++ .../bifurcation3d/eulerEuler/definitions.mk | 30 + .../particles/bifurcation3d/eulerEuler/module.mk | 29 + .../particles/bifurcation3d/eulerLagrange/Makefile | 105 ++++ .../bifurcation3d/eulerLagrange/bifurcation3d.cpp | 515 ++++++++++++++++ .../bifurcation3d/eulerLagrange/definitions.mk | 30 + .../bifurcation3d/eulerLagrange/module.mk | 29 + examples/particles/dkt2d/Makefile | 115 ++++ examples/particles/dkt2d/definitions.mk | 30 + examples/particles/dkt2d/dkt.p | 22 + examples/particles/dkt2d/dkt2d.cpp | 267 +++++++++ examples/particles/dkt2d/module.mk | 29 + examples/particles/magneticParticles3d/Makefile | 105 ++++ .../particles/magneticParticles3d/definitions.mk | 30 + .../magneticParticles3d/magneticParticles3d.cpp | 644 +++++++++++++++++++++ .../magneticParticles3d/magneticParticles3d.xml | 47 ++ examples/particles/magneticParticles3d/module.mk | 29 + examples/particles/settlingCube3d/Makefile | 105 ++++ examples/particles/settlingCube3d/definitions.mk | 30 + examples/particles/settlingCube3d/module.mk | 29 + .../particles/settlingCube3d/settlingCube3D.cpp | 294 ++++++++++ 23 files changed, 3039 insertions(+) create mode 100644 examples/particles/bifurcation3d/bifurcation3d.stl create mode 100644 examples/particles/bifurcation3d/eulerEuler/Makefile create mode 100644 examples/particles/bifurcation3d/eulerEuler/bifurcation3d.cpp create mode 100644 examples/particles/bifurcation3d/eulerEuler/definitions.mk create mode 100644 examples/particles/bifurcation3d/eulerEuler/module.mk create mode 100644 examples/particles/bifurcation3d/eulerLagrange/Makefile create mode 100644 examples/particles/bifurcation3d/eulerLagrange/bifurcation3d.cpp create mode 100644 examples/particles/bifurcation3d/eulerLagrange/definitions.mk create mode 100644 examples/particles/bifurcation3d/eulerLagrange/module.mk create mode 100644 examples/particles/dkt2d/Makefile create mode 100644 examples/particles/dkt2d/definitions.mk create mode 100644 examples/particles/dkt2d/dkt.p create mode 100644 examples/particles/dkt2d/dkt2d.cpp create mode 100644 examples/particles/dkt2d/module.mk create mode 100644 examples/particles/magneticParticles3d/Makefile create mode 100644 examples/particles/magneticParticles3d/definitions.mk create mode 100644 examples/particles/magneticParticles3d/magneticParticles3d.cpp create mode 100644 examples/particles/magneticParticles3d/magneticParticles3d.xml create mode 100644 examples/particles/magneticParticles3d/module.mk create mode 100755 examples/particles/settlingCube3d/Makefile create mode 100755 examples/particles/settlingCube3d/definitions.mk create mode 100644 examples/particles/settlingCube3d/module.mk create mode 100644 examples/particles/settlingCube3d/settlingCube3D.cpp (limited to 'examples/particles') diff --git a/examples/particles/bifurcation3d/bifurcation3d.stl b/examples/particles/bifurcation3d/bifurcation3d.stl new file mode 100644 index 0000000..929e3a4 Binary files /dev/null and b/examples/particles/bifurcation3d/bifurcation3d.stl differ diff --git a/examples/particles/bifurcation3d/eulerEuler/Makefile b/examples/particles/bifurcation3d/eulerEuler/Makefile new file mode 100644 index 0000000..a953954 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerEuler/Makefile @@ -0,0 +1,105 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + +########################################################################### +## definitions + +include definitions.mk + +include $(ROOT)/global.mk + +OBJECTS := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.o)) +DEPS := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.d)) + +########################################################################### +## all + +all : depend compile link + + +########################################################################### +## dependencies + +depend : $(DEPS) + +$(PWD)/%.d : %.cpp + @echo Create dependencies for $< + @$(SHELL) -ec '$(CXX) -M $(CXXFLAGS) $(IDIR) $< \ + | sed -e "s!$*\.o!$(PWD)\/$*\.o!1" > .tmpfile; \ + cp -f .tmpfile $@;' + +########################################################################### +## compile + +compile : $(OBJECTS) + +$(PWD)/%.o: %.cpp + @echo Compile $< + $(CXX) $(CXXFLAGS) $(IDIR) -c $< -o $@ + +########################################################################### +## clean + +clean : cleanrub cleanobj cleandep + +cleanrub: + @echo Clean rubbish files + @rm -f *~ core .tmpfile tmp/*.* $(OUTPUT) + @rm -f tmp/vtkData/*.* tmp/vtkData/data/*.* tmp/imageData/*.* tmp/gnuplotData/*.* tmp/gnuplotData/data/*.* + +cleanobj: + @echo Clean object files + @rm -f $(OBJECTS) + +cleandep: + @echo Clean dependencies files + @rm -f $(DEPS) + +cleanbuild: + @cd $(ROOT); \ + $(MAKE) cleanlib; + +########################################################################### +## update lib + +$(ROOT)/$(LIBDIR)/lib$(LIB).a : + @cd $(ROOT); \ + $(MAKE) all + +########################################################################### +## link + +link: $(OUTPUT) + +$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a + @echo Link $@ + $(CXX) $(foreach file, $(SRC), $(file:.cpp=.o)) $(LDFLAGS) -L$(ROOT)/$(LIBDIR) -l$(LIB) -lz -o $@ + +########################################################################### +## include dependencies + +ifneq "$(strip $(wildcard *.d))" "" + include $(foreach file,$(DEPS),$(file)) +endif + +########################################################################### +########################################################################### diff --git a/examples/particles/bifurcation3d/eulerEuler/bifurcation3d.cpp b/examples/particles/bifurcation3d/eulerEuler/bifurcation3d.cpp new file mode 100644 index 0000000..69ac959 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerEuler/bifurcation3d.cpp @@ -0,0 +1,420 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2011-2016 Robin Trunk, Mathias J. Krause + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * This program 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 2 + * of the License, or (at your option) any later version. + * + * This program 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 this program; if not, write to the Free + * Software Foundation, Inc., 51 Franklin Street, Fifth Floor, + * Boston, MA 02110-1301, USA. + */ + +/* bifurcation3d.cpp: + * This example examines a steady particulate flow past a bifurcation. At the inlet, + * an inflow condition with grad_n u = 0 and rho = 1 is implemented. + * At both outlets, a Poiseuille profile is imposed on the velocity. + * After a start time, particles are put into the bifurcation by imposing + * a inflow condition with rho = 1 on the second euler phase at the inlet. + * The particles are simulated as a continuum with a advection-diffusion equation + * and experience a stokes drag force. + * + * A publication using the same geometry can be found here: + * http://link.springer.com/chapter/10.1007/978-3-642-36961-2_5 + * * + */ + +#include "olb3D.h" +#include "olb3D.hh" // use only generic version! + +using namespace std; +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; + +typedef double T; +#define NSDESCRIPTOR D3Q19<> +#define ADDESCRIPTOR D3Q7 + +const T Re = 50; // Reynolds number +const int N = 19; // resolution of the model +const int iTperiod = 100; // amount of timesteps when new boundary conditions are reset and results are visualized +const T diffusion = 1.e-6; // diffusion coefficient for advection-diffusion equation +const T radius = 1.5e-04; // particles radius +const T partRho = 998.2; // particles density +const T maxPhysT = 10.; // max. simulation time in s, SI unit + +// center of inflow and outflow regions [m] +Vector inletCenter( T(), T(), 0.0786395 ); +Vector outletCenter0( -0.0235929682287551, -0.000052820468762797, -0.021445708949909 ); +Vector outletCenter1( 0.0233643529416147, 0.00000212439067050152, -0.0211994104877918 ); + +// radii of inflow and outflow regions [m] +T inletRadius = 0.00999839; +T outletRadius0 = 0.007927; +T outletRadius1 = 0.00787134; + +// normals of inflow and outflow regions +Vector inletNormal( T(), T(), T( -1 ) ); +Vector outletNormal0( 0.505126, -0.04177, 0.862034 ); +Vector outletNormal1( -0.483331, -0.0102764, 0.875377 ); + +void prepareGeometry( UnitConverter const& converter, + IndicatorF3D& indicator, STLreader& stlReader, + SuperGeometry3D& superGeometry ) +{ + + OstreamManager clout( std::cout, "prepareGeometry" ); + + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename( 0, 2, indicator ); + superGeometry.rename( 2, 1, stlReader ); + + superGeometry.clean(); + + // rename the material at the inlet + IndicatorCircle3D inletCircle( inletCenter, inletNormal, inletRadius ); + IndicatorCylinder3D inlet( inletCircle, 2 * converter.getConversionFactorLength() ); + superGeometry.rename( 2, 3, 1, inlet ); + + // rename the material at the outlet0 + IndicatorCircle3D outletCircle0( outletCenter0, outletNormal0, 0.95*outletRadius0 ); + IndicatorCylinder3D outlet0( outletCircle0, 4 * converter.getConversionFactorLength() ); + superGeometry.rename( 2, 4, outlet0 ); + + // rename the material at the outlet1 + IndicatorCircle3D outletCircle1( outletCenter1, outletNormal1, 0.95*outletRadius1 ); + IndicatorCylinder3D outlet1( outletCircle1, 4 * converter.getConversionFactorLength() ); + superGeometry.rename( 2, 5, outlet1 ); + + superGeometry.clean(); + superGeometry.innerClean( true ); + superGeometry.checkForErrors(); + + superGeometry.print(); + + clout << "Prepare Geometry ... OK" << std::endl; + return; +} + +void prepareLattice( SuperLattice3D& sLatticeNS, + SuperLattice3D& sLatticeAD, + UnitConverter const& converter, + Dynamics& bulkDynamics, + Dynamics& bulkDynamicsAD, + sOnLatticeBoundaryCondition3D& bc, + sOnLatticeBoundaryCondition3D& bcAD, + SuperGeometry3D& superGeometry, + T omegaAD ) +{ + + OstreamManager clout( std::cout, "prepareLattice" ); + clout << "Prepare Lattice ..." << std::endl; + + const T omega = converter.getLatticeRelaxationFrequency(); + + // Material=0 --> do nothing + sLatticeNS.defineDynamics( superGeometry, 0, &instances::getNoDynamics() ); + sLatticeAD.defineDynamics( superGeometry, 0, &instances::getNoDynamics() ); + + // Material=1 --> bulk dynamics + // Material=3 --> bulk dynamics (inflow) + auto inflowIndicator = superGeometry.getMaterialIndicator({1, 3}); + sLatticeNS.defineDynamics( inflowIndicator, &bulkDynamics ); + sLatticeAD.defineDynamics( inflowIndicator, &bulkDynamicsAD ); + + // Material=2 --> bounce-back / + sLatticeNS.defineDynamics( superGeometry, 2, &instances::getBounceBack() ); + sLatticeAD.defineDynamics( superGeometry, 2, &instances::getZeroDistributionDynamics() ); + + // Material=4,5 -->bulk dynamics / do-nothing (outflow) + auto outflowIndicator = superGeometry.getMaterialIndicator({4, 5}); + sLatticeNS.defineDynamics( outflowIndicator, &bulkDynamics ); + sLatticeAD.defineDynamics( outflowIndicator, &instances::getNoDynamics() ); + + // Setting of the boundary conditions + bc.addPressureBoundary( superGeometry, 3, omega ); + bc.addVelocityBoundary( outflowIndicator, omega ); + bcAD.addZeroDistributionBoundary( superGeometry, 2 ); + bcAD.addTemperatureBoundary( superGeometry, 3, omegaAD ); + bcAD.addConvectionBoundary( outflowIndicator ); + bcAD.addExtFieldBoundary( superGeometry.getMaterialIndicator({2, 3, 4, 5}), ADDESCRIPTOR::index() ); + + // Initial conditions + AnalyticalConst3D rho1( 1. ); + AnalyticalConst3D rho0( 1.e-8 ); + std::vector velocity( 3,T() ); + AnalyticalConst3D u0( velocity ); + + // Initialize all values of distribution functions to their local equilibrium + sLatticeNS.defineRhoU( superGeometry.getMaterialIndicator({1, 2, 3, 4, 5}), rho1, u0 ); + sLatticeNS.iniEquilibrium( superGeometry.getMaterialIndicator({1, 2, 3, 4, 5}), rho1, u0 ); + sLatticeAD.defineRho( superGeometry, 3, rho1 ); + sLatticeAD.iniEquilibrium( superGeometry.getMaterialIndicator({1, 2, 4, 5}), rho0, u0 ); + + // Lattice initialize + sLatticeNS.initialize(); + sLatticeAD.initialize(); + + clout << "Prepare Lattice ... OK" << std::endl; + return; +} + +void setBoundaryValues( SuperLattice3D& sLatticeNS, + UnitConverter const& converter, int iT, + T maxPhysT, SuperGeometry3D& superGeometry ) +{ + + OstreamManager clout( std::cout, "setBoundaryValues" ); + + // No of time steps for smooth start-up + int iTmaxStart = converter.getLatticeTime( 0.8*maxPhysT ); + // Set inflow velocity + T maxVelocity = converter.getCharLatticeVelocity() * 3. / 4. * std::pow( + inletRadius, 2 ) / std::pow( outletRadius0, 2 ); + if ( iT % iTperiod == 0 ) { + if ( iT <= iTmaxStart ) { + SinusStartScale startScale( iTmaxStart, T( 1 ) ); + int iTvec[1] = { iT }; + T frac[1] = { T( 0 ) }; + startScale( frac, iTvec ); + maxVelocity *= frac[0]; + } + + CirclePoiseuille3D poiseuilleU4( outletCenter0[0], outletCenter0[1], + outletCenter0[2], outletNormal0[0], + outletNormal0[1], outletNormal0[2], + outletRadius0 * 0.95, -maxVelocity ); + + CirclePoiseuille3D poiseuilleU5( outletCenter1[0], outletCenter1[1], + outletCenter1[2], outletNormal1[0], + outletNormal1[1], outletNormal1[2], + outletRadius1 * 0.95, -maxVelocity ); + + sLatticeNS.defineU( superGeometry, 4, poiseuilleU4 ); + sLatticeNS.defineU( superGeometry, 5, poiseuilleU5 ); + } +} + +void getResults( SuperLattice3D& sLatticeNS, + SuperLattice3D& sLatticeAD, + UnitConverter const& converter, int iT, + SuperGeometry3D& superGeometry, + Timer& timer ) +{ + + OstreamManager clout( std::cout, "getResults" ); + SuperVTMwriter3D vtmWriter( "bifurcation3d_fluid" ); + SuperVTMwriter3D vtmWriterAD( "bifurcation3d_particle" ); + SuperLatticePhysVelocity3D velocity( sLatticeNS, converter ); + SuperLatticeVelocity3D latticeVelocity( sLatticeNS); + + SuperLatticePhysPressure3D pressure( sLatticeNS, converter ); + SuperLatticeDensity3D particles( sLatticeAD ); + SuperLatticePhysExternal3D extField(sLatticeAD, + converter.getConversionFactorVelocity(), + ADDESCRIPTOR::index(), + ADDESCRIPTOR::size()); + + vtmWriter.addFunctor( velocity ); + vtmWriter.addFunctor( pressure ); + vtmWriterAD.addFunctor( particles ); + vtmWriterAD.addFunctor( extField ); + + if ( iT == 0 ) { + SuperLatticeGeometry3D geometry( sLatticeNS, superGeometry ); + SuperLatticeCuboid3D cuboid( sLatticeNS ); + SuperLatticeRank3D rank( sLatticeNS ); + vtmWriter.write( geometry ); + vtmWriter.write( cuboid ); + vtmWriter.write( rank ); + vtmWriter.createMasterFile(); + vtmWriterAD.createMasterFile(); + + // Print some output of the chosen simulation setup + clout << "N=" << N << "; maxTimeSteps=" << converter.getLatticeTime( maxPhysT ) + << "; noOfCuboid=" << superGeometry.getCuboidGeometry().getNc() << "; Re=" << Re + << "; St=" << ( 2.*partRho*radius*radius*converter.getCharPhysVelocity() ) / ( 9.*converter.getPhysViscosity()*converter.getPhysDensity()*converter.getCharPhysLength() ) + << std::endl; + } + + if ( iT % iTperiod == 0 ) { + // Writes the vtk files + vtmWriter.write( iT ); + vtmWriterAD.write( iT ); + + // GIF Writer + SuperEuklidNorm3D normVel( velocity ); + HyperplaneLattice3D gifLattice( + superGeometry.getCuboidGeometry(), + Hyperplane3D() + .centeredIn(superGeometry.getCuboidGeometry().getMotherCuboid()) + .normalTo({0, -1, 0}), + 600); + BlockReduction3D2D planeReductionVelocity( normVel, gifLattice, BlockDataSyncMode::ReduceOnly ); + BlockReduction3D2D planeReductionParticles( particles, gifLattice, BlockDataSyncMode::ReduceOnly ); + // write output as JPEG + heatmap::write(planeReductionVelocity, iT); + heatmap::write(planeReductionParticles, iT); + + // Writes output on the console + timer.update( iT ); + timer.printStep(); + sLatticeNS.getStatistics().print( iT, iT * converter.getCharLatticeVelocity() / T(converter.getResolution()) ); + + // preparation for flux computations + const std::vector materials = { 1, 3, 4, 5 }; + IndicatorCircle3D inlet( inletCenter + 2. * converter.getConversionFactorLength() * inletNormal, + inletNormal, + inletRadius + 2. * converter.getConversionFactorLength() ); + IndicatorCircle3D outlet0( outletCenter0 + 2. * converter.getConversionFactorLength() * outletNormal0, + outletNormal0, + outletRadius0 + 2. * converter.getConversionFactorLength() ); + IndicatorCircle3D outlet1( outletCenter1 + 2. * converter.getConversionFactorLength() * outletNormal1, + outletNormal1, + outletRadius1 + 2. * converter.getConversionFactorLength() ); + + // Flux of the fluid at the inlet and outlet regions + SuperPlaneIntegralFluxVelocity3D vFluxInflow( sLatticeNS, converter, superGeometry, inlet, materials ); + vFluxInflow.print( "inflow", "ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxInflow( sLatticeNS, converter, superGeometry, inlet, materials ); + pFluxInflow.print( "inflow", "N", "Pa" ); + SuperPlaneIntegralFluxVelocity3D vFluxOutflow0( sLatticeNS, converter, superGeometry, outlet0, materials ); + vFluxOutflow0.print( "outflow0", "ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxOutflow0( sLatticeNS, converter, superGeometry, outlet0, materials ); + pFluxOutflow0.print( "outflow0", "N", "Pa" ); + SuperPlaneIntegralFluxVelocity3D vFluxOutflow1( sLatticeNS, converter, superGeometry, outlet1, materials ); + vFluxOutflow1.print( "outflow1", "ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxOutflow1( sLatticeNS, converter, superGeometry, outlet1, materials ); + pFluxOutflow1.print( "outflow1", "N", "Pa" ); + + int input[4] = {0}; + T mFlux[5] = {0.}, mFlux0[5] = {0.}, mFlux1[5] = {0.}; + // Flux of particles at the inlet and outlet regions: Inflow, Outflow0 and Outlfow1 + SuperPlaneIntegralFluxMass3D mFluxInflow(latticeVelocity,particles, + superGeometry, converter.getConversionFactorMass(), + converter.getConversionFactorTime(), inlet, materials); + SuperPlaneIntegralFluxMass3D mFluxOutflow0(latticeVelocity, particles, + superGeometry, converter.getConversionFactorMass(), + converter.getConversionFactorTime(),outlet0, materials); + SuperPlaneIntegralFluxMass3D mFluxOutflow1(latticeVelocity, particles, + superGeometry, converter.getConversionFactorMass(), + converter.getConversionFactorTime(), outlet1, materials); + + mFluxInflow( mFlux,input ); + mFluxOutflow0( mFlux0,input ); + mFluxOutflow1( mFlux1,input ); + + // Since more diffusion is added to ensure stability the computed escaperate falls short of the real value, + // therefore it is scaled by the factor 1.4 computed by a simulation without drag force. This value is computed + // for this specific setup. For further information see R.Trunk, T.Henn, W.Dörfler, H.Nirschl, M.J.Krause, + // "Inertial Dilute Particulate Fluid Flow Simulations with an Euler-Euler Lattice Boltzmann Method" + T escr = -( mFlux0[0]+mFlux1[0] )/mFlux[0]*1.4; + clout << "escapeRate=" << escr << "; captureRate="<< 1-escr << std::endl; + } +} + +int main( int argc, char* argv[] ) +{ + + // === 1st Step: Initialization === + + olbInit( &argc, &argv ); + singleton::directories().setOutputDir( "./tmp/" ); + OstreamManager clout( std::cout, "main" ); + + UnitConverterFromResolutionAndRelaxationTime const converter( + int {N}, // resolution: number of voxels per charPhysL + (T) 0.557646, // latticeRelaxationTime: relaxation time, have to be greater than 0.5! + (T) inletRadius*2., // charPhysLength: reference length of simulation geometry + (T) Re*1.5e-5/( inletRadius*2 ), // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__ + (T) 1.5e-5, // physViscosity: physical kinematic viscosity in __m^2 / s__ + (T) 1.225 // physDensity: physical density in __kg / m^3__ + ); + // Prints the converter log as console output + converter.print(); + // Writes the converter log in a file + converter.write("bifurcation3d"); + + // compute relaxation parameter to solve the advection-diffusion equation in the lattice Boltzmann context + T omegaAD = converter.getLatticeRelaxationFrequencyFromDiffusivity( diffusion ); + + // === 2nd Step: Prepare Geometry === + + STLreader stlReader( "../bifurcation3d.stl",converter.getConversionFactorLength() ); + IndicatorLayer3D extendedDomain( stlReader,converter.getConversionFactorLength() ); + + // Instantiation of an empty cuboidGeometry + int noOfCuboids = std::max( 20, singleton::mpi().getSize() ); + CuboidGeometry3D cuboidGeometry( extendedDomain, converter.getConversionFactorLength(), + noOfCuboids, "weight" ); + cuboidGeometry.printExtended(); + clout << "min / max ratio (volume) = " << (T) cuboidGeometry.getMinLatticeVolume() + / cuboidGeometry.getMaxLatticeVolume() << endl; + clout << "min / max ratio (weight) = " << (T) cuboidGeometry.getMinLatticeWeight() + / cuboidGeometry.getMaxLatticeWeight() << endl; + + // Instantiation of an empty loadBalancer + HeuristicLoadBalancer loadBalancer( cuboidGeometry ); + // Default instantiation of superGeometry + SuperGeometry3D superGeometry( cuboidGeometry, loadBalancer, 2 ); + + prepareGeometry( converter, extendedDomain, stlReader, superGeometry ); + + // === 3rd Step: Prepare Lattice === + SuperLattice3D sLatticeNS( superGeometry ); + SuperLattice3D sLatticeAD( superGeometry ); + SuperExternal3D sExternal( superGeometry, sLatticeAD, sLatticeAD.getOverlap() ); + + BGKdynamics bulkDynamics( converter.getLatticeRelaxationFrequency(), + instances::getBulkMomenta() ); + ParticleAdvectionDiffusionBGKdynamics bulkDynamicsAD ( omegaAD, + instances::getBulkMomenta() ); + + sOnLatticeBoundaryCondition3D sBoundaryCondition( sLatticeNS ); + createInterpBoundaryCondition3D( sBoundaryCondition ); + + sOnLatticeBoundaryCondition3D sBoundaryConditionAD( sLatticeAD ); + createAdvectionDiffusionBoundaryCondition3D( sBoundaryConditionAD ); + + prepareLattice( sLatticeNS, sLatticeAD, converter, bulkDynamics, bulkDynamicsAD, + sBoundaryCondition, sBoundaryConditionAD, superGeometry, omegaAD ); + + AdvectionDiffusionParticleCouplingGenerator3D coupling( + ADDESCRIPTOR::index() ); + + AdvDiffDragForce3D dragForce( converter,radius,partRho ); + coupling.addForce( dragForce ); + sLatticeNS.addLatticeCoupling( superGeometry, 1, coupling, sLatticeAD ); + + // === 4th Step: Main Loop with Timer === + Timer timer( converter.getLatticeTime( maxPhysT ), + superGeometry.getStatistics().getNvoxel() ); + timer.start(); + + for ( int iT = 0; iT <= converter.getLatticeTime( maxPhysT ); ++iT ) { + getResults( sLatticeNS, sLatticeAD, converter, iT, superGeometry, timer ); + setBoundaryValues( sLatticeNS, converter, iT, maxPhysT, superGeometry ); + sLatticeNS.executeCoupling(); + sExternal.communicate(); + sLatticeNS.collideAndStream(); + sLatticeAD.collideAndStream(); + } + timer.stop(); + timer.printSummary(); + +} diff --git a/examples/particles/bifurcation3d/eulerEuler/definitions.mk b/examples/particles/bifurcation3d/eulerEuler/definitions.mk new file mode 100644 index 0000000..3ad56e8 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerEuler/definitions.mk @@ -0,0 +1,30 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + + +########################################################################### +########################################################################### +## DEFINITIONS TO BE CHANGED + +ROOT := ../../../.. +SRC := bifurcation3d.cpp +OUTPUT := bifurcation3d diff --git a/examples/particles/bifurcation3d/eulerEuler/module.mk b/examples/particles/bifurcation3d/eulerEuler/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerEuler/module.mk @@ -0,0 +1,29 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2017 Markus Mohrhard +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + +current_dir := $(dir $(word $(words $(MAKEFILE_LIST)),$(MAKEFILE_LIST))) + +include global.mk +include rules.mk +include $(addsuffix definitions.mk, $(current_dir)) + +$(eval $(call sample,$(current_dir)$(OUTPUT),$(addprefix $(current_dir), $(SRC)))) diff --git a/examples/particles/bifurcation3d/eulerLagrange/Makefile b/examples/particles/bifurcation3d/eulerLagrange/Makefile new file mode 100644 index 0000000..a953954 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerLagrange/Makefile @@ -0,0 +1,105 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + +########################################################################### +## definitions + +include definitions.mk + +include $(ROOT)/global.mk + +OBJECTS := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.o)) +DEPS := $(foreach file, $(SRC), $(PWD)/$(file:.cpp=.d)) + +########################################################################### +## all + +all : depend compile link + + +########################################################################### +## dependencies + +depend : $(DEPS) + +$(PWD)/%.d : %.cpp + @echo Create dependencies for $< + @$(SHELL) -ec '$(CXX) -M $(CXXFLAGS) $(IDIR) $< \ + | sed -e "s!$*\.o!$(PWD)\/$*\.o!1" > .tmpfile; \ + cp -f .tmpfile $@;' + +########################################################################### +## compile + +compile : $(OBJECTS) + +$(PWD)/%.o: %.cpp + @echo Compile $< + $(CXX) $(CXXFLAGS) $(IDIR) -c $< -o $@ + +########################################################################### +## clean + +clean : cleanrub cleanobj cleandep + +cleanrub: + @echo Clean rubbish files + @rm -f *~ core .tmpfile tmp/*.* $(OUTPUT) + @rm -f tmp/vtkData/*.* tmp/vtkData/data/*.* tmp/imageData/*.* tmp/gnuplotData/*.* tmp/gnuplotData/data/*.* + +cleanobj: + @echo Clean object files + @rm -f $(OBJECTS) + +cleandep: + @echo Clean dependencies files + @rm -f $(DEPS) + +cleanbuild: + @cd $(ROOT); \ + $(MAKE) cleanlib; + +########################################################################### +## update lib + +$(ROOT)/$(LIBDIR)/lib$(LIB).a : + @cd $(ROOT); \ + $(MAKE) all + +########################################################################### +## link + +link: $(OUTPUT) + +$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a + @echo Link $@ + $(CXX) $(foreach file, $(SRC), $(file:.cpp=.o)) $(LDFLAGS) -L$(ROOT)/$(LIBDIR) -l$(LIB) -lz -o $@ + +########################################################################### +## include dependencies + +ifneq "$(strip $(wildcard *.d))" "" + include $(foreach file,$(DEPS),$(file)) +endif + +########################################################################### +########################################################################### diff --git a/examples/particles/bifurcation3d/eulerLagrange/bifurcation3d.cpp b/examples/particles/bifurcation3d/eulerLagrange/bifurcation3d.cpp new file mode 100644 index 0000000..43cbcb5 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerLagrange/bifurcation3d.cpp @@ -0,0 +1,515 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2011-2016 Thomas Henn, Mathias J. Krause, + * Marie-Luise Maier + * E-mail contact: info@openlb.net + * The most recent release of OpenLB can be downloaded at + * + * + * This program 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 2 + * of the License, or (at your option) any later version. + * + * This program 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 this program; if not, write to the Free + * Software Foundation, Inc., 51 Franklin Street, Fifth Floor, + * Boston, MA 02110-1301, USA. + */ + +/* bifurcation3d.cpp: + * This example examines a steady particulate flow past a bifurcation. At the inlet, + * an inflow condition with grad_n u = 0 and rho = 1 is implemented. + * At both outlets, a Poiseuille profile is imposed on the velocity. + * After a start time, particles are put into the bifurcation at the + * inlet and experience a stokes drag force. + * + * A publication using the same geometry can be found here: + * http://link.springer.com/chapter/10.1007/978-3-642-36961-2_5 + * * + */ + +#include "olb3D.h" +#ifndef OLB_PRECOMPILED // Unless precompiled version is used, +#include "olb3D.hh" // include full template code; +#endif + +using namespace std; +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; + +typedef double T; +#define DESCRIPTOR D3Q19<> +#define PARTICLE Particle3D + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +const T Re = 50; // Reynolds number +const int N = 19; // resolution of the model +const T radius = 1.5e-4; // particles radius +const T partRho = 998.2; //particles density + +const T fluidMaxPhysT = T( 5 ); // max. fluid simulation time in s, SI unit +const T particleMaxPhysT = T( 10 ); // max. particle simulation time in s, SI unit + +const int noOfParticles = 1000; // total number of inserted particles + +// center of inflow and outflow regions [m] +Vector inletCenter( T(), T(), 0.0786395 ); +Vector outletCenter0( -0.0235929682287551, -0.000052820468762797, + -0.021445708949909 ); +Vector outletCenter1( 0.0233643529416147, 0.00000212439067050152, + -0.0211994104877918 ); + +// radii of inflow and outflow regions [m] +T inletRadius = 0.00999839; +T outletRadius0 = 0.007927; +T outletRadius1 = 0.00787134; + +// normals of inflow and outflow regions +Vector inletNormal( T(), T(), T( -1 ) ); +Vector outletNormal0( 0.505126, -0.04177, 0.862034 ); +Vector outletNormal1( -0.483331, -0.0102764, 0.875377 ); + +void prepareGeometry( UnitConverter const& converter, + IndicatorF3D& indicator, STLreader& stlReader, + SuperGeometry3D& superGeometry ) +{ + + OstreamManager clout( std::cout, "prepareGeometry" ); + + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename( 0, 2, indicator ); + superGeometry.rename( 2, 1, stlReader ); + + superGeometry.clean(); + + // rename the material at the inlet + IndicatorCircle3D inletCircle( inletCenter, inletNormal, + inletRadius ); + IndicatorCylinder3D inlet( inletCircle, + 2 * converter.getConversionFactorLength() ); + superGeometry.rename( 2, 3, 1, inlet ); + + // rename the material at the outlet0 + IndicatorCircle3D outletCircle0( outletCenter0, outletNormal0, + 0.95 * outletRadius0 ); + IndicatorCylinder3D outlet0( outletCircle0, + 4 * converter.getConversionFactorLength() ); + superGeometry.rename( 2, 4, outlet0 ); + + // rename the material at the outlet1 + IndicatorCircle3D outletCircle1( outletCenter1, outletNormal1, + 0.95 * outletRadius1 ); + IndicatorCylinder3D outlet1( outletCircle1, + 4 * converter.getConversionFactorLength() ); + superGeometry.rename( 2, 5, outlet1 ); + + superGeometry.clean(); + superGeometry.innerClean( true ); + superGeometry.checkForErrors(); + + superGeometry.print(); + + clout << "Prepare Geometry ... OK" << std::endl; + return; +} + +void prepareLattice( SuperLattice3D& sLattice, + UnitConverter const& converter, Dynamics& + bulkDynamics, + sOnLatticeBoundaryCondition3D& bc, + sOffLatticeBoundaryCondition3D& offBc, + STLreader& stlReader, SuperGeometry3D& superGeometry ) +{ + + OstreamManager clout( std::cout, "prepareLattice" ); + clout << "Prepare Lattice ..." << std::endl; + + const T omega = converter.getLatticeRelaxationFrequency(); + + // Material=0 -->do nothing + sLattice.defineDynamics( superGeometry, 0, + &instances::getNoDynamics() ); + + // Material=1 -->bulk dynamics + sLattice.defineDynamics( superGeometry, 1, &bulkDynamics ); + + // Material=2 -->bounce back + sLattice.defineDynamics( superGeometry, 2, + &instances::getBounceBack() ); + + // Material=3 -->bulk dynamics (inflow) + sLattice.defineDynamics( superGeometry, 3, &bulkDynamics ); + + // Material=4 -->bulk dynamics (outflow) + sLattice.defineDynamics( superGeometry, 4, &bulkDynamics ); + sLattice.defineDynamics( superGeometry, 5, &bulkDynamics ); + + // Setting of the boundary conditions + bc.addPressureBoundary( superGeometry, 3, omega ); + bc.addVelocityBoundary( superGeometry, 4, omega ); + bc.addVelocityBoundary( superGeometry, 5, omega ); + + clout << "Prepare Lattice ... OK" << std::endl; + return; +} + +// Generates a slowly increasing sinuidal inflow for the first iTMax timesteps +void setBoundaryValues( SuperLattice3D& sLattice, + UnitConverter const& converter, int iT, T maxPhysT, + SuperGeometry3D& superGeometry ) +{ + + OstreamManager clout( std::cout, "setBoundaryValues" ); + + // No of time steps for smooth start-up + int iTmaxStart = converter.getLatticeTime( 0.8*maxPhysT ); + int iTperiod = 100; // amount of timesteps when new boundary conditions are reset + + if ( iT == 0 ) { + + AnalyticalConst3D rhoF( 1 ); + std::vector velocity( 3, T() ); + AnalyticalConst3D uF( velocity ); + + sLattice.iniEquilibrium( superGeometry, 1, rhoF, uF ); + sLattice.iniEquilibrium( superGeometry, 2, rhoF, uF ); + sLattice.iniEquilibrium( superGeometry, 3, rhoF, uF ); + sLattice.iniEquilibrium( superGeometry, 4, rhoF, uF ); + sLattice.iniEquilibrium( superGeometry, 5, rhoF, uF ); + + sLattice.defineRhoU( superGeometry, 1, rhoF, uF ); + sLattice.defineRhoU( superGeometry, 2, rhoF, uF ); + sLattice.defineRhoU( superGeometry, 3, rhoF, uF ); + sLattice.defineRhoU( superGeometry, 4, rhoF, uF ); + sLattice.defineRhoU( superGeometry, 5, rhoF, uF ); + + // Make the lattice ready for simulation + sLattice.initialize(); + } + + else if ( iT <= iTmaxStart && iT % iTperiod == 0 ) { + SinusStartScale startScale( iTmaxStart, T( 1 ) ); + int iTvec[1] = { iT }; + T frac[1] = { T( 0 ) }; + startScale( frac, iTvec ); + T maxVelocity = frac[0] * converter.getCharLatticeVelocity() * 3. / 4. + * std::pow( inletRadius, 2 ) / std::pow( outletRadius0, 2 ); + + CirclePoiseuille3D poiseuilleU4( outletCenter0[0], outletCenter0[1], + outletCenter0[2], outletNormal0[0], + outletNormal0[1], outletNormal0[2], + outletRadius0 * 0.95, -maxVelocity ); + + CirclePoiseuille3D poiseuilleU5( outletCenter1[0], outletCenter1[1], + outletCenter1[2], outletNormal1[0], + outletNormal1[1], outletNormal1[2], + outletRadius1 * 0.95, -maxVelocity ); + + sLattice.defineU( superGeometry, 4, poiseuilleU4 ); + sLattice.defineU( superGeometry, 5, poiseuilleU5 ); + } +} + +// Computes the pressure drop between voxels before and after the cylinder +bool getResults( SuperLattice3D& sLattice, + UnitConverter const& converter, int iT, int iTperiod, + SuperGeometry3D& superGeometry, + Timer& fluidTimer, STLreader& stlReader, + SuperParticleSystem3D& supParticleSystem, + T radii, T partRho, Timer& particleTimer, + SuperParticleSysVtuWriter& supParticleWriter, + bool fluidExists) +{ + + OstreamManager clout( std::cout, "getResults" ); + SuperVTMwriter3D vtmWriter( "bifurcation3d" ); + SuperVTMwriter3D vtmWriterStartTime( "startingTimeBifurcation3d" ); + + SuperLatticePhysVelocity3D velocity( sLattice, converter ); + SuperLatticePhysPressure3D pressure( sLattice, converter ); + vtmWriter.addFunctor( velocity ); + vtmWriter.addFunctor( pressure ); + + vtmWriterStartTime.addFunctor( velocity ); + vtmWriterStartTime.addFunctor( pressure ); + + int fluidMaxT = converter.getLatticeTime( fluidMaxPhysT ); + + if ( iT == 0 ) { + SuperLatticeGeometry3D geometry( sLattice, superGeometry ); + SuperLatticeCuboid3D cuboid( sLattice ); + SuperLatticeRank3D rank( sLattice ); + vtmWriter.write( geometry ); + vtmWriter.write( cuboid ); + vtmWriter.write( rank ); + vtmWriter.createMasterFile(); + vtmWriterStartTime.createMasterFile(); + + + // Print some output of the chosen simulation setup + clout << "N=" << N <<"; maxTimeSteps(fluid)=" + << converter.getLatticeTime( fluidMaxPhysT ) << "; noOfCuboid=" + << superGeometry.getCuboidGeometry().getNc() << "; Re=" << Re + << "; noOfParticles=" << noOfParticles << "; maxTimeSteps(particle)=" + << converter.getLatticeTime( particleMaxPhysT ) + << "; St=" << ( 2.*partRho*radius*radius*converter.getCharPhysVelocity() ) / ( 9.*converter.getPhysViscosity()*converter.getPhysDensity()*converter.getCharPhysLength() ) << std::endl; + } + + // Writes the vtk and gif files + if ( iT % iTperiod == 0 ) { + if ( !fluidExists && iT <= fluidMaxT ) { + vtmWriterStartTime.write(iT); + SuperEuklidNorm3D normVel( velocity ); + BlockReduction3D2D planeReduction( normVel, {0, -1, 0}, 600, BlockDataSyncMode::ReduceOnly ); + // write output as JPEG + heatmap::write(planeReduction, iT); + } + if (iT > fluidMaxT) { + // only write vtk-files after the fluid calculation is finished + vtmWriter.write(iT - fluidMaxT); + } + } + + // Writes output on the console for the fluid phase + if (iT < converter.getLatticeTime( fluidMaxPhysT ) && iT%iTperiod == 0 ) { + + // Timer statics + fluidTimer.update( iT ); + fluidTimer.printStep(); + + // Lattice statistics + sLattice.getStatistics().print( iT, converter.getPhysTime( iT ) ); + + // Flux at the inlet and outlet regions + const std::vector materials = { 1, 3, 4, 5 }; + + IndicatorCircle3D inlet( + inletCenter + 2. * converter.getConversionFactorLength() * inletNormal, + inletNormal, inletRadius + 2. * converter.getConversionFactorLength() ); + SuperPlaneIntegralFluxVelocity3D vFluxInflow( sLattice, converter, superGeometry, inlet, materials ); + vFluxInflow.print( "inflow", "ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxInflow( sLattice, converter, superGeometry, inlet, materials ); + pFluxInflow.print( "inflow", "N", "Pa" ); + + IndicatorCircle3D outlet0( + outletCenter0 + 2. * converter.getConversionFactorLength() * outletNormal0, + outletNormal0, outletRadius0 + 2. * converter.getConversionFactorLength() ); + SuperPlaneIntegralFluxVelocity3D vFluxOutflow0( sLattice, converter, superGeometry, outlet0, materials ); + vFluxOutflow0.print( "outflow0", "ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxOutflow0( sLattice, converter, superGeometry, outlet0, materials ); + pFluxOutflow0.print( "outflow0", "N", "Pa" ); + + IndicatorCircle3D outlet1( + outletCenter1 + 2. * converter.getConversionFactorLength() * outletNormal1, + outletNormal1, outletRadius1 + 2. * converter.getConversionFactorLength() ); + SuperPlaneIntegralFluxVelocity3D vFluxOutflow1( sLattice, converter, superGeometry, outlet1, materials ); + vFluxOutflow1.print( "outflow1", "ml/s" ); + SuperPlaneIntegralFluxPressure3D pFluxOutflow1( sLattice, converter, superGeometry, outlet1, materials ); + pFluxOutflow1.print( "outflow1", "N", "Pa" ); + } + + // Writes output on the console for the fluid phase + if ( iT >= converter.getLatticeTime( fluidMaxPhysT ) && + (iT%iTperiod == 0 || iT == converter.getLatticeTime( fluidMaxPhysT )) ) { + + particleTimer.print( iT - fluidMaxT ); + + // console output number of particles at different material numbers mat + supParticleSystem.print( {1,2,3,4,5} ); + // console output of escape (E), capture (C) rate for material numbers mat + supParticleSystem.captureEscapeRate( {4,5} ); + + // only write vtk-files after the fluid calculation is finished + supParticleWriter.write( iT - fluidMaxT ); + + // true as long as certain amount of active particles + if ( supParticleSystem.globalNumOfActiveParticles() < 0.001 * noOfParticles + && iT > 0.9*converter.getLatticeTime( fluidMaxPhysT + particleMaxPhysT ) ) { + return false; + } + } + return true; +} + +int main( int argc, char* argv[] ) +{ + + // === 1st Step: Initialization === + + olbInit( &argc, &argv ); + + singleton::directories().setOutputDir( "./tmp/" ); + OstreamManager clout( std::cout, "main" ); + + UnitConverterFromResolutionAndRelaxationTime const converter( + int {N}, // resolution: number of voxels per charPhysL + (T) 0.557646, // latticeRelaxationTime: relaxation time, have to be greater than 0.5! + (T) inletRadius*2., // charPhysLength: reference length of simulation geometry + (T) Re*1.5e-5/( inletRadius*2 ), // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__ + (T) 1.5e-5, // physViscosity: physical kinematic viscosity in __m^2 / s__ + (T) 1.225 // physDensity: physical density in __kg / m^3__ + ); + // Prints the converter log as console output + converter.print(); + // Writes the converter log in a file + converter.write("bifurcation3d"); + + // === 2nd Step: Prepare Geometry === + STLreader stlReader( "../bifurcation3d.stl", converter.getConversionFactorLength() ); + IndicatorLayer3D extendedDomain( stlReader, + converter.getConversionFactorLength() ); + + // Instantiation of an empty cuboidGeometry + int noOfCuboids = std::max( 16, 4 * singleton::mpi().getSize() ); + + CuboidGeometry3D cuboidGeometry( extendedDomain, converter.getConversionFactorLength(), + noOfCuboids ); + + // Instantiation of an empty loadBalancer + HeuristicLoadBalancer loadBalancer( cuboidGeometry ); + // Default instantiation of superGeometry + SuperGeometry3D superGeometry( cuboidGeometry, loadBalancer, 2 ); + + prepareGeometry( converter, extendedDomain, stlReader, superGeometry ); + + // === 3rd Step: Prepare Lattice === + + SuperLattice3D sLattice( superGeometry ); + + BGKdynamics bulkDynamics( converter.getLatticeRelaxationFrequency(), + instances::getBulkMomenta() ); + + sOnLatticeBoundaryCondition3D sBoundaryCondition( sLattice ); + createInterpBoundaryCondition3D( sBoundaryCondition ); + + sOffLatticeBoundaryCondition3D sOffBoundaryCondition( + sLattice ); + createBouzidiBoundaryCondition3D( sOffBoundaryCondition ); + + prepareLattice( sLattice, converter, bulkDynamics, sBoundaryCondition, + sOffBoundaryCondition, stlReader, superGeometry ); + + // === 3.1 Step: Particles === + clout << "Prepare Particles ..." << std::endl; + + // SuperParticleSystems3D + SuperParticleSystem3D supParticleSystem( superGeometry ); + // define which properties are to be written in output data + SuperParticleSysVtuWriter supParticleWriter( supParticleSystem, + "particles", SuperParticleSysVtuWriter::particleProperties:: + velocity + | SuperParticleSysVtuWriter::particleProperties::mass + | SuperParticleSysVtuWriter::particleProperties::radius + | SuperParticleSysVtuWriter::particleProperties::active ); + + SuperLatticeInterpPhysVelocity3D getVel( sLattice, converter ); + + auto stokesDragForce = make_shared + < StokesDragForce3D + > ( getVel, converter ); + + // material numbers where particles should be reflected + std::set boundMaterial = { 2, 4, 5}; + auto materialBoundary = make_shared + < MaterialBoundary3D + > ( superGeometry, boundMaterial ); + + supParticleSystem.addForce( stokesDragForce ); + supParticleSystem.addBoundary( materialBoundary ); + supParticleSystem.setOverlap( 2. * converter.getConversionFactorLength() ); + + // particles generation at inlet3 + Vector c( inletCenter ); + c[2] = 0.074; + IndicatorCircle3D inflowCircle( c, inletNormal, inletRadius - + converter.getConversionFactorLength() * 2.5 ); + IndicatorCylinder3D inletCylinder( inflowCircle, 0.01 * + converter.getConversionFactorLength() ); + supParticleSystem.addParticle( inletCylinder, 4. / 3. * M_PI * + std::pow( radius, 3 ) * partRho, radius, noOfParticles ); + + clout << "Prepare Particles ... OK" << std::endl; + + // === 4th Step: Main Loop with Timer === + + Timer fluidTimer( converter.getLatticeTime( fluidMaxPhysT ), + superGeometry.getStatistics().getNvoxel() ); + + Timer particleTimer( converter.getLatticeTime( particleMaxPhysT ), + noOfParticles ); + fluidTimer.start(); + + int iT = 0; + // amount of timesteps when getResults rewrites data + int iTperiod = converter.getLatticeTime( .2 ); + + bool fluidExists = true; + + // checks whether there is already data of the fluid from an earlier calculation + if ( !( sLattice.load( "fluidSolution" ) ) ) { + + fluidExists = false; + + // if there is no data available, it is generated + for ( ; iT <= converter.getLatticeTime( fluidMaxPhysT ); ++iT ) { + + // during run up time boundary values are set, collide and stream step, + // results of fluid, afterwards only particles are simulated + setBoundaryValues( sLattice, converter, iT, fluidMaxPhysT, superGeometry ); + sLattice.collideAndStream(); + + getResults( sLattice, converter, iT, iTperiod, superGeometry, fluidTimer, stlReader, + supParticleSystem, radius, partRho, particleTimer, + supParticleWriter, fluidExists ); + } + + fluidTimer.stop(); + fluidTimer.printSummary(); + + sLattice.communicate(); + // calculated results are written in a file + sLattice.save( "fluidSolution" ); + } + + // if there exists already data of the fluid from an earlier calculation, this is used + else { + + iT = converter.getLatticeTime( fluidMaxPhysT ); + getResults( sLattice, converter, iT, + iTperiod, superGeometry, fluidTimer, stlReader, + supParticleSystem, radius, partRho, particleTimer, + supParticleWriter, fluidExists ); + + } + + // after the fluid calculation, particle simulation starts + supParticleSystem.setVelToFluidVel( getVel ); + particleTimer.start(); + + for ( ; iT <= converter.getLatticeTime( fluidMaxPhysT + particleMaxPhysT ); ++iT ) { + // particles simulation starts after run up time is over + supParticleSystem.simulate( converter.getConversionFactorTime() ); + + if ( !getResults( sLattice, converter, iT, + iTperiod, superGeometry, fluidTimer, + stlReader, supParticleSystem, radius, partRho, + particleTimer, supParticleWriter, fluidExists) ) { + break; + } + } + particleTimer.stop(); + particleTimer.printSummary(); +} diff --git a/examples/particles/bifurcation3d/eulerLagrange/definitions.mk b/examples/particles/bifurcation3d/eulerLagrange/definitions.mk new file mode 100644 index 0000000..3ad56e8 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerLagrange/definitions.mk @@ -0,0 +1,30 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + + +########################################################################### +########################################################################### +## DEFINITIONS TO BE CHANGED + +ROOT := ../../../.. +SRC := bifurcation3d.cpp +OUTPUT := bifurcation3d diff --git a/examples/particles/bifurcation3d/eulerLagrange/module.mk b/examples/particles/bifurcation3d/eulerLagrange/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/particles/bifurcation3d/eulerLagrange/module.mk @@ -0,0 +1,29 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2017 Markus Mohrhard +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + +current_dir := $(dir $(word $(words $(MAKEFILE_LIST)),$(MAKEFILE_LIST))) + +include global.mk +include rules.mk +include $(addsuffix definitions.mk, $(current_dir)) + +$(eval $(call sample,$(current_dir)$(OUTPUT),$(addprefix $(current_dir), $(SRC)))) diff --git a/examples/particles/dkt2d/Makefile b/examples/particles/dkt2d/Makefile new file mode 100644 index 0000000..910c343 --- /dev/null +++ b/examples/particles/dkt2d/Makefile @@ -0,0 +1,115 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + + +########################################################################### +########################################################################### +## DEFINITIONS TO BE CHANGED + +ROOT := ../../.. +SRC := +OUTPUT := dkt2d + +########################################################################### +## definitions + +include $(ROOT)/global.mk + +OBJECTS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).o) +DEPS := $(foreach file, $(SRC) $(OUTPUT), $(PWD)/$(file).d) + +########################################################################### +## all + +all : depend compile updatelib link + + +########################################################################### +## dependencies + +depend : $(DEPS) + +$(PWD)/%.d : %.cpp + @echo Create dependencies for $< + @$(SHELL) -ec '$(CXX) -M $(CXXFLAGS) $(IDIR) $< \ + | sed -e "s!$*\.o!$(PWD)\/$*\.o!1" > .tmpfile; \ + cp -f .tmpfile $@;' + +########################################################################### +## compile + +compile : $(OBJECTS) + +$(PWD)/%.o: %.cpp + @echo Compile $< + $(CXX) $(CXXFLAGS) $(IDIR) -c $< -o $@ + +########################################################################### +## clean + +clean : cleanrub cleanobj cleandep + +cleanrub: + @echo Clean rubbish files + @rm -f *~ core .tmpfile tmp/*.* $(OUTPUT) + @rm -f *.pdf *.dat + @rm -f tmp/vtkData/*.* tmp/vtkData/data/*.* tmp/imageData/*.* + @rm -f tmp/gnuplotData/*.* tmp/gnuplotData/data/*.* + +cleanobj: + @echo Clean object files + @rm -f $(OBJECTS) + +cleandep: + @echo Clean dependencies files + @rm -f $(DEPS) + +cleanbuild: + @echo Clean olb main + @cd $(ROOT); \ + $(MAKE) cleanbuild; + +########################################################################### +## update lib + +updatelib : + @cd $(ROOT); \ + $(MAKE) all; + +########################################################################### +## link + +link: $(OUTPUT) + +$(OUTPUT): $(OBJECTS) $(ROOT)/$(LIBDIR)/lib$(LIB).a + @echo Link $@ + $(CXX) $(foreach file, $(SRC), $(file).o) $@.o $(LDFLAGS) -L$(ROOT)/$(LIBDIR) -l$(LIB) -lz -o $@ + +########################################################################### +## include dependencies + +ifneq "$(strip $(wildcard *.d))" "" + include $(foreach file,$(DEPS),$(file)) +endif + +########################################################################### +########################################################################### diff --git a/examples/particles/dkt2d/definitions.mk b/examples/particles/dkt2d/definitions.mk new file mode 100644 index 0000000..e68e06b --- /dev/null +++ b/examples/particles/dkt2d/definitions.mk @@ -0,0 +1,30 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# E-mail contact: info@openlb.net +# The most recent release of OpenLB can be downloaded at +# +# +# This program 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 2 +# of the License, or (at your option) any later version. +# +# This program 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 this program; if not, write to the Free +# Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +# Boston, MA 02110-1301, USA. + + +########################################################################### +########################################################################### +## DEFINITIONS TO BE CHANGED + +ROOT := ../../.. +SRC := dkt2d.cpp +OUTPUT := dkt2d diff --git a/examples/particles/dkt2d/dkt.p b/examples/particles/dkt2d/dkt.p new file mode 100644 index 0000000..a7ebb8f --- /dev/null +++ b/examples/particles/dkt2d/dkt.p @@ -0,0 +1,22 @@ +reset +set grid +set terminal pdf enhanced +set xrange [0:6] +set xtics 0,1,6 +set xlabel "Time [s]" + +set output "x_coordinates.pdf" +set ylabel "horizontal position [m]" +set yrange [0:0.02] +set ytics 0,0.002,0.02 +set key left bottom +plot "gnuplot.dat" using 1:4 with lines title "Leading Particle", \ + "gnuplot.dat" using 1:5 with lines title "Trailing Particle" + +set output "y_coordinates.pdf" +set ylabel "vertical position [m]" +set yrange [0:0.08] +set ytics 0,0.005,0.08 +set key left bottom +plot "gnuplot.dat" using 1:2 with lines title "Leading Particle", \ + "gnuplot.dat" using 1:3 with lines title "Trailing Particle" diff --git a/examples/particles/dkt2d/dkt2d.cpp b/examples/particles/dkt2d/dkt2d.cpp new file mode 100644 index 0000000..e90f49c --- /dev/null +++ b/examples/particles/dkt2d/dkt2d.cpp @@ -0,0 +1,267 @@ +/* dkt2d.cpp: + * The case examines the settling of two circles under gravity + * in a surrounding fluid. The rectangular domain is limited + * by no-slip boundary conditions. + * For the calculation of forces a DNS approach is chosen + * which also leads to a back-coupling of the particle on the fluid, + * inducing a flow. + * The simulation is based on the homogenised lattice Boltzmann approach + * (HLBM) introduced by Krause et al. in "Particle flow simulations + * with homogenised lattice Boltzmann methods". + * The drafting-kissing-tumbling benchmark case is e.g. described + * in "Drafting, kissing and tumbling process of two particles + * with different sizes" by Wang et al. +