dynamicsProcessor3D.hh

00001 /* This file is part of the Palabos library.
00002  *
00003  * Copyright (C) 2011 FlowKit Sarl
00004  * Avenue de Chailly 23
00005  * 1012 Lausanne, Switzerland
00006  * E-mail contact: contact@flowkit.com
00007  *
00008  * The most recent release of Palabos can be downloaded at 
00009  * <http://www.palabos.org/>
00010  *
00011  * The library Palabos is free software: you can redistribute it and/or
00012  * modify it under the terms of the GNU Affero General Public License as
00013  * published by the Free Software Foundation, either version 3 of the
00014  * License, or (at your option) any later version.
00015  *
00016  * The library is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU Affero General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU Affero General Public License
00022  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00023 */
00024 
00028 #ifndef DYNAMICS_PROCESSOR_3D_HH
00029 #define DYNAMICS_PROCESSOR_3D_HH
00030 
00031 #include "basicDynamics/dynamicsProcessor3D.h"
00032 #include "core/cell.h"
00033 #include "latticeBoltzmann/geometricOperationTemplates.h"
00034 #include "latticeBoltzmann/latticeTemplates.h"
00035 #include "atomicBlock/blockLattice3D.h"
00036 #include "multiGrid/multiGridUtil.h"
00037 #include "core/plbProfiler.h"
00038 
00039 namespace plb {
00040 
00041 /* ************* Class ExternalRhoJcollideAndStream3D ******************* */
00042 
00043 template<typename T, template<typename U> class Descriptor>
00044 void ExternalRhoJcollideAndStream3D<T,Descriptor>::collide (
00045         BlockLattice3D<T,Descriptor>& lattice, Box3D const& domain,
00046         ScalarField3D<T> const& rhoBarField, Dot3D const& offset1,
00047         TensorField3D<T,3> const& jField, Dot3D const& offset2, BlockStatistics& stat )
00048 {
00049     for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
00050         for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
00051             for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
00052                 Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
00053                 T rhoBar            = rhoBarField.get(iX+offset1.x, iY+offset1.y, iZ+offset1.z);
00054                 Array<T,3> const& j = jField.get(iX+offset2.x, iY+offset2.y, iZ+offset2.z);
00055                 cell.getDynamics().collide(cell, rhoBar, j, T(), stat);
00056                 cell.revert();
00057             }
00058         }
00059     }
00060 }
00061 
00062 template<typename T, template<typename U> class Descriptor>
00063 void ExternalRhoJcollideAndStream3D<T,Descriptor>::bulkCollideAndStream (
00064         BlockLattice3D<T,Descriptor>& lattice, Box3D const& domain,
00065         ScalarField3D<T> const& rhoBarField, Dot3D const& offset1,
00066         TensorField3D<T,3> const& jField, Dot3D const& offset2, BlockStatistics& stat )
00067 {
00068     for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
00069         for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
00070             for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
00071                 Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
00072                 T rhoBar            = rhoBarField.get(iX+offset1.x, iY+offset1.y, iZ+offset1.z);
00073                 Array<T,3> const& j = jField.get(iX+offset2.x, iY+offset2.y, iZ+offset2.z);
00074                 cell.getDynamics().collide(cell, rhoBar, j, T(), stat);
00075                 latticeTemplates<T,Descriptor>::swapAndStream3D(lattice.grid, iX, iY, iZ);
00076             }
00077         }
00078     }
00079 }
00080 
00081 
00082 template<typename T, template<typename U> class Descriptor>
00083 void ExternalRhoJcollideAndStream3D<T,Descriptor>::boundaryStream (
00084         BlockLattice3D<T,Descriptor>& lattice,
00085         Box3D const& bound, Box3D const& domain )
00086 {
00087     // Make sure domain is contained within bound
00088     PLB_PRECONDITION( contained(domain, bound) );
00089 
00090     for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
00091         for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
00092             for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
00093                 for (plint iPop=1; iPop<=Descriptor<T>::q/2; ++iPop) {
00094                     plint nextX = iX + Descriptor<T>::c[iPop][0];
00095                     plint nextY = iY + Descriptor<T>::c[iPop][1];
00096                     plint nextZ = iZ + Descriptor<T>::c[iPop][2];
00097                     if ( nextX>=bound.x0 && nextX<=bound.x1 && nextY>=bound.y0 && nextY<=bound.y1 &&
00098                          nextZ>=bound.z0 && nextZ<=bound.z1 )
00099                     {
00100                         std::swap(lattice.grid[iX][iY][iZ][iPop+Descriptor<T>::q/2],
00101                                   lattice.grid[nextX][nextY][nextZ][iPop]);
00102                     }
00103                 }
00104             }
00105         }
00106     }
00107 }
00108 
00109 template<typename T, template<typename U> class Descriptor>
00110 void ExternalRhoJcollideAndStream3D<T,Descriptor>::processGenericBlocks (
00111         Box3D domain, std::vector<AtomicBlock3D*> atomicBlocks )
00112 {
00113     global::timer("collideAndStream").start();
00114     BlockLattice3D<T,Descriptor>& lattice =
00115         dynamic_cast<BlockLattice3D<T,Descriptor>&>(*atomicBlocks[0]);
00116     ScalarField3D<T> const& rhoBarField =
00117         dynamic_cast<ScalarField3D<T> const&>(*atomicBlocks[1]);
00118     TensorField3D<T,3> const& jField =
00119         dynamic_cast<TensorField3D<T,3> const&>(*atomicBlocks[2]);
00120 
00121     BlockStatistics stat;
00122     stat.subscribeAverage();
00123     stat.subscribeAverage();
00124     stat.subscribeMax();
00125 
00126     static const plint vicinity = Descriptor<T>::vicinity;
00127     Box3D extDomain(domain.enlarge(vicinity));
00128 
00129     Dot3D offset1 = computeRelativeDisplacement(lattice, rhoBarField);
00130     Dot3D offset2 = computeRelativeDisplacement(lattice, jField);
00131 
00132     global::profiler().start("collStream");
00133     global::profiler().increment("collStreamCells", extDomain.nCells());
00134 
00135     // First, do the collision on cells within a boundary envelope of width
00136     // equal to the range of the lattice vectors (e.g. 1 for D2Q9)
00137     collide(lattice,
00138             Box3D(extDomain.x0,extDomain.x0+vicinity-1,
00139                   extDomain.y0,extDomain.y1,
00140                   extDomain.z0,extDomain.z1),
00141             rhoBarField, offset1, jField, offset2, stat);
00142     collide(lattice,
00143             Box3D(extDomain.x1-vicinity+1,extDomain.x1,
00144                   extDomain.y0,extDomain.y1,
00145                   extDomain.z0,extDomain.z1),
00146             rhoBarField, offset1, jField, offset2, stat);
00147     collide(lattice,
00148             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00149                   extDomain.y0,extDomain.y0+vicinity-1,
00150                   extDomain.z0,extDomain.z1),
00151             rhoBarField, offset1, jField, offset2, stat);
00152     collide(lattice,
00153             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00154                   extDomain.y1-vicinity+1,extDomain.y1,
00155                   extDomain.z0,extDomain.z1),
00156             rhoBarField, offset1, jField, offset2, stat);
00157     collide(lattice,
00158             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00159                   extDomain.y0+vicinity,extDomain.y1-vicinity,
00160                   extDomain.z0,extDomain.z0+vicinity-1),
00161             rhoBarField, offset1, jField, offset2, stat);
00162     collide(lattice,
00163             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00164                   extDomain.y0+vicinity,extDomain.y1-vicinity,
00165                   extDomain.z1-vicinity+1,extDomain.z1),
00166             rhoBarField, offset1, jField, offset2, stat);
00167 
00168     // Then, do the efficient collideAndStream algorithm in the bulk,
00169     // excluding the envelope (this is efficient because there is no
00170     // if-then-else statement within the loop, given that the boundary
00171     // region is excluded)
00172     bulkCollideAndStream(lattice,
00173                          Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00174                                extDomain.y0+vicinity,extDomain.y1-vicinity,
00175                                extDomain.z0+vicinity,extDomain.z1-vicinity),
00176                          rhoBarField, offset1, jField, offset2, stat);
00177 
00178     // Finally, do streaming in the boundary envelope to conclude the
00179     // collision-stream cycle
00180     boundaryStream(lattice, extDomain, Box3D(extDomain.x0,extDomain.x0+vicinity-1,
00181                                              extDomain.y0,extDomain.y1,
00182                                              extDomain.z0,extDomain.z1) );
00183     boundaryStream(lattice, extDomain, Box3D(extDomain.x1-vicinity+1,extDomain.x1,
00184                                              extDomain.y0,extDomain.y1,
00185                                              extDomain.z0,extDomain.z1) );
00186     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00187                                              extDomain.y0,extDomain.y0+vicinity-1,
00188                                              extDomain.z0,extDomain.z1) );
00189     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00190                                              extDomain.y1-vicinity+1,extDomain.y1,
00191                                              extDomain.z0,extDomain.z1) );
00192     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00193                                              extDomain.y0+vicinity,extDomain.y1-vicinity,
00194                                              extDomain.z0,extDomain.z0+vicinity-1) );
00195     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00196                                              extDomain.y0+vicinity,extDomain.y1-vicinity,
00197                                              extDomain.z1-vicinity+1,extDomain.z1) );
00198     global::profiler().stop("collStream");
00199     global::timer("collideAndStream").stop();
00200 }
00201 
00202 template<typename T, template<typename U> class Descriptor>
00203 ExternalRhoJcollideAndStream3D<T,Descriptor>*
00204     ExternalRhoJcollideAndStream3D<T,Descriptor>::clone() const
00205 {
00206     return new ExternalRhoJcollideAndStream3D<T,Descriptor>(*this);
00207 }
00208 
00209 template<typename T, template<typename U> class Descriptor>
00210 void ExternalRhoJcollideAndStream3D<T,Descriptor>::getTypeOfModification (
00211         std::vector<modif::ModifT>& modified) const
00212 {
00213     modified[0] = modif::staticVariables;
00214     modified[1] = modif::nothing;
00215     modified[2] = modif::nothing;
00216 }
00217 
00218 
00219 
00220 /* ************* Class OnLinkExternalRhoJcollideAndStream3D ******************* */
00221 
00222 template<typename T, template<typename U> class Descriptor>
00223 void OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>::collide (
00224         BlockLattice3D<T,Descriptor>& lattice, Box3D const& domain,
00225         ScalarField3D<T> const& rhoBarField, Dot3D const& offset1,
00226         TensorField3D<T,3> const& jField, Dot3D const& offset2, BlockStatistics& stat )
00227 {
00228     for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
00229         for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
00230             for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
00231                 Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
00232                 T rhoBar            = rhoBarField.get(iX+offset1.x, iY+offset1.y, iZ+offset1.z);
00233                 Array<T,3> const& j = jField.get(iX+offset2.x, iY+offset2.y, iZ+offset2.z);
00234                 cell.getDynamics().collide(cell, rhoBar, j, T(), stat);
00235                 cell.revert();
00236             }
00237         }
00238     }
00239 }
00240 
00241 template<typename T, template<typename U> class Descriptor>
00242 void onLinkSwapAndStream3D( Cell<T,Descriptor> ***grid,
00243                             plint iX, plint iY, plint iZ)
00244 {
00245     static int bbId = BounceBack<T,Descriptor>().getId();
00246     const plint half = Descriptor<T>::q/2;
00247     for (plint iPop=1; iPop<=half; ++iPop) {
00248         plint nextX = iX + Descriptor<T>::c[iPop][0];
00249         plint nextY = iY + Descriptor<T>::c[iPop][1];
00250         plint nextZ = iZ + Descriptor<T>::c[iPop][2];
00251         if (grid[iX][iY][iZ].getDynamics().getId()==bbId ||
00252             grid[nextX][nextY][nextZ].getDynamics().getId()==bbId)
00253         {
00254             std::swap(grid[iX][iY][iZ][iPop],grid[iX][iY][iZ][iPop+half]);
00255         }
00256         else {
00257             T fTmp                          = grid[iX][iY][iZ][iPop];
00258             grid[iX][iY][iZ][iPop]          = grid[iX][iY][iZ][iPop+half];
00259             grid[iX][iY][iZ][iPop+half]     = grid[nextX][nextY][nextZ][iPop];
00260             grid[nextX][nextY][nextZ][iPop] = fTmp;
00261         }
00262      }
00263 }
00264 
00265 template<typename T, template<typename U> class Descriptor>
00266 void OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>::bulkCollideAndStream (
00267         BlockLattice3D<T,Descriptor>& lattice, Box3D const& domain,
00268         ScalarField3D<T> const& rhoBarField, Dot3D const& offset1,
00269         TensorField3D<T,3> const& jField, Dot3D const& offset2, BlockStatistics& stat )
00270 {
00271     for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
00272         for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
00273             for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
00274                 Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
00275                 T rhoBar            = rhoBarField.get(iX+offset1.x, iY+offset1.y, iZ+offset1.z);
00276                 Array<T,3> const& j = jField.get(iX+offset2.x, iY+offset2.y, iZ+offset2.z);
00277                 cell.getDynamics().collide(cell, rhoBar, j, T(), stat);
00278                 onLinkSwapAndStream3D<T,Descriptor>(lattice.grid, iX, iY, iZ);
00279             }
00280         }
00281     }
00282 }
00283 
00284 template<typename T, template<typename U> class Descriptor>
00285 void OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>::boundaryStream (
00286         BlockLattice3D<T,Descriptor>& lattice,
00287         Box3D const& bound, Box3D const& domain )
00288 {
00289     // Make sure domain is contained within bound
00290     PLB_PRECONDITION( contained(domain, bound) );
00291 
00292     int bbId = BounceBack<T,Descriptor>().getId();
00293     for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
00294         for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
00295             for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
00296                 for (plint iPop=1; iPop<=Descriptor<T>::q/2; ++iPop) {
00297                     plint nextX = iX + Descriptor<T>::c[iPop][0];
00298                     plint nextY = iY + Descriptor<T>::c[iPop][1];
00299                     plint nextZ = iZ + Descriptor<T>::c[iPop][2];
00300                     if ( nextX>=bound.x0 && nextX<=bound.x1 && nextY>=bound.y0 && nextY<=bound.y1 &&
00301                          nextZ>=bound.z0 && nextZ<=bound.z1 )
00302                     {
00303                         if (lattice.grid[iX][iY][iZ].getDynamics().getId()!=bbId &&
00304                             lattice.grid[nextX][nextY][nextZ].getDynamics().getId()!=bbId)
00305                         {
00306                             std::swap(lattice.grid[iX][iY][iZ][iPop+Descriptor<T>::q/2],
00307                                       lattice.grid[nextX][nextY][nextZ][iPop]);
00308                         }
00309                     }
00310                 }
00311             }
00312         }
00313     }
00314 }
00315 
00316 template<typename T, template<typename U> class Descriptor>
00317 void OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>::processGenericBlocks (
00318         Box3D domain, std::vector<AtomicBlock3D*> atomicBlocks )
00319 {
00320     global::timer("collideAndStream").start();
00321     BlockLattice3D<T,Descriptor>& lattice =
00322         dynamic_cast<BlockLattice3D<T,Descriptor>&>(*atomicBlocks[0]);
00323     ScalarField3D<T> const& rhoBarField =
00324         dynamic_cast<ScalarField3D<T> const&>(*atomicBlocks[1]);
00325     TensorField3D<T,3> const& jField =
00326         dynamic_cast<TensorField3D<T,3> const&>(*atomicBlocks[2]);
00327 
00328     BlockStatistics stat;
00329     stat.subscribeAverage();
00330     stat.subscribeAverage();
00331     stat.subscribeMax();
00332 
00333     static const plint vicinity = Descriptor<T>::vicinity;
00334     Box3D extDomain(domain.enlarge(vicinity));
00335 
00336     Dot3D offset1 = computeRelativeDisplacement(lattice, rhoBarField);
00337     Dot3D offset2 = computeRelativeDisplacement(lattice, jField);
00338 
00339     global::profiler().start("collStream");
00340     global::profiler().increment("collStreamCells", extDomain.nCells());
00341 
00342     // First, do the collision on cells within a boundary envelope of width
00343     // equal to the range of the lattice vectors (e.g. 1 for D2Q9)
00344     collide(lattice,
00345             Box3D(extDomain.x0,extDomain.x0+vicinity-1,
00346                   extDomain.y0,extDomain.y1,
00347                   extDomain.z0,extDomain.z1),
00348             rhoBarField, offset1, jField, offset2, stat);
00349     collide(lattice,
00350             Box3D(extDomain.x1-vicinity+1,extDomain.x1,
00351                   extDomain.y0,extDomain.y1,
00352                   extDomain.z0,extDomain.z1),
00353             rhoBarField, offset1, jField, offset2, stat);
00354     collide(lattice,
00355             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00356                   extDomain.y0,extDomain.y0+vicinity-1,
00357                   extDomain.z0,extDomain.z1),
00358             rhoBarField, offset1, jField, offset2, stat);
00359     collide(lattice,
00360             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00361                   extDomain.y1-vicinity+1,extDomain.y1,
00362                   extDomain.z0,extDomain.z1),
00363             rhoBarField, offset1, jField, offset2, stat);
00364     collide(lattice,
00365             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00366                   extDomain.y0+vicinity,extDomain.y1-vicinity,
00367                   extDomain.z0,extDomain.z0+vicinity-1),
00368             rhoBarField, offset1, jField, offset2, stat);
00369     collide(lattice,
00370             Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00371                   extDomain.y0+vicinity,extDomain.y1-vicinity,
00372                   extDomain.z1-vicinity+1,extDomain.z1),
00373             rhoBarField, offset1, jField, offset2, stat);
00374 
00375     // Then, do the efficient collideAndStream algorithm in the bulk,
00376     // excluding the envelope (this is efficient because there is no
00377     // if-then-else statement within the loop, given that the boundary
00378     // region is excluded)
00379     bulkCollideAndStream(lattice,
00380                          Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00381                                extDomain.y0+vicinity,extDomain.y1-vicinity,
00382                                extDomain.z0+vicinity,extDomain.z1-vicinity),
00383                          rhoBarField, offset1, jField, offset2, stat);
00384 
00385     // Finally, do streaming in the boundary envelope to conclude the
00386     // collision-stream cycle
00387     boundaryStream(lattice, extDomain, Box3D(extDomain.x0,extDomain.x0+vicinity-1,
00388                                              extDomain.y0,extDomain.y1,
00389                                              extDomain.z0,extDomain.z1) );
00390     boundaryStream(lattice, extDomain, Box3D(extDomain.x1-vicinity+1,extDomain.x1,
00391                                              extDomain.y0,extDomain.y1,
00392                                              extDomain.z0,extDomain.z1) );
00393     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00394                                              extDomain.y0,extDomain.y0+vicinity-1,
00395                                              extDomain.z0,extDomain.z1) );
00396     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00397                                              extDomain.y1-vicinity+1,extDomain.y1,
00398                                              extDomain.z0,extDomain.z1) );
00399     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00400                                              extDomain.y0+vicinity,extDomain.y1-vicinity,
00401                                              extDomain.z0,extDomain.z0+vicinity-1) );
00402     boundaryStream(lattice, extDomain, Box3D(extDomain.x0+vicinity,extDomain.x1-vicinity,
00403                                              extDomain.y0+vicinity,extDomain.y1-vicinity,
00404                                              extDomain.z1-vicinity+1,extDomain.z1) );
00405     global::profiler().stop("collStream");
00406     global::timer("collideAndStream").stop();
00407 }
00408 
00409 template<typename T, template<typename U> class Descriptor>
00410 OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>*
00411     OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>::clone() const
00412 {
00413     return new OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>(*this);
00414 }
00415 
00416 template<typename T, template<typename U> class Descriptor>
00417 void OnLinkExternalRhoJcollideAndStream3D<T,Descriptor>::getTypeOfModification (
00418         std::vector<modif::ModifT>& modified) const
00419 {
00420     modified[0] = modif::staticVariables;
00421     modified[1] = modif::nothing;
00422     modified[2] = modif::nothing;
00423 }
00424 
00425 }  // namespace plb
00426 
00427 #endif  // DYNAMICS_PROCESSOR_3D_HH