multiGridGenerator3D.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 
00025 /* Main author: Daniel Lagrava
00026  **/
00027 
00032 #ifndef MULTI_GRID_GENERATOR_3D_HH
00033 #define MULTI_GRID_GENERATOR_3D_HH
00034 
00035 #include "multiGrid/multiGridGenerator3D.h"
00036 #include "multiGrid/domainDivision3D.h"
00037 
00038 namespace plb {
00039     
00040 enum {cornerLL=0,cornerUL,cornerLR,cornerUR}; // to access the corners generated with reduce()
00041 enum {borderUp=0,borderLeft,borderDown,borderRight}; // to access the borders generated with reduce()
00042 
00044 // Dot3D const deltaBulk[3][2] = {
00045 //     {Dot3D(0,1,0), Dot3D(0,0,1)}, // YZ
00046 //     {Dot3D(1,0,0), Dot3D(0,0,1)}, // XZ
00047 //     {Dot3D(1,0,0), Dot3D(0,1,0)}  // XY
00048 // };
00049 // 
00050 // Dot3D const deltaCorner[4][3][2] = {
00051 //     {{Dot3D(0,1,0), Dot3D(0,0,1)}, {Dot3D(1,0,0),Dot3D(0,0,1)},{Dot3D(1,0,0),Dot3D(0,1,0)}}, // cornerLL
00052 //     {{Dot3D(0,1,0), Dot3D(0,0,-1)}, {Dot3D(1,0,0),Dot3D(0,0,-1)},{Dot3D(1,0,0),Dot3D(0,-1,0)}}, // cornerUL
00053 //     {{Dot3D(0,-1,0), Dot3D(0,0,1)}, {Dot3D(-1,0,0),Dot3D(0,0,1)},{Dot3D(-1,0,0),Dot3D(0,1,0)}}, // cornerLR
00054 //     {{Dot3D(0,-1,0), Dot3D(0,0,-1)}, {Dot3D(-1,0,0),Dot3D(0,0,-1)},{Dot3D(-1,0,0),Dot3D(0,-1,0)}}, // cornerUR
00055 // };
00056 // 
00057 // Dot3D const deltaBoundary[4][3][2] = {
00058 //     {{Dot3D(0,1,0), Dot3D(0,0,-1)}, {Dot3D(1,0,0),Dot3D(0,0,-1)},{Dot3D(1,0,0),Dot3D(0,-1,0)}}, // borderUp
00059 //     {{Dot3D(0,0,1), Dot3D(0,1,0)}, {Dot3D(0,0,1),Dot3D(1,0,0)},{Dot3D(0,1,0),Dot3D(1,0,0)}}, // borderLeft
00060 //     {{Dot3D(0,1,0), Dot3D(0,0,1)}, {Dot3D(1,0,0),Dot3D(0,0,1)},{Dot3D(1,0,0),Dot3D(0,1,0)}}, // borderDown
00061 //     {{Dot3D(0,0,1), Dot3D(0,-1,0)}, {Dot3D(0,0,1),Dot3D(-1,0,0)},{Dot3D(0,1,0),Dot3D(-1,0,0)}}, // borderRight
00062 // };
00063 
00064 
00065 
00066 Box3D reduce( Box3D const& fineGridInterface, std::vector<Box3D>& excess,
00067                     std::vector<Box3D>& corner, plint const& direction ){
00068     corner.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x0,
00069                    fineGridInterface.y0,fineGridInterface.y0,
00070                    fineGridInterface.z0,fineGridInterface.z0)); // corner lower left
00071     Box3D result;
00072     if (direction==0){ // plane YZ
00073         PLB_ASSERT( fineGridInterface.x0 == fineGridInterface.x1 );
00074         result = Box3D(fineGridInterface.x0,fineGridInterface.x1,
00075                      fineGridInterface.y0+1,fineGridInterface.y1-1,
00076                      fineGridInterface.z0+1,fineGridInterface.z1-1);
00077         excess.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x1,
00078                                 fineGridInterface.y0+1,fineGridInterface.y1-1,
00079                                 fineGridInterface.z1,fineGridInterface.z1)); // border up
00080         excess.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x1,
00081                                 fineGridInterface.y0,fineGridInterface.y0,
00082                                 fineGridInterface.z0+1,fineGridInterface.z1-1)); // border left
00083         excess.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x1,
00084                                 fineGridInterface.y0+1,fineGridInterface.y1-1,
00085                                 fineGridInterface.z0,fineGridInterface.z0)); // border down
00086         excess.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x1,
00087                                 fineGridInterface.y1,fineGridInterface.y1,
00088                                 fineGridInterface.z0+1,fineGridInterface.z1-1)); // border right
00089         corner.push_back(Box3D(fineGridInterface.x1,fineGridInterface.x1,
00090                    fineGridInterface.y0,fineGridInterface.y0,
00091                    fineGridInterface.z1,fineGridInterface.z1)); // corner upper left
00092         corner.push_back(Box3D(fineGridInterface.x1,fineGridInterface.x1,
00093                    fineGridInterface.y1,fineGridInterface.y1,
00094                    fineGridInterface.z0,fineGridInterface.z0)); // corner lower right 
00095     }
00096     if (direction==1){ // plane XZ
00097         PLB_ASSERT( fineGridInterface.y0 == fineGridInterface.y1 );
00098         result = Box3D(fineGridInterface.x0+1,fineGridInterface.x1-1,
00099                      fineGridInterface.y0,fineGridInterface.y1,
00100                      fineGridInterface.z0+1,fineGridInterface.z1-1);
00101         excess.push_back(Box3D(fineGridInterface.x0+1,fineGridInterface.x1-1,
00102                                 fineGridInterface.y0,fineGridInterface.y1,
00103                                 fineGridInterface.z1,fineGridInterface.z1)); // border up
00104         excess.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x0,
00105                                 fineGridInterface.y0,fineGridInterface.y1,
00106                                 fineGridInterface.z0+1,fineGridInterface.z1-1)); // border left
00107         excess.push_back(Box3D(fineGridInterface.x0+1,fineGridInterface.x1-1,
00108                                 fineGridInterface.y0,fineGridInterface.y1,
00109                                 fineGridInterface.z0,fineGridInterface.z0)); // border down
00110         excess.push_back(Box3D(fineGridInterface.x1,fineGridInterface.x1,
00111                                 fineGridInterface.y0,fineGridInterface.y1,
00112                                 fineGridInterface.z0+1,fineGridInterface.z1-1)); // border right
00113         corner.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x0,
00114                    fineGridInterface.y0,fineGridInterface.y1,
00115                    fineGridInterface.z1,fineGridInterface.z1)); // corner upper left
00116         corner.push_back(Box3D(fineGridInterface.x1,fineGridInterface.x1,
00117                    fineGridInterface.y1,fineGridInterface.y1,
00118                    fineGridInterface.z0,fineGridInterface.z0)); // corner lower right
00119     }
00120     if (direction==2){ // plane XY
00121         PLB_ASSERT( fineGridInterface.z0 == fineGridInterface.z1 );
00122         result = Box3D(fineGridInterface.x0+1,fineGridInterface.x1-1,
00123                      fineGridInterface.y0+1,fineGridInterface.y1-1,
00124                      fineGridInterface.z0,fineGridInterface.z1);
00125         excess.push_back(Box3D(fineGridInterface.x0+1,fineGridInterface.x1-1,
00126                                 fineGridInterface.y1,fineGridInterface.y1,
00127                                 fineGridInterface.z0,fineGridInterface.z1)); // border up
00128         excess.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x0,
00129                                 fineGridInterface.y0+1,fineGridInterface.y1-1,
00130                                 fineGridInterface.z0,fineGridInterface.z1)); // border left
00131         excess.push_back(Box3D(fineGridInterface.x0+1,fineGridInterface.x1-1,
00132                                 fineGridInterface.y0,fineGridInterface.y0,
00133                                 fineGridInterface.z0,fineGridInterface.z1)); // border down
00134         excess.push_back(Box3D(fineGridInterface.x1,fineGridInterface.x1,
00135                                 fineGridInterface.y0+1,fineGridInterface.y1-1,
00136                                 fineGridInterface.z0,fineGridInterface.z1)); // border right
00137         
00138         corner.push_back(Box3D(fineGridInterface.x0,fineGridInterface.x0,
00139                    fineGridInterface.y1,fineGridInterface.y1,
00140                    fineGridInterface.z1,fineGridInterface.z1)); // corner upper left
00141         corner.push_back(Box3D(fineGridInterface.x1,fineGridInterface.x1,
00142                    fineGridInterface.y0,fineGridInterface.y0,
00143                    fineGridInterface.z1,fineGridInterface.z1)); // corner lower right
00144     }
00145 
00146     corner.push_back(Box3D(fineGridInterface.x1,fineGridInterface.x1,
00147                    fineGridInterface.y1,fineGridInterface.y1,
00148                    fineGridInterface.z1,fineGridInterface.z1)); // corner upper right
00149 
00150     return result;
00151 }
00152 
00154 template<typename T, template<typename U> class Descriptor>
00155 std::vector<MultiBlockLattice3D<T,Descriptor>*> generateLattices(
00156                 MultiGridManagement3D management,
00157                 std::vector<Dynamics<T,Descriptor>*> backgroundDynamics,
00158                 std::vector<BlockCommunicator3D*> communicators,
00159                 std::vector<CombinedStatistics*> combinedStatistics )
00160 {
00161     // get the MPI id's (if available)
00162     std::vector<std::vector<plint> > mpiProcesses = management.getMpiProcesses();
00163     std::vector<std::vector<Box3D> > bulks = management.getBulks();
00164     std::vector<std::vector<Box3D> > uniqueBulks = management.getUniqueBulks();
00165     PLB_PRECONDITION (bulks.size() == uniqueBulks.size());
00166     PLB_PRECONDITION (mpiProcesses.size()==0 || mpiProcesses.size()==bulks.size());
00167     // for cubic interpolation 
00168     plint envelopeWidth = 2; 
00169     // allocation of the result vector
00170     std::vector<MultiBlockLattice3D<T,Descriptor>*> multiBlocks(bulks.size());
00171     for (pluint iLevel=0; iLevel<bulks.size(); ++iLevel) {
00172         // create a thread attribution 
00173         ExplicitThreadAttribution* threadAttribution = new ExplicitThreadAttribution;
00174         PLB_ASSERT(mpiProcesses.size()==0 || mpiProcesses[iLevel].size()==bulks[iLevel].size());
00175         SparseBlockStructure3D geometry(management.getBoundingBox(iLevel));
00176         for (pluint iBlock=0; iBlock<bulks[iLevel].size(); ++iBlock) {
00177             plint mpiProcess = mpiProcesses.size() > 0 ? mpiProcesses[iLevel][iBlock] : 0;
00178             plint blockId = geometry.nextIncrementalId();
00179             Box3D bulk = bulks[iLevel][iBlock];
00180             Box3D uniqueBulk = bulk;
00181             geometry.addBlock(bulk, uniqueBulk, blockId); //TODO: add uniqueBulk
00182             threadAttribution->addBlock(blockId, mpiProcess);
00183         }
00184         MultiBlockManagement3D
00185             blockManagement( geometry, threadAttribution,
00186                              envelopeWidth, iLevel );
00187         multiBlocks[iLevel] = new MultiBlockLattice3D<T,Descriptor> (
00188             blockManagement,
00189             communicators[iLevel],
00190             combinedStatistics[iLevel],
00191             defaultMultiBlockPolicy3D().getMultiCellAccess<T,Descriptor>(),
00192             backgroundDynamics[iLevel] );
00193         multiBlocks[iLevel] -> periodicity().toggleAll(false);
00194     }
00195     return multiBlocks;
00196 }
00197 
00198 template<typename T, template<typename U> class Descriptor>
00199 std::vector<MultiBlockLattice3D<T,Descriptor>*> generateLattices(
00200                 MultiGridManagement3D management,
00201                 std::vector<Dynamics<T,Descriptor>*> backgroundDynamics)
00202 {
00203     return generateLattices 
00204         (management, backgroundDynamics,
00205          defaultMultiGridPolicy3D().getBlockCommunicator(management.getNumLevels()),
00206          defaultMultiGridPolicy3D().getCombinedStatistics(management.getNumLevels()));
00207 }
00208 
00209 
00210 template<typename T, template<typename U> class Descriptor>
00211 void createInterfaces(std::vector<MultiBlockLattice3D<T,Descriptor>*>& multiBlocks,
00212                       MultiGridManagement3D management)
00213 {
00214     std::vector<std::vector<Box3D> > fineGridInterfaces = management.getFineInterface();
00215     std::vector<std::vector<Box3D> > coarseGridInterfaces = management.getCoarseInterface();
00216     
00217     // Create the part of the coupling in which a lattice acts as a coarse lattice.
00218     for (pluint iCoarse=0; iCoarse<multiBlocks.size()-1; ++iCoarse) {
00219         for (pluint iInterface=0; iInterface<coarseGridInterfaces[iCoarse].size(); ++iInterface) {
00220             createCoarseGridInterface( iCoarse, coarseGridInterfaces[iCoarse][iInterface], multiBlocks );
00221         }
00222     }
00223     // Create the part of the coupling in which a lattice acts as a fine lattice.
00224     for (pluint iFine=1; iFine<multiBlocks.size(); ++iFine) {
00225         for (pluint iInterface=0; iInterface<fineGridInterfaces[iFine].size(); ++iInterface) {
00226             createFineGridInterface( iFine-1, fineGridInterfaces[iFine][iInterface], multiBlocks );
00227         }
00228     }
00229 }
00230 
00231 
00232 template<typename T, template<typename U> class Descriptor>
00233 void createCoarseGridInterface (
00234         plint coarseLevel, Box3D coarseGridInterface, std::vector<MultiBlockLattice3D<T,Descriptor>*>& multiBlocks )
00235 {
00236     PLB_PRECONDITION( coarseGridInterface.getNx()==1 ||
00237                       coarseGridInterface.getNy()==1 ||
00238                       coarseGridInterface.getNz()==1 );
00239                       
00240     pluint fineLevel = coarseLevel + 1;
00241     PLB_PRECONDITION( fineLevel < multiBlocks.size() );
00242 
00243     MultiBlockLattice3D<T,Descriptor>& coarseLattice = *multiBlocks[coarseLevel];
00244     MultiBlockLattice3D<T,Descriptor>& fineLattice   = *multiBlocks[fineLevel];
00245 
00246     // Number of time steps executed by the fine grid during a coarse iteration.
00247     plint numTimeSteps = 2; 
00248     plint executionTime = numTimeSteps-1;
00249     // Add a data processor which copies and rescales from the fine to the coarse lattice.
00250     plint scalingOrder=1;
00251     integrateProcessingFunctional (
00252             new CopyFineToCoarse3D<T, Descriptor,Descriptor> (
00253                 new ConvectiveRescaleEngine<T,Descriptor>(scalingOrder), numTimeSteps, executionTime),
00254             coarseGridInterface.multiply(2),
00255             fineLattice, coarseLattice );
00256 }
00257 
00258 
00259 template<typename T, template<typename U> class Descriptor>
00260 void createFineGridInterface (
00261         plint coarseLevel, Box3D fineGridInterface, std::vector<MultiBlockLattice3D<T,Descriptor>*>& multiBlocks )
00262 {
00263     plint scalingOrder=1;
00264     RescaleEngine<T,Descriptor>* rescaleEngine = new ConvectiveRescaleEngine<T,Descriptor>(scalingOrder);
00265             
00266     PLB_PRECONDITION( fineGridInterface.getNx()==1 || fineGridInterface.getNy()==1 || fineGridInterface.getNz()==1);
00267     pluint fineLevel = coarseLevel + 1;
00268     PLB_PRECONDITION( fineLevel < multiBlocks.size() );
00269 
00270     MultiBlockLattice3D<T,Descriptor>& coarseLattice = *multiBlocks[coarseLevel];
00271     MultiBlockLattice3D<T,Descriptor>& fineLattice   = *multiBlocks[fineLevel];
00272 
00273     plint direction;
00274     if (fineGridInterface.getNx()==1) direction = 0;
00275     if (fineGridInterface.getNy()==1) direction = 1;
00276     if (fineGridInterface.getNz()==1) direction = 2;
00277 
00278     // Number of time steps executed by the fine grid during a coarse iteration.
00279     plint numTimeSteps = 2; 
00280     plint executionTime = numTimeSteps-1;
00281 
00282     // Instantiate specific time-interpolation dynamics on the interface of the fine grid.
00283     setCompositeDynamics (
00284             fineLattice,
00285             fineGridInterface.multiply(2),
00286             new FineGridBoundaryDynamics<T,Descriptor> (
00287                 new NoDynamics<T,Descriptor>, fineLattice.getTimeCounter(), numTimeSteps,
00288                 rescaleEngine->getDecompositionOrder() ) );
00289 
00290     
00291     // The coarse processors are at processing level 1, because they need to access information
00292     //   from within the coarse envelope.
00293     plint processorLevelCoarse=1;
00294    
00295     //TODO add interpolation processors
00296     
00297     // Add a data processor which imposes a time-cyclic behavior on the fine grid
00298     //   boundary-dynamics.
00299     integrateProcessingFunctional (
00300             new Copy_t1_to_t0_3D<T, Descriptor>(numTimeSteps, executionTime),
00301             fineGridInterface.multiply(2),
00302             fineLattice );
00303 
00304     delete rescaleEngine;
00305 }
00306 
00307 
00308 
00309 template<typename T>
00310 std::vector<MultiScalarField3D<T>*> generateScalarFields(
00311                         MultiGridManagement3D const& management,
00312                         std::vector<BlockCommunicator3D*> communicators,
00313                         std::vector<CombinedStatistics*> combinedStatistics )
00314 {
00315     std::vector<std::vector<Box3D> > bulks = management.getBulks(); 
00316     std::vector<std::vector<Box3D> > uniqueBulks = management.getUniqueBulks(); 
00317     std::vector<std::vector<plint> > mpiProcesses = management.getMpiProcesses();
00318     PLB_PRECONDITION (bulks.size() == uniqueBulks.size());
00319     PLB_PRECONDITION (mpiProcesses.size()==0 || mpiProcesses.size()==bulks.size());
00320     plint envelopeWidth=1;
00321     std::vector<MultiScalarField3D<T>*> scalars(bulks.size());
00322     for (plint iLevel=0; iLevel<(plint)bulks.size(); ++iLevel){
00323         PLB_ASSERT(mpiProcesses.size()==0 || mpiProcesses[iLevel].size()==bulks[iLevel].size());
00324         ExplicitThreadAttribution* threadAttribution = new ExplicitThreadAttribution;
00325         SparseBlockStructure3D geometry(management.getBoundingBox(iLevel));
00326         for (pluint iBlock=0; iBlock<bulks[iLevel].size(); ++iBlock) {
00327             plint mpiProcess = mpiProcesses.size() > 0 ? mpiProcesses[iLevel][iBlock] : 0;
00328             plint blockId = geometry.nextIncrementalId();
00329             geometry.addBlock(bulks[iLevel][iBlock], bulks[iLevel][iBlock], blockId);
00330             threadAttribution->addBlock(blockId, mpiProcess);
00331         }
00332         MultiBlockManagement3D blockManagement (
00333                 geometry, threadAttribution,
00334                 envelopeWidth, iLevel );
00335         scalars[iLevel] = new MultiScalarField3D<T>(
00336                                 blockManagement,communicators[iLevel],
00337                                 combinedStatistics[iLevel], 
00338                                 defaultMultiBlockPolicy3D().getMultiScalarAccess<T>());  
00339         
00340     }
00341     
00342     return scalars;
00343 }
00344 
00345 
00346 template<typename T, int nDim>
00347 std::vector<MultiTensorField3D<T,nDim> *> generateTensorFields (
00348         MultiGridManagement3D const& management,
00349         std::vector<BlockCommunicator3D*> communicators,
00350         std::vector<CombinedStatistics*> combinedStatistics )
00351 {
00352     std::vector<MultiTensorField3D<T,nDim>* > fields(management.getNumLevels());
00353     std::vector<std::vector<Box3D> > bulks = management.getBulks(); 
00354     std::vector<std::vector<plint> > ids = management.getMpiProcesses();
00355     PLB_PRECONDITION (ids.size()==0 || ids.size()==bulks.size());
00356     plint envelopeWidth=1;
00357     fields.resize(bulks.size());
00358     for (plint iLevel=0; iLevel<(plint)bulks.size(); ++iLevel){
00359         ExplicitThreadAttribution* threadAttribution = new ExplicitThreadAttribution;
00360         PLB_ASSERT(ids.size()==0 || ids[iLevel].size()==bulks[iLevel].size());
00361         SparseBlockStructure3D geometry(management.getBoundingBox(iLevel));
00362         for (pluint iBlock=0; iBlock<bulks[iLevel].size(); ++iBlock) {
00363             plint mpiProcess = ids.size() > 0 ? ids[iLevel][iBlock] : 0;
00364             plint blockId = geometry.nextIncrementalId();
00365             geometry.addBlock(bulks[iLevel][iBlock], blockId);
00366             threadAttribution->addBlock(blockId, mpiProcess);
00367         }
00368        
00369         MultiBlockManagement3D blockManagement (
00370                 geometry, threadAttribution,
00371                 envelopeWidth, iLevel );
00372         fields[iLevel] = new MultiTensorField3D<T,nDim>(
00373                                     blockManagement,   
00374                                     communicators[iLevel],
00375                                     combinedStatistics[iLevel], 
00376                                     defaultMultiBlockPolicy3D().getMultiTensorAccess<T,nDim>());  
00377         
00378     }
00379     
00380     return fields;
00381 }
00382 
00383 
00384 
00385 } // namespace plb
00386 
00387 #endif  // MULTI_GRID_GENERATOR_3D_HH
00388