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/settlingCube3d/Makefile | 105 ++++++++ examples/particles/settlingCube3d/definitions.mk | 30 +++ examples/particles/settlingCube3d/module.mk | 29 ++ .../particles/settlingCube3d/settlingCube3D.cpp | 294 +++++++++++++++++++++ 4 files changed, 458 insertions(+) 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/settlingCube3d') diff --git a/examples/particles/settlingCube3d/Makefile b/examples/particles/settlingCube3d/Makefile new file mode 100755 index 0000000..a953954 --- /dev/null +++ b/examples/particles/settlingCube3d/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/settlingCube3d/definitions.mk b/examples/particles/settlingCube3d/definitions.mk new file mode 100755 index 0000000..0a286c7 --- /dev/null +++ b/examples/particles/settlingCube3d/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 := settlingCube3D.cpp +OUTPUT := settlingCube3D diff --git a/examples/particles/settlingCube3d/module.mk b/examples/particles/settlingCube3d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/particles/settlingCube3d/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/settlingCube3d/settlingCube3D.cpp b/examples/particles/settlingCube3d/settlingCube3D.cpp new file mode 100644 index 0000000..f075927 --- /dev/null +++ b/examples/particles/settlingCube3d/settlingCube3D.cpp @@ -0,0 +1,294 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2006-2015 Fabian Klemens, Jonas Latt, Mathias J. Krause + * Vojtech Cvrcek, Peter Weisbrod + * 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. + */ + +/* settlingCube3d.cpp: + * The case examines the settling of a cubical silica particle + * under the influence of gravity. + * The object is surrounded by water in a rectangular domain + * limited by no-slip boundary conditions. + * For the calculation of forces an 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 in "Particle flow simulations with homogenised + * lattice Boltzmann methods" by Krause et al. + * and extended in "Towards the simulation of arbitrarily shaped 3D particles + * using a homogenised lattice Boltzmann method" by Trunk et al. + * for the simulation of 3D particles. + * + * This example demonstrates the usage of HLBM in the OpenLB framework. + */ + +#include "olb3D.h" +#include "olb3D.hh" // use generic version only! + +#include +#include +#include +#include + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; +using namespace std; + +typedef double T; +#define DESCRIPTOR D3Q19 + +#define WriteVTK + +// Discretization Settings +int res = 30; +T const charLatticeVelocity = 0.01; + +// Time Settings +T const maxPhysT = 0.5; // max. simulation time in s +T const iTwrite = 0.02; // write out intervall in s + +// Domain Settings +T const lengthX = 0.01; +T const lengthY = 0.01; +T const lengthZ = 0.05; + +// Fluid Settings +T const physDensity = 1000; +T const physViscosity = 1E-5; + +//Particle Settings +T centerX = lengthX*.5; +T centerY = lengthY*.5; +T centerZ = lengthZ*.9; +T const cubeDensity = 2500; +T const cubeEdgeLength = 0.0025; +Vector cubeCenter = {centerX,centerY,centerZ}; +Vector cubeOrientation = {0.,15.,0.}; +Vector cubeVelocity = {0.,0.,0.}; +Vector externalAcceleration = {.0, .0, -9.81 * (1. - physDensity / cubeDensity)}; + +// Characteristic Quantities +T const charPhysLength = lengthX; +T const charPhysVelocity = 0.15; // Assumed maximal velocity + + +// Prepare geometry +void prepareGeometry(UnitConverter const& converter, + SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout, "prepareGeometry"); + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename(0, 2); + superGeometry.rename(2, 1, 1, 1, 1); + + superGeometry.clean(); + superGeometry.innerClean(); + + superGeometry.checkForErrors(); + superGeometry.getStatistics().print(); + clout << "Prepare Geometry ... OK" << std::endl; + return; +} + + +// Set up the geometry of the simulation +void prepareLattice( + SuperLattice3D& sLattice, UnitConverter const& converter, + Dynamics& designDynamics, + sOnLatticeBoundaryCondition3D& sBoundaryCondition, + SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout, "prepareLattice"); + clout << "Prepare Lattice ..." << std::endl; + clout << "setting Velocity Boundaries ..." << std::endl; + + /// Material=0 -->do nothing + sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics()); + sLattice.defineDynamics(superGeometry, 1, &designDynamics); + sLattice.defineDynamics(superGeometry, 2, &instances::getBounceBack()); + + clout << "Prepare Lattice ... OK" << std::endl; +} + + +//Set Boundary Values +void setBoundaryValues(SuperLattice3D& sLattice, + UnitConverter const& converter, int iT, + SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout, "setBoundaryValues"); + + if (iT == 0) { + AnalyticalConst3D zero(0.); + AnalyticalConst3D one(1.); + sLattice.defineField(superGeometry.getMaterialIndicator({0,1,2}), one); + // Set initial condition + AnalyticalConst3D ux(0.); + AnalyticalConst3D uy(0.); + AnalyticalConst3D uz(0.); + AnalyticalConst3D rho(1.); + AnalyticalComposed3D u(ux, uy, uz); + + //Initialize all values of distribution functions to their local equilibrium + sLattice.defineRhoU(superGeometry, 1, rho, u); + sLattice.iniEquilibrium(superGeometry, 1, rho, u); + + // Make the lattice ready for simulation + sLattice.initialize(); + } +} + + +/// Computes the pressure drop between the voxels before and after the cylinder +void getResults(SuperLattice3D& sLattice, + UnitConverter const& converter, int iT, + SuperGeometry3D& superGeometry, Timer& timer, + ParticleDynamics3D particleDynamics) +{ + OstreamManager clout(std::cout, "getResults"); + +#ifdef WriteVTK + SuperVTMwriter3D vtkWriter("sedimentation"); + SuperLatticePhysVelocity3D velocity(sLattice, converter); + SuperLatticePhysPressure3D pressure(sLattice, converter); + SuperLatticePhysExternalPorosity3D externalPor(sLattice, converter); + vtkWriter.addFunctor(velocity); + vtkWriter.addFunctor(pressure); + vtkWriter.addFunctor(externalPor); + + if (iT == 0) { + /// Writes the converter log file + SuperLatticeGeometry3D geometry(sLattice, superGeometry); + SuperLatticeCuboid3D cuboid(sLattice); + SuperLatticeRank3D rank(sLattice); + vtkWriter.write(geometry); + vtkWriter.write(cuboid); + vtkWriter.write(rank); + vtkWriter.createMasterFile(); + } + + if (iT % converter.getLatticeTime(iTwrite) == 0) { + vtkWriter.write(iT); + } +#endif + + /// Writes output on the console + if (iT % converter.getLatticeTime(iTwrite) == 0) { + timer.update(iT); + timer.printStep(); + sLattice.getStatistics().print(iT, converter.getPhysTime(iT)); + particleDynamics.print(); + } +} + +int main(int argc, char* argv[]) +{ + /// === 1st Step: Initialization === + olbInit(&argc, &argv); + singleton::directories().setOutputDir("./tmp/"); + OstreamManager clout(std::cout, "main"); + + UnitConverterFromResolutionAndLatticeVelocity converter( + (int) res, //resolution + ( T ) charLatticeVelocity, //charLatticeVelocity + ( T ) charPhysLength, //charPhysLength + ( T ) charPhysVelocity, //charPhysVelocity + ( T ) physViscosity, //physViscosity + ( T ) physDensity //physDensity + ); + converter.print(); + + /// === 2rd Step: Prepare Geometry === + /// Instantiation of a cuboidGeometry with weights + std::vector extend(3, T()); + extend[0] = lengthX; + extend[1] = lengthY; + extend[2] = lengthZ; + std::vector origin(3, T()); + IndicatorCuboid3D cuboid(extend, origin); + +#ifdef PARALLEL_MODE_MPI + CuboidGeometry3D cuboidGeometry(cuboid, converter.getConversionFactorLength(), singleton::mpi().getSize()); +#else + CuboidGeometry3D cuboidGeometry(cuboid, converter.getConversionFactorLength(), 7); +#endif + cuboidGeometry.print(); + + HeuristicLoadBalancer loadBalancer(cuboidGeometry); + SuperGeometry3D superGeometry(cuboidGeometry, loadBalancer, 2); + prepareGeometry(converter, superGeometry); + + /// === 3rd Step: Prepare Lattice === + SuperLattice3D sLattice(superGeometry); + + PorousParticleBGKdynamics designDynamics(converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta()); + + sOnLatticeBoundaryCondition3D sBoundaryCondition(sLattice); + createLocalBoundaryCondition3D(sBoundaryCondition); + + prepareLattice(sLattice, converter, designDynamics, sBoundaryCondition, superGeometry); + + /// === 4th Step: Main Loop with Timer === + Timer timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel()); + timer.start(); + + // Create Particle Dynamics + ParticleDynamics3D particleDynamics(sLattice, converter, superGeometry, lengthX, lengthY, lengthZ, externalAcceleration); + + // Create Cube Indicator + T epsilon = 0.5*converter.getConversionFactorLength(); + + //Cube indicator + SmoothIndicatorCuboid3D particleIndicator(cubeCenter, cubeEdgeLength, cubeEdgeLength, cubeEdgeLength, epsilon, cubeOrientation, cubeDensity, cubeVelocity); + + //Sphere indicator + //SmoothIndicatorSphere3D particleIndicator(cubeCenter, 0.5*cubeEdgeLength, epsilon, cubeDensity, cubeVelocity); + + //Cylinder indicator + //SmoothIndicatorCylinder3D particleIndicator(cubeCenter, { 1, 0, 0 }, 0.5*cubeEdgeLength, cubeEdgeLength, epsilon, cubeOrientation, cubeDensity, cubeVelocity); + + SuperExternal3D superExtPorosity(superGeometry, sLattice, sLattice.getOverlap()); + SuperExternal3D superExtNumerator(superGeometry, sLattice, sLattice.getOverlap()); + SuperExternal3D superExtDenominator(superGeometry, sLattice, sLattice.getOverlap()); + particleDynamics.addParticle( particleIndicator ); + particleDynamics.print(); + + /// === 5th Step: Definition of Initial and Boundary Conditions === + setBoundaryValues(sLattice, converter, 0, superGeometry); + + clout << "MaxIT: " << converter.getLatticeTime(maxPhysT) << std::endl; + for (int iT = 0; iT < converter.getLatticeTime(maxPhysT)+10; ++iT) { + particleDynamics.simulateTimestep("verlet"); + getResults(sLattice, converter, iT, superGeometry, timer, particleDynamics); + sLattice.collideAndStream(); + superExtPorosity.communicate(); + superExtNumerator.communicate(); + superExtDenominator.communicate(); + } + + timer.stop(); + timer.printSummary(); +} -- cgit v1.2.3