OpenLB 1.7
Loading...
Searching...
No Matches
shanChenForcedSingleComponentPostProcessor3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani,
4 * Jonas Latt, 2013 Mathias J. Krause
5 * E-mail contact: info@openlb.net
6 * The most recent release of OpenLB can be downloaded at
7 * <http://www.openlb.net/>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public
20 * License along with this program; if not, write to the Free
21 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 * Boston, MA 02110-1301, USA.
23*/
24
25#ifndef SHAN_CHEN_FORCED_SINGLE_COMPONENT_POST_PROCESSOR_3D_HH
26#define SHAN_CHEN_FORCED_SINGLE_COMPONENT_POST_PROCESSOR_3D_HH
27
30#include "core/util.h"
32
33namespace olb {
34
36
37
38template<typename T, typename DESCRIPTOR>
39ShanChenForcedSingleComponentPostProcessor3D <T,DESCRIPTOR>::
40ShanChenForcedSingleComponentPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
41 T G_, std::vector<T> rho0_, AnalyticalF<1,T,T>& iP_,
42 std::vector<BlockStructureD<3>*> partners_)
43 : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), G(G_), rho0(rho0_), interactionPotential(iP_), partners(partners_)
44{
45 this->getName() = "ShanChenForcedSingleComponentPostProcessor3D";
46}
47
48template<typename T, typename DESCRIPTOR>
49ShanChenForcedSingleComponentPostProcessor3D <T,DESCRIPTOR>::
50ShanChenForcedSingleComponentPostProcessor3D(T G_, std::vector<T> rho0_, AnalyticalF<1,T,T>& iP_,
51 std::vector<BlockStructureD<3>*> partners_)
52 : x0(0), x1(0), y0(0), y1(0), z0(0), z1(0), G(G_), rho0(rho0_), interactionPotential(iP_), partners(partners_)
53{
54 this->getName() = "ShanChenForcedSingleComponentPostProcessor3D";
55}
56
57template<typename T, typename DESCRIPTOR>
60 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_ )
61{
62 typedef DESCRIPTOR L;
63
64 int newX0, newX1, newY0, newY1, newZ0, newZ1;
65 if ( util::intersect ( x0, x1, y0, y1, z0, z1,
66 x0_, x1_, y0_, y1_, z0_, z1_,
67 newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
68
69 auto& rhoField = blockLattice.template getField<RHO_CACHE>()[0];
70
71 // Compute density and velocity on every site of first lattice, and store result
72 // in external scalars; envelope cells are included, because they are needed
73 // to compute the interaction potential in what follows.
74 for (int iX=newX0-1; iX<=newX1+1; ++iX) {
75 for (int iY=newY0-1; iY<=newY1+1; ++iY) {
76 for (int iZ=newZ0-1; iZ<=newZ1+1; ++iZ) {
77 Cell<T,DESCRIPTOR> cell = blockLattice.get(iX,iY,iZ);
78 rhoField[cell.getCellId()] = cell.computeRho()*rho0[0];
79 }
80 }
81 }
82
83 for (int iX=newX0; iX<=newX1; ++iX) {
84 for (int iY=newY0; iY<=newY1; ++iY) {
85 for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
86 Cell<T,DESCRIPTOR> blockCell = blockLattice.get(iX,iY,iZ);
87
88 auto j = blockCell.template getField<descriptors::VELOCITY>();
89 lbm<DESCRIPTOR>::computeJ(blockCell,j);
90 blockCell.template setField<descriptors::VELOCITY>(j);
91
92 T blockOmega = blockCell.getDynamics()->getParameters(blockLattice).template getOrFallback<descriptors::OMEGA>(0);
93
94 // Computation of the common velocity, shared among the two populations
95 T rhoTot = rhoField[blockCell.getCellId()]*blockOmega;
96
97 Vector<T, 3> uTot;
98 auto blockU = blockCell.template getFieldPointer<descriptors::VELOCITY>(); // contains precomputed value rho*u
99 uTot = (blockU*rho0[0]*blockOmega) / rhoTot;
100
101 // Computation of the interaction potential
102 Vector<T, 3> rhoBlockContribution;
103 T psi;
104 interactionPotential(&psi, &rhoField[blockCell.getCellId()]);
105 for (int iPop = 0; iPop < L::q; ++iPop) {
106 int nextX = iX + descriptors::c<L>(iPop,0);
107 int nextY = iY + descriptors::c<L>(iPop,1);
108 int nextZ = iZ + descriptors::c<L>(iPop,2);
109 T blockRho;
110 interactionPotential(&blockRho, &rhoField[blockLattice.getCellId(nextX, nextY, nextZ)]);
111 rhoBlockContribution += psi * blockRho * descriptors::c<L>(iPop)* descriptors::t<T,L>(iPop);
112 }
113
114 // Computation and storage of the final velocity, consisting
115 // of u and the momentum difference due to interaction
116 // potential plus external force
117 auto externalBlockForce = blockCell.template getFieldPointer<descriptors::EXTERNAL_FORCE>();
118
119 blockCell.template setField<descriptors::VELOCITY>(uTot);
120 blockCell.template setField<descriptors::FORCE>(externalBlockForce
121 - G*rhoBlockContribution/rhoField[blockCell.getCellId()]);
122 }
123 }
124 }
125 }
126}
127
128template<typename T, typename DESCRIPTOR>
131{
132 processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
133}
134
135
137
138template<typename T, typename DESCRIPTOR>
140 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T G_, std::vector<T> rho0_, AnalyticalF<1,T,T>& iP_)
141 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_), G(G_), rho0(rho0_), interactionPotential(iP_)
142{ }
143
144template<typename T, typename DESCRIPTOR>
146 T G_, std::vector<T> rho0_, AnalyticalF<1,T,T>& iP_)
147 : LatticeCouplingGenerator3D<T,DESCRIPTOR>(0, 0, 0, 0, 0, 0), G(G_), rho0(rho0_), interactionPotential(iP_)
148{ }
149
150template<typename T, typename DESCRIPTOR>
152 std::vector<BlockStructureD<3>*> partners) const
153{
155 this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, G, rho0, interactionPotential, partners);
156}
157
158template<typename T, typename DESCRIPTOR>
163
164
165
166} // namespace olb
167
168#endif
AnalyticalF are applications from DD to XD, where X is set by the constructor.
Platform-abstracted block lattice for external access and inter-block interaction.
Cell< T, DESCRIPTOR > get(CellID iCell)
Get Cell interface for index iCell.
Base of a regular block.
CellID getCellId(LatticeR< D > latticeR) const
Get 1D cell ID.
Highest-level interface to Cell data.
Definition cell.h:148
Dynamics< T, DESCRIPTOR > * getDynamics()
Get a pointer to the dynamics.
Definition cell.hh:174
T computeRho() const
Compute particle density on the cell.
Definition cell.hh:206
std::size_t getCellId() const
Return memory ID of the currently represented cell.
Definition cell.hh:51
std::string & getName()
read and write access to name
PostProcessor3D< T, DESCRIPTOR > * generate(std::vector< BlockStructureD< 3 > * > partners) const override
LatticeCouplingGenerator3D< T, DESCRIPTOR > * clone() const override
ShanChenForcedSingleComponentGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, T G_, std::vector< T > rho0_, AnalyticalF< 1, T, T > &iP_)
LatticeCouplingGenerator for NS coupling.
Multiphysics class for coupling between different lattices.
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
Plain old scalar vector.
Definition vector.h:47
bool intersect(int x0, int x1, int y0, int y1, int x0_, int x1_, int y0_, int y1_, int &newX0, int &newX1, int &newY0, int &newY1)
Definition util.h:89
Top level namespace for all of OpenLB.
static void computeJ(CELL &cell, J &j) any_platform
Computation of momentum.
Definition lbm.h:197
Set of functions commonly used in LB computations – header file.