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 --- ...AdvectionDiffusionMRTCouplingPostProcessor3D.hh | 281 +++++++++++++++++++++ 1 file changed, 281 insertions(+) create mode 100644 src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh (limited to 'src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh') diff --git a/src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh b/src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh new file mode 100644 index 0000000..a32d22f --- /dev/null +++ b/src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh @@ -0,0 +1,281 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani, Jonas Latt + * 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. + */ + +#ifndef NAVIER_STOKES_INTO_ADVECTION_DIFFUSION_MRT_COUPLING_POST_PROCESSOR_3D_HH +#define NAVIER_STOKES_INTO_ADVECTION_DIFFUSION_MRT_COUPLING_POST_PROCESSOR_3D_HH + +#include "latticeDescriptors.h" +#include "navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.h" +#include "core/blockLattice3D.h" +#include "core/util.h" +#include "core/finiteDifference3D.h" +#include "advectionDiffusionForces.hh" +#include "advectionDiffusionForces.h" + +namespace olb { + +//===================================================================================== +//============== NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D =========== +//===================================================================================== +template +NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D:: +NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, + T gravity_, T T0_, T deltaTemp_, std::vector dir_, + std::vector partners_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), + dir(dir_), partners(partners_) +{ + // we normalize the direction of force vector + T normDir = T(); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + normDir += dir[iD]*dir[iD]; + } + normDir = sqrt(normDir); + for (unsigned iD = 0; iD < dir.size(); ++iD) { + dir[iD] /= normDir; + } +} + +template +void NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + typedef DESCRIPTOR L; + enum {x,y,z}; + enum { + velOffset = descriptors::AdvectionDiffusionMRTD3Q7Descriptor::template index(), + forceOffset = DESCRIPTOR::template index() + }; + + BlockLattice3D *tPartner = + dynamic_cast *>(partners[0]); + + int newX0, newX1, newY0, newY1, newZ0, newZ1; + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + // Velocity coupling + T *u = tPartner->get(iX,iY,iZ).template getFieldPointer(); + blockLattice.get(iX,iY,iZ).computeU(u); + + //coupling between the temperature and navier stokes. + + T *force = blockLattice.get(iX,iY,iZ).template getFieldPointer(); + T temperature = tPartner->get(iX,iY,iZ).computeRho(); + T rho = blockLattice.get(iX,iY,iZ).computeRho(); + for (unsigned iD = 0; iD < L::d; ++iD) { + force[iD] = gravity * rho * (temperature - T0) / deltaTemp * dir[iD]; + } + } + } + } + } +} + +template +void NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +//===================================================================================== +//============== StokesDragCouplingPostProcessor3D =========== +//===================================================================================== + +template +AdvectionDiffusionParticleMRTCouplingPostProcessor3D:: +AdvectionDiffusionParticleMRTCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int iC_, + int offset_, std::vector partners_, + std::vector > > forces_) + : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), iC(iC_), offset(offset_), forces(forces_) +{ + BlockLattice3D *partnerLattice = + dynamic_cast *>(partners_[0]); + adCell = &partnerLattice->get(x0,y0,z0); + vel = partnerLattice->get(x0,y0,z0)[offset]; + vel_new = partnerLattice->get(x0,y0,z0)[offset]; + velXp = partnerLattice->get(x0+1,y0,z0)[offset]; + velXn = partnerLattice->get(x0-1,y0,z0)[offset]; + velYp = partnerLattice->get(x0,y0+1,z0)[offset]; + velYn = partnerLattice->get(x0,y0-1,z0)[offset]; + velZp = partnerLattice->get(x0,y0,z0+1)[offset]; + velZn = partnerLattice->get(x0,y0,z0-1)[offset]; +} + +template +void AdvectionDiffusionParticleMRTCouplingPostProcessor3D:: +processSubDomain(BlockLattice3D& blockLattice, + int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) +{ + int newX0, newX1, newY0, newY1, newZ0, newZ1; + + int off = (par) ? 3 : 0; + int off2 = (par) ? 0 : 3; + + if ( util::intersect ( + x0, x1, y0, y1, z0, z1, + x0_, x1_, y0_, y1_, z0_, z1_, + newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) { + + for (int iX=newX0; iX<=newX1; ++iX) { + for (int iY=newY0; iY<=newY1; ++iY) { + for (int iZ=newZ0; iZ<=newZ1; ++iZ) { + + int latticeR[4] = {iC, iX, iY, iZ}; + T velGrad[3] = {0.,0.,0.}; + T forceValue[3] = {0.,0.,0.}; + T velF[3] = {0.,0.,0.}; + + nsCell = &(blockLattice.get(iX,iY,iZ)); + //.computeU(velF); + + if (forces.begin() != forces.end()) { + // calculating upwind Gradient + if (vel[0+off]<0.) { + velGrad[0] = vel[0+off]*(velXp[0+off]-vel[0+off]); + velGrad[1] = vel[0+off]*(velXp[1+off]-vel[1+off]); + velGrad[2] = vel[0+off]*(velXp[2+off]-vel[2+off]); + } else { + velGrad[0] = vel[0+off]*(vel[0+off]-velXn[0+off]); + velGrad[1] = vel[0+off]*(vel[1+off]-velXn[1+off]); + velGrad[2] = vel[0+off]*(vel[2+off]-velXn[2+off]); + } + if (vel[1+off]<0.) { + velGrad[0] += vel[1+off]*(velYp[0+off]-vel[0+off]); + velGrad[1] += vel[1+off]*(velYp[1+off]-vel[1+off]); + velGrad[2] += vel[1+off]*(velYp[2+off]-vel[2+off]); + } else { + velGrad[0] += vel[1+off]*(vel[0+off]-velYn[0+off]); + velGrad[1] += vel[1+off]*(vel[1+off]-velYn[1+off]); + velGrad[2] += vel[1+off]*(vel[2+off]-velYn[2+off]); + } + if (vel[2+off]<0.) { + velGrad[0] += vel[2+off]*(velZp[0+off]-vel[0+off]); + velGrad[1] += vel[2+off]*(velZp[1+off]-vel[1+off]); + velGrad[2] += vel[2+off]*(velZp[2+off]-vel[2+off]); + } else { + velGrad[0] += vel[2+off]*(vel[0+off]-velZn[0+off]); + velGrad[1] += vel[2+off]*(vel[1+off]-velZn[1+off]); + velGrad[2] += vel[2+off]*(vel[2+off]-velZn[2+off]); + } + + for (AdvectionDiffusionForce3D& f : forces) { + f.applyForce(forceValue, nsCell, adCell, &vel[off], latticeR); + } + + // compute new particle velocity + for (int i=0; i < DESCRIPTOR::d; i++) { + vel_new[i+off2] = vel[i+off] + forceValue[i] - velGrad[i]; + } + } else { + nsCell->computeU(velF); + T scaleFactor = T(1); + for (int i = 0; i < DESCRIPTOR::d; i++) { + vel_new[i + off2] = velF[i] * scaleFactor; + } + } + } + } + } + } + par = !par; +} + +template +void AdvectionDiffusionParticleMRTCouplingPostProcessor3D:: +process(BlockLattice3D& blockLattice) +{ + processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1); +} + +// LatticeCouplingGenerator for advectionDiffusion coupling + +template +NavierStokesAdvectionDiffusionMRTCouplingGenerator3D:: +NavierStokesAdvectionDiffusionMRTCouplingGenerator3D(int x0_, int x1_, int y0_, int y1_,int z0_, int z1_, + T gravity_, T T0_, T deltaTemp_, std::vector dir_) + : LatticeCouplingGenerator3D(x0_, x1_, y0_, y1_, z0_, z1_), + gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_) +{ } + +template +PostProcessor3D* NavierStokesAdvectionDiffusionMRTCouplingGenerator3D::generate ( + std::vector partners) const +{ + return new NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D( + this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, gravity, T0, deltaTemp, dir,partners); +} + +template +LatticeCouplingGenerator3D* NavierStokesAdvectionDiffusionMRTCouplingGenerator3D::clone() const +{ + return new NavierStokesAdvectionDiffusionMRTCouplingGenerator3D(*this); +} + +// LatticeCouplingGenerator for one-way advectionDiffusion coupling with Stokes drag + +template +AdvectionDiffusionParticleMRTCouplingGenerator3D:: +AdvectionDiffusionParticleMRTCouplingGenerator3D(int offset_) + : LatticeCouplingGenerator3D(0, 0, 0, 0, 0, 0), offset(offset_) +{ } + +template +PostProcessor3D* AdvectionDiffusionParticleMRTCouplingGenerator3D::generate ( + std::vector partners) const +{ + return new AdvectionDiffusionParticleMRTCouplingPostProcessor3D(this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, this->iC, offset, partners, ADforces); +} + +template +LatticeCouplingGenerator3D* AdvectionDiffusionParticleMRTCouplingGenerator3D::clone() const +{ + return new AdvectionDiffusionParticleMRTCouplingGenerator3D(*this); +} + +template +void AdvectionDiffusionParticleMRTCouplingGenerator3D::addForce( + AdvectionDiffusionForce3D &force) +{ + ADforces.push_back(force); +} + + +} // namespace olb + +#endif -- cgit v1.2.3