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/magneticParticles3d/Makefile | 105 ++++ .../particles/magneticParticles3d/definitions.mk | 30 + .../magneticParticles3d/magneticParticles3d.cpp | 644 +++++++++++++++++++++ .../magneticParticles3d/magneticParticles3d.xml | 47 ++ examples/particles/magneticParticles3d/module.mk | 29 + 5 files changed, 855 insertions(+) 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 (limited to 'examples/particles/magneticParticles3d') diff --git a/examples/particles/magneticParticles3d/Makefile b/examples/particles/magneticParticles3d/Makefile new file mode 100644 index 0000000..a953954 --- /dev/null +++ b/examples/particles/magneticParticles3d/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/magneticParticles3d/definitions.mk b/examples/particles/magneticParticles3d/definitions.mk new file mode 100644 index 0000000..02bf7e1 --- /dev/null +++ b/examples/particles/magneticParticles3d/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 := magneticParticles3d.cpp +OUTPUT := magneticParticles3d diff --git a/examples/particles/magneticParticles3d/magneticParticles3d.cpp b/examples/particles/magneticParticles3d/magneticParticles3d.cpp new file mode 100644 index 0000000..673733d --- /dev/null +++ b/examples/particles/magneticParticles3d/magneticParticles3d.cpp @@ -0,0 +1,644 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2018 Sascha Janz, 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. + */ + +/* magneticParticles3d.cpp: + * High-gradient magnetic separation is a method + * to separate ferromagnetic particles from a suspension. + * The simulation shows the deposition of magnetic particles + * on a single magnetized wire and models + * the magnetic separation step of the complete process. + */ + + +#include "olb3D.h" +#ifndef OLB_PRECOMPILED // Unless precompiled version is used +#include "olb3D.hh" // Include full template code +#endif + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; +using namespace std; + +#define DESCRIPTOR D3Q19<> + +#define PARTICLE MagneticParticle3D + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +typedef double T; + +int iT = 0; + +// simulation geometry dimensions +T geometryLengthX = 3.e-3 ; // length x-dircetion [m] +T geometryLengthY = 2.e-3 ; // length y-dircetion [m] +T geometryLengthZ = 1.e-3 ; // length z-dircetion [m] + +// magnetic particle paramters +T pRadius = 1.e-5 ; // particle radius [m] +T pDensity = 3250.; // particle density [kg / m^3] +T pVolume = 4. / 3. * M_PI * std::pow(pRadius, 3.); // particle volume [m^3] +T pMass = pVolume * pDensity; // particle mass [kg] +T pMag = 67. * pDensity * 0.06; // particle sat. magnetization [A / m] +T elastModulus = 2.e4; // elastic modulus [Pa] +T poissonRatio = 0.3; // poisson's ratio [-] +T shearModulus = elastModulus / (2. * (1. + poissonRatio)); // shear modulus [Pa] + +// magnetic wire paramters +T wRadius = 5.e-4; // wire radius [m] +T wLength = geometryLengthZ; // wire length [m] +T wDensity = 7874.; // wire (iron) density [kg / m^3] +T wMag = 215. * wDensity * 8.5e-3; // wire sat. magnetization [A / m] +// mass-specific sat. magnetization iron = 215. [A m^2 / kg] + +// Rayleigh time step +// can be used as an orientation value for particle time step sizes +T eta = 0.8766 + 0.1631 * poissonRatio; // [-] +T t_Rayleigh = (M_PI * pRadius / eta) * std::sqrt(pDensity / shearModulus); // [s] + +void prepareGeometry(UnitConverter const& converter, IndicatorF3D& wire, + IndicatorF3D& indicator, IndicatorF3D& extendedDomain, + SuperGeometry3D& superGeometry) +{ + + OstreamManager clout(std::cout, "prepareGeometry"); + clout << "Prepare Geometry ..." << std::endl; + + superGeometry.rename(0, 2, extendedDomain); + superGeometry.rename(2, 1, indicator); + superGeometry.rename(1, 5, wire); + + T minR = superGeometry.getStatistics().getMinPhysR(2)[0]; + minR -= 0.5 * converter.getConversionFactorLength(); + std::vector < T > center = superGeometry.getStatistics().getCenterPhysR(2); + std::vector physExtend = superGeometry.getStatistics().getPhysExtend(1); + + Vector origin = { minR, center[1], center[2] }; + Vector extend = { converter.getConversionFactorLength(), physExtend[1], physExtend[2] }; + + IndicatorCuboid3D inlet(extend[0], extend[1], extend[2], origin); + origin[0] += superGeometry.getStatistics().getPhysExtend(2)[0]; + IndicatorCuboid3D outlet(extend[0], extend[1], extend[2], origin); + + // rename the material at the inlet + superGeometry.rename(2, 3, inlet); + // rename the material at the outlet + superGeometry.rename(2, 4, outlet); + + superGeometry.clean(); + superGeometry.innerClean(); + superGeometry.checkForErrors(); + superGeometry.print(); + + clout << "Prepare Geometry ... OK" << std::endl; +} + +/// Set up the geometry of the simulation +void prepareLattice(SuperLattice3D& sLattice, + UnitConverter const& converter, IndicatorF3D& wire, + Dynamics& bulkDynamics, + sOnLatticeBoundaryCondition3D& sOnBC, + sOffLatticeBoundaryCondition3D& offBc, + SuperGeometry3D& superGeometry) +{ + OstreamManager clout(std::cout, "prepareLattice"); + clout << "Prepare Lattice ..." << std::endl; + + const T omega = (1. / converter.getLatticeRelaxationTime()); + + /// 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()); + + sLattice.defineDynamics(superGeometry, 3, &bulkDynamics); + sLattice.defineDynamics(superGeometry, 4, &bulkDynamics); + + /// Material=5 -->do nothing + sLattice.defineDynamics(superGeometry, 5, + &instances::getNoDynamics()); + offBc.addZeroVelocityBoundary(superGeometry, 5, wire); + + // boundary conditions for fluid + + // inlet + sOnBC.addVelocityBoundary(superGeometry, 3, omega); + // outlet + sOnBC.addPressureBoundary(superGeometry, 4, omega); + + // initialisation + AnalyticalConst3D roh(1.); + AnalyticalConst3D u0(0., 0., 0.); + + sLattice.defineRhoU(superGeometry, 1, roh, u0); + sLattice.iniEquilibrium(superGeometry, 1, roh, u0); + sLattice.defineRhoU(superGeometry, 3, roh, u0); + sLattice.iniEquilibrium(superGeometry, 3, roh, u0); + sLattice.defineRhoU(superGeometry, 4, roh, u0); + sLattice.iniEquilibrium(superGeometry, 4, roh, u0); + + sLattice.initialize(); + + clout << "Prepare Lattice ... OK" << std::endl; +} + +void setBoundaryValues(SuperLattice3D& sLattice, + UnitConverter const& converter, SuperParticleSystem3D& spSys, + int iT, int itStartScaleT, SuperGeometry3D& superGeometry, bool outNS, bool outLP) +{ + + OstreamManager clout(std::cout, "setBoundaryValues"); + std::vector < T > maxVelocity(3, T()); + T distanceToBoundary = converter.getConversionFactorLength() / 2.; + + if (outNS && iT <= itStartScaleT) { + + SinusStartScale startScale(itStartScaleT, T(1)); + int help[1] = { iT }; + T frac[3] = { T() }; + startScale(frac, help); + + maxVelocity[0] = converter.getCharPhysVelocity() * frac[0]; + + // outlet + RectanglePoiseuille3D u3(superGeometry, 3, maxVelocity, + distanceToBoundary, distanceToBoundary, distanceToBoundary); + + sLattice.defineU(superGeometry, 3, u3); + } + + // to reset and adopt boundaries in case of loaded fluid data + if (outLP && iT == 0) { + + maxVelocity[0] = converter.getCharPhysVelocity(); + + // inlet + RectanglePoiseuille3D u3(superGeometry, 3, maxVelocity, + distanceToBoundary, distanceToBoundary, distanceToBoundary); + + clout << "BC for loaded fluid data is reset on mat 3 " << std::endl; + + sLattice.defineU(superGeometry, 3, u3); + } +} + +void getResults(SuperGeometry3D& superGeometry, + SuperLattice3D& sLattice, + UnitConverter const& converter, + AnalyticalF3D& magForce, AnalyticalF3D& magField, + SuperParticleSystem3D& spSys, + SuperParticleSysVtuWriterMag& particleOut, Timer& timer, int& iT, + int& itConsoleOutputFluid, int& itVtkOutputFluid, + int& itConsoleOutputMagParticles, int& itVtkOutputMagParticles, + bool& outNS, bool& outLP) +{ + + OstreamManager clout(std::cout, "getResults"); + + // start vtm objects data until fluid reaches equilibrium + SuperVTMwriter3D vtmWriterFluidStart("fluidReachingEquilibrium"); + + SuperLatticePhysVelocity3D velocity(sLattice, converter); + SuperLatticeDensity3D density(sLattice); + SuperLatticeFfromAnalyticalF3D superMagPForceOne(magForce, + sLattice); + SuperLatticeFfromAnalyticalF3D superMagPFieldOne (magField, sLattice); + + vtmWriterFluidStart.addFunctor(velocity); + vtmWriterFluidStart.addFunctor(density); + vtmWriterFluidStart.addFunctor(superMagPForceOne); + + if (outNS && iT == 0) { + SuperLatticeGeometry3D geometry(sLattice, superGeometry); + SuperLatticeCuboid3D cuboid(sLattice); + SuperLatticeRank3D rank(sLattice); + + vtmWriterFluidStart.write(geometry); + vtmWriterFluidStart.write(cuboid); + vtmWriterFluidStart.write(rank); + + vtmWriterFluidStart.createMasterFile(); + + converter.write("magnetics 1. loop"); + + clout << "fluid computation" << std::endl; + } + if (outLP && iT == 0) { + + converter.write("magnetics 2. loop"); + clout << "magParticles computation" << std::endl; + } + + // === fluid + // console output fluid + if (outNS && iT % itConsoleOutputFluid == 0) { + timer.update(iT); + timer.print( iT ); + // output for latticeStatistics + sLattice.getStatistics().print(iT, converter.getPhysTime(iT)); + } + // vtk output fluid + if (outNS && iT % itVtkOutputFluid == 0) { + vtmWriterFluidStart.write(iT); + } + + // === particles + // Writes the console output files of the magnetic particles + if (outLP && iT % itConsoleOutputMagParticles == 0) { + /// Lattice statistics console output + timer.print( iT ); + } + // Writes the pvd files of the magnetic particles + if (outLP && iT % itVtkOutputMagParticles == 0) { + particleOut.write(iT); + } +} + +int main(int argc, char* argv[]) +{ + + /// === 1st Step: Initialization === + olbInit(&argc, &argv); + singleton::directories().setOutputDir("./tmp/"); + OstreamManager clout(std::cout, "main"); + + string fName("magneticParticles3d.xml"); + XMLreader config(fName); + + // converter contains parameters of fluid simulation + UnitConverter* converter = createUnitConverter(config); + + clout << "fluid converter: ..." << std::endl; + converter->print(); // gives overview of converter into console + + // particles time step size + T physDeltaTParticles = t_Rayleigh * 0.4; // [s] + // particles relaxation time + T tau_particles = (physDeltaTParticles * 3 * converter->getPhysViscosity() / + std::pow(converter->getConversionFactorLength(), 2.) + 0.5); + + // converter contains parameters of particle simulation + UnitConverterFromResolutionAndRelaxationTime const converterParticles( + int {converter->getResolution()}, // resolution: number of voxels per charPhysL + (T) tau_particles, // latticeRelaxationTime: relaxation time, has to be greater than 0.5! + (T) converter->getCharPhysLength(), // (hydraulic raius [m]) charPhysLength: reference length of simulation geometry + (T) converter->getCharPhysVelocity(), // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__ + (T) converter->getPhysViscosity(), // physViscosity: physical kinematic viscosity in __m^2 / s__ + (T) converter->getPhysDensity() // physDensity: physical density in __kg / m^3__ + ); + + clout << "particle converter: ..." << std::endl; + converterParticles.print(); + + // name of particle data output file + string filename = "particleData.txt"; + std::ofstream outputFile; + + bool multiOutput = false; + config["Output"]["MultiOutput"].read(multiOutput); + clout.setMultiOutput(multiOutput); + + std::string olbdir, outputdir; + config["Application"]["OlbDir"].read(olbdir); + config["Output"]["OutputDir"].read(outputdir); + singleton::directories().setOlbDir(olbdir); + singleton::directories().setOutputDir(outputdir); + + /// Instantiation of a cuboidGeometry with weights +#ifdef PARALLEL_MODE_MPI + const int noOfCuboids = singleton::mpi().getSize(); +#else + const int noOfCuboids = 7; +#endif + + // read physical simulation times from xml-file + T physFluidNST; // time for fluid simulation + config["Application"]["SimulationTimes"]["PhysFluidNST"].read(physFluidNST); + + T physStartScT; // factor * physFluidNST; + config["Application"]["SimulationTimes"]["PhysStartScT"].read(physStartScT); + + T physParticleT; + config["Application"]["SimulationTimes"]["PhysParticleT"].read(physParticleT); + + // convert physical times to lattice time steps + int itFluidNST = converter->getLatticeTime(physFluidNST); + int itStartScaleT = converter->getLatticeTime(physStartScT); + int itParticleT = converterParticles.getLatticeTime(physParticleT); + + int itConsoleOutputFluid = itFluidNST / 10. , itVtkOutputFluid = itFluidNST / 4.; + int itConsoleOutputMagParticles = itParticleT / 75., itVtkOutputMagParticles = itParticleT / 249.; + + clout << "Fluid: physFluidNST = " << physFluidNST << "s itFluidNST = " + << itFluidNST << std::endl; + + clout << "Particles: physParticleT = " << physParticleT << "s itParticleT = " + << itParticleT << std::endl; + + clout << "noOfCuboids = " << noOfCuboids << std::endl; + + /// === 2nd Step: Prepare Geometry === + + // A: Cuboid as Channel + std::vector < T > center(3, T()); + std::vector length(3, T()); + length[0] = geometryLengthX; + length[1] = geometryLengthY; + length[2] = geometryLengthZ; + + IndicatorCuboid3D cuboid(length, center); + IndicatorLayer3D extendedDomainCuboid(cuboid, converter->getConversionFactorLength()); + + std::vector centerW(3, T()); + centerW[0] = geometryLengthX - 0.675e-3; + centerW[1] = geometryLengthY / 2; + centerW[2] = geometryLengthZ / 2.; + + std::vector normalW(3, T()); + normalW[0] = T(0); + normalW[1] = T(0); + normalW[2] = T(1); + wLength += 2 * converter->getConversionFactorLength(); + IndicatorCylinder3D wire(centerW, normalW, wRadius, wLength); + + CuboidGeometry3D cuboidGeometry(extendedDomainCuboid, + converter->getConversionFactorLength(), noOfCuboids); + + /// Instantiation of a loadBalancer + HeuristicLoadBalancer loadBalancer(cuboidGeometry); + + /// Instantiation of a superGeometry + SuperGeometry3D superGeometry(cuboidGeometry, loadBalancer, 2); + + prepareGeometry(*converter, wire, cuboid, extendedDomainCuboid, + superGeometry); + + /// === 3rd Step: Prepare Lattice === + + SuperLattice3D sLattice(superGeometry); + + BGKdynamics bulkDynamics((1. / converter->getLatticeRelaxationTime()), + instances::getBulkMomenta()); + + // choose between local and non-local boundary condition + sOnLatticeBoundaryCondition3D sOnBoundaryCondition(sLattice); + createInterpBoundaryCondition3D(sOnBoundaryCondition); + + // for the velocity field around the wire + sOffLatticeBoundaryCondition3D sOffBoundaryCondition(sLattice); + createBouzidiBoundaryCondition3D(sOffBoundaryCondition); + + // gives dynamics to cells + prepareLattice(sLattice, *converter, wire, bulkDynamics, sOnBoundaryCondition, + sOffBoundaryCondition, superGeometry); + + /// === 4th Step: Prepare Lagrange Particles + + SuperParticleSystem3D spSys(cuboidGeometry, loadBalancer, + superGeometry); + + // contact detection + NanoflannContact nanoflannContact( + *(spSys.getParticleSystems()[0]), 6. * pRadius); + spSys.setContactDetection(nanoflannContact); + + // magnetic force field from wire + std::vector < T > origin = superGeometry.getStatistics().getCenterPhysR(1); + origin[0] = geometryLengthX - 0.675e-3; + origin[2] = superGeometry.getStatistics().getMinPhysR(1)[2]; + + std::vector orientation(3, T()); + orientation[0] = 1.; + T degree = 0.; + + // car2cyl helps to know which cylinder coordinates fit to the cartesion ones + CartesianToCylinder3D car2cyl(origin, degree, orientation); + + /// Forces of Lagrange particles + + // Object of magnetic force to know magnetic value at any place in geometry + // but needs to be called by operator function. + // inhomogeneous magnetic field round cylinder + MagneticForceFromCylinder3D magForceFromCylOnParticle(car2cyl, wLength, + wRadius, 4., wMag, pMag, pRadius); + MagneticFieldFromCylinder3D magFieldFromCylOnParticle(car2cyl, wLength, + wRadius, 4., wMag); + + // external magnetic force from the wire + auto magneticPForceOne = make_shared + < MagneticForceForMagP3D + > (magForceFromCylOnParticle, magFieldFromCylOnParticle, 1.e0); + spSys.addForce(magneticPForceOne); + + // interparticular magnetic force + auto interpMagF = make_shared + < InterpMagForceForMagP3D + > (1.e0, 1.e0); + spSys.addForce(interpMagF); + + // dynamic viscosity + T dynVisc = converter->getPhysViscosity() * converter->getPhysDensity(); + + // damping force for the particle rotation + auto rotDampingF = make_shared + < LinearDampingForceForMagDipoleMoment3D + > (dynVisc, 1.e0); + spSys.addForce(rotDampingF); + + // mechanic contact force + auto hertzMindlinContF = make_shared + < HertzMindlinDeresiewicz3D + > (shearModulus, shearModulus, poissonRatio, poissonRatio, 1.e0, 1.e0, false); + spSys.addForce(hertzMindlinContF); + + // stokes drag force + SuperLatticeInterpPhysVelocity3D getVel( sLattice, *converter ); + auto stokesDragF = make_shared + < StokesDragForce3D + > (getVel, converterParticles); + spSys.addForce(stokesDragF); + + /// Boundarys + + T dT = converter->getConversionFactorTime() ; + std::set wireBMaterial = { 5, 4 }; + std::set reflBMat = { 2, 3 }; + + // particle boundary for deposited particles + auto wireBoundary = make_shared + < WireBoundaryForMagP3D + > (superGeometry, wireBMaterial); + spSys.addBoundary(wireBoundary); + + // particle boundary + auto materialreflectBoundary = make_shared + < SimpleReflectBoundary3D + > (dT, superGeometry, reflBMat); + spSys.addBoundary(materialreflectBoundary); + + // cuboid overlap + spSys.setOverlap(2.*converter->getConversionFactorLength()); + + // particles vtu output + std::string particleOutputName = "vtuOutMagP"; + std::bitset < 9 > properties; + properties.set(); + + SuperParticleSysVtuWriterMag particleOut(spSys, particleOutputName, properties); + + clout << "Number of Forces: " << spSys.numOfForces()[0] << std::endl; + + /// Paramters for magnetic particle generation + std::vector pPos = { 0., 0., 0. }; // position + std::vector pVel = { 0., 0., 0. }; // velocity + std::vector pAVel = { 0., 0., 0. }; // angular velocity + std::vector pTrq = { 0., 0., 0. }; // torque + std::vector pDMoment = { -1., 0., 0. }; // orientation mag. dipole moment + if(!util::nearZero(util::norm(pDMoment))){util::normalize(pDMoment);} + else {clout << "Norm of pDMoment near zero!" << endl; exit(0);} + + // magnetic particle magPartTemplate is used as copy template + PARTICLE magPartTemplate(pPos, pVel, pMass, pRadius, 0, pDMoment, pAVel, pTrq, pMag, 1); + + // intialization geometry particles + std::vector particlesOrigin(3, T(0)); + particlesOrigin = superGeometry.getStatistics().getMinPhysR(1); + std::vector particlesExtend(3, T(0)); + + particlesOrigin[0] = 1.e-4; + particlesOrigin[1] = (superGeometry.getStatistics().getMaxPhysR(1)[1]) * (1./2. - 1./6.); + particlesOrigin[2] = (superGeometry.getStatistics().getMaxPhysR(1)[2]) * (1./2. - 1./6.); + particlesExtend[0] = (superGeometry.getStatistics().getMaxPhysR(1)[0]) * 1./15. ; + particlesExtend[1] = (superGeometry.getStatistics().getMaxPhysR(1)[1]) * 2./6. ; + particlesExtend[2] = (superGeometry.getStatistics().getMaxPhysR(1)[2]) * 2./6. ; + IndicatorCuboid3D cuboidPart(particlesExtend, particlesOrigin); + + clout << "starting simulation" << endl; + + /// === 5th Step: Main Loop with Timer === + + // is set on true for the appropriate loop of fluid / particle simulation + bool outNS = false; + bool outLP = false; + + // === Calculating Navier-Stokes Fluid + + // if there is no computed fluid data in a external data files + // compute fluid data in data file, else load it and jump to next loop + if (!(sLattice.load("water"))) { + + clout << std::endl; + clout << "Fluid data has to be computed " << std::endl; + clout << std::endl; + + outNS = true; + + Timer timerFluid(itFluidNST, superGeometry.getStatistics().getNvoxel()); + timerFluid.start(); + + for (int iT = 0; iT < itFluidNST; ++iT) { + + setBoundaryValues(sLattice, *converter, spSys, iT, itStartScaleT, + superGeometry, outNS, outLP); + + sLattice.collideAndStream(); + + getResults(superGeometry, sLattice, *converter, + magForceFromCylOnParticle, magFieldFromCylOnParticle, spSys, + particleOut, timerFluid, iT, itConsoleOutputFluid, itVtkOutputFluid, + itConsoleOutputMagParticles, itVtkOutputMagParticles, outNS, outLP); + } + + outNS = false; + + // save computed fluid data in external data files + clout << "save fluid data " << std::endl; + + sLattice.communicate(); + sLattice.save("water"); + + timerFluid.stop(); + timerFluid.printSummary(); + + } else { + sLattice.communicate(); + sLattice.collideAndStream(); + sLattice.getStatistics().print(iT, converter->getPhysTime(iT)); + } + + // === Calculating Particles + + clout << std::endl; + clout << "particles computation" << std::endl; + clout << std::endl; + + outLP = true; + + // number of Particles in initialization geometry + int noOfPartX = 4 ; + int noOfPartY = 12 ; + int noOfPartZ = 6 ; + int noOfPart = noOfPartX * noOfPartY * noOfPartZ; + + Timer timerParticles(itParticleT, noOfPart); + timerParticles.start(); + + setBoundaryValues(sLattice, converterParticles, spSys, iT, itStartScaleT, + superGeometry, outNS, outLP); + + for (int iT = 0; iT < itParticleT; ++iT) { + + // particles initialization + if (iT == 0) { + + spSys.addParticleEquallyDistributed(cuboidPart, noOfPartX, noOfPartY, noOfPartZ, magPartTemplate); + spSys.setParticlesPosRandom(0.75 * pRadius, 0.75 * pRadius, 0.55 * pRadius); + } + + getResults(superGeometry, sLattice, converterParticles, + magForceFromCylOnParticle, magFieldFromCylOnParticle, spSys, + particleOut, timerParticles, iT, itConsoleOutputFluid, itVtkOutputFluid, + itConsoleOutputMagParticles, itVtkOutputMagParticles, + outNS, outLP); + + // particles simulation + spSys.simulate(converterParticles.getConversionFactorTime()); + + } + + timerParticles.stop(); + timerParticles.printSummary(); + + clout << "End of simulation" << std::endl; +} diff --git a/examples/particles/magneticParticles3d/magneticParticles3d.xml b/examples/particles/magneticParticles3d/magneticParticles3d.xml new file mode 100644 index 0000000..8235e9f --- /dev/null +++ b/examples/particles/magneticParticles3d/magneticParticles3d.xml @@ -0,0 +1,47 @@ + + + sedimentation + 3 + ../../../ + + 0.75 + 50 + + + .001 + 0.002 + Kinematic viscosity + 1.e-6 + 1000. + 0. + + + .12 + .1 + 1.55 + + + + ./tmp/ + false + + magneticParticles3d.log + 1. + + + magneticParticles3dvtk + 0.01 + + + magneticParticles3dimage + 0.01 + + + 0 + 0 + 0.1 + 0.1 + true + + + diff --git a/examples/particles/magneticParticles3d/module.mk b/examples/particles/magneticParticles3d/module.mk new file mode 100644 index 0000000..1190482 --- /dev/null +++ b/examples/particles/magneticParticles3d/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)))) -- cgit v1.2.3