smagorinskyDynamics.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 /* Orestis Malaspinas designed some of the classes and concepts contained
00026  * in this file. */
00027 
00028 #ifndef SMAGORINSKY_DYNAMICS_HH
00029 #define SMAGORINSKY_DYNAMICS_HH
00030 
00031 #include "complexDynamics/smagorinskyDynamics.h"
00032 #include "core/util.h"
00033 #include "core/latticeStatistics.h"
00034 #include "core/dynamicsIdentifiers.h"
00035 #include "latticeBoltzmann/momentTemplates.h"
00036 #include "latticeBoltzmann/dynamicsTemplates.h"
00037 #include "latticeBoltzmann/geometricOperationTemplates.h"
00038 #include "latticeBoltzmann/externalForceTemplates.h"
00039 #include <cmath>
00040 
00041 namespace plb {
00042 
00043 template<typename T, template<typename U> class Descriptor>
00044 struct SmagoOperations {
00045     static T computePrefactor(T omega0, T cSmago) {
00046         return (T)0.5 * util::sqr(cSmago*omega0*Descriptor<T>::invCs2);
00047     }
00048     static T computeOmega(T omega0, T preFactor, T rhoBar, Array<T,SymmetricTensor<T,Descriptor>::n> const& PiNeq)
00049     {
00050         T PiNeqNormSqr = SymmetricTensor<T,Descriptor>::tensorNormSqr(PiNeq);
00051         T PiNeqNorm    = sqrt(PiNeqNormSqr);
00052         T alpha        = preFactor * Descriptor<T>::invRho(rhoBar);
00053         T linearTerm   = alpha*PiNeqNorm;
00054         T squareTerm   = (T)2*alpha*alpha*PiNeqNormSqr;
00055         // In the following formula, the square-root appearing in the explicit form of
00056         //   omega is developed to second-order.
00057         return omega0*(1-linearTerm+squareTerm);
00058     }
00059 };
00060 
00061 template<typename T, template<typename U> class Descriptor>
00062 int SmagorinskyDynamics<T,Descriptor>::id =
00063     meta::registerGeneralDynamics<T,Descriptor,SmagorinskyDynamics<T,Descriptor> >("Smagorinsky_Generic");
00064 
00065 template<typename T, template<typename U> class Descriptor>
00066 SmagorinskyDynamics<T,Descriptor>::SmagorinskyDynamics(Dynamics<T,Descriptor>* baseDynamics_,
00067                                                        T omega0_, T cSmago_, bool automaticPrepareCollision)
00068     : OmegaFromPiDynamics<T,Descriptor>(baseDynamics_, automaticPrepareCollision),
00069       omega0(omega0_),
00070       cSmago(cSmago_),
00071       preFactor(SmagoOperations<T,Descriptor>::computePrefactor(omega0, cSmago))
00072 { }
00073 
00074 template<typename T, template<typename U> class Descriptor>
00075 SmagorinskyDynamics<T,Descriptor>::SmagorinskyDynamics(HierarchicUnserializer& unserializer)
00076     : OmegaFromPiDynamics<T,Descriptor>(0, false),
00077       omega0(T()),
00078       cSmago(T()),
00079       preFactor(T())
00080 {
00081     unserialize(unserializer);
00082 }
00083 
00084 
00085 template<typename T, template<typename U> class Descriptor>
00086 T SmagorinskyDynamics<T,Descriptor>::getOmegaFromPiAndRhoBar(Array<T,SymmetricTensor<T,Descriptor>::n> const& PiNeq, T rhoBar) const
00087 {
00088     return SmagoOperations<T,Descriptor>::computeOmega(omega0, preFactor, rhoBar, PiNeq);
00089 }
00090 
00091 template<typename T, template<typename U> class Descriptor>
00092 SmagorinskyDynamics<T,Descriptor>* SmagorinskyDynamics<T,Descriptor>::clone() const {
00093     return new SmagorinskyDynamics<T,Descriptor>(*this);
00094 }
00095  
00096 template<typename T, template<typename U> class Descriptor>
00097 int SmagorinskyDynamics<T,Descriptor>::getId() const {
00098     return id;
00099 }
00100 
00101 template<typename T, template<typename U> class Descriptor>
00102 void SmagorinskyDynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00103 {
00104     serializer.addValue(omega0);
00105     serializer.addValue(cSmago);
00106     OmegaFromPiDynamics<T,Descriptor>::serialize(serializer);
00107 }
00108 
00109 template<typename T, template<typename U> class Descriptor>
00110 void SmagorinskyDynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00111 {
00112     unserializer.readValue(omega0);
00113     unserializer.readValue(cSmago);
00114     preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0, cSmago);
00115     OmegaFromPiDynamics<T,Descriptor>::unserialize(unserializer);
00116 }
00117 
00118 template<typename T, template<typename U> class Descriptor>
00119 void SmagorinskyDynamics<T,Descriptor>::setOmega(T omega0_)
00120 {
00121     // Just to be sure to avoid an undefined state, also reset the value of
00122     // omega in the baseDynamics, although the actual value of omega in the
00123     // baseDynamics is omega=omega0+deltaOmega and is recomputed before each
00124     // collision.
00125     omega0 = omega0_;
00126     preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0, cSmago);
00127     this->getBaseDynamics().setOmega(omega0);
00128 }
00129 
00130 template<typename T, template<typename U> class Descriptor>
00131 T SmagorinskyDynamics<T,Descriptor>::getOmega() const {
00132     return omega0;
00133 }
00134 
00135 
00136 /* *************** Class SmagorinskyBGKdynamics ************************************** */
00137 
00138 template<typename T, template<typename U> class Descriptor>
00139 int SmagorinskyBGKdynamics<T,Descriptor>::id =
00140     meta::registerGeneralDynamics<T,Descriptor,SmagorinskyBGKdynamics<T,Descriptor> >("BGK_Smagorinsky");
00141 
00144 template<typename T, template<typename U> class Descriptor>
00145 SmagorinskyBGKdynamics<T,Descriptor>::SmagorinskyBGKdynamics (
00146         T omega0_, T cSmago_ )
00147     : IsoThermalBulkDynamics<T,Descriptor>(omega0_),
00148       omega0(omega0_),
00149       cSmago(cSmago_),
00150       preFactor(SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago))
00151 { }
00152 
00153 template<typename T, template<typename U> class Descriptor>
00154 SmagorinskyBGKdynamics<T,Descriptor>::SmagorinskyBGKdynamics (
00155         HierarchicUnserializer& unserializer )
00156     : IsoThermalBulkDynamics<T,Descriptor>(T()),
00157       omega0(T()),
00158       cSmago(T()),
00159       preFactor(T())
00160 {
00161     unserialize(unserializer);
00162 }
00163 
00164 template<typename T, template<typename U> class Descriptor>
00165 void SmagorinskyBGKdynamics<T,Descriptor>::setOmega(T omega0_)
00166 {
00167     omega0 = omega0_;
00168     preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago);
00169 }
00170 
00171 template<typename T, template<typename U> class Descriptor>
00172 T SmagorinskyBGKdynamics<T,Descriptor>::getOmega() const {
00173     return omega0;
00174 }
00175 
00176 template<typename T, template<typename U> class Descriptor>
00177 SmagorinskyBGKdynamics<T,Descriptor>* SmagorinskyBGKdynamics<T,Descriptor>::clone() const {
00178     return new SmagorinskyBGKdynamics<T,Descriptor>(*this);
00179 }
00180  
00181 template<typename T, template<typename U> class Descriptor>
00182 int SmagorinskyBGKdynamics<T,Descriptor>::getId() const {
00183     return id;
00184 }
00185 
00186 template<typename T, template<typename U> class Descriptor>
00187 void SmagorinskyBGKdynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00188 {
00189     serializer.addValue(omega0);
00190     serializer.addValue(cSmago);
00191     IsoThermalBulkDynamics<T,Descriptor>::serialize(serializer);
00192 }
00193 
00194 template<typename T, template<typename U> class Descriptor>
00195 void SmagorinskyBGKdynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00196 {
00197     unserializer.readValue(omega0);
00198     unserializer.readValue(cSmago);
00199     preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago);
00200     IsoThermalBulkDynamics<T,Descriptor>::unserialize(unserializer);
00201 }
00202 
00203 template<typename T, template<typename U> class Descriptor>
00204 void SmagorinskyBGKdynamics<T,Descriptor>::collide (
00205         Cell<T,Descriptor>& cell,
00206         BlockStatistics& statistics )
00207 {
00208     T rhoBar;
00209     Array<T,Descriptor<T>::d> j;
00210     Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
00211     momentTemplates<T,Descriptor>::compute_rhoBar_j_PiNeq(cell, rhoBar, j, PiNeq);
00212     T omega = SmagoOperations<T,Descriptor>::computeOmega (
00213             omega0, preFactor, rhoBar, PiNeq );
00214     T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, omega);
00215     if (cell.takesStatistics()) {
00216         gatherStatistics(statistics, rhoBar, uSqr);
00217     }
00218 }
00219 
00220 template<typename T, template<typename U> class Descriptor>
00221 T SmagorinskyBGKdynamics<T,Descriptor>::computeEquilibrium (
00222         plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
00223         T jSqr, T thetaBar ) const
00224 {
00225     T invRho = Descriptor<T>::invRho(rhoBar);
00226     return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
00227 }
00228 
00229 
00230 /* *************** Class GuoExternalForceSmagorinskyBGKdynamics ************************************** */
00231 
00232 
00233 template<typename T, template<typename U> class Descriptor>
00234 int GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::id =
00235     meta::registerGeneralDynamics<T,Descriptor,GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor> >("BGK_Smagorinsky_Guo");
00236 
00239 template<typename T, template<typename U> class Descriptor>
00240 GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::GuoExternalForceSmagorinskyBGKdynamics (
00241         T omega0_, T cSmago_ )
00242     : IsoThermalBulkDynamics<T,Descriptor>(omega0_),
00243       omega0(omega0_),
00244       cSmago(cSmago_)
00245 { }
00246 
00247 template<typename T, template<typename U> class Descriptor>
00248 GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::GuoExternalForceSmagorinskyBGKdynamics(HierarchicUnserializer& unserializer)
00249     : IsoThermalBulkDynamics<T,Descriptor>(T()),
00250       omega0(T()),
00251       cSmago(T())
00252 {
00253     unserialize(unserializer);
00254 }
00255 
00256 template<typename T, template<typename U> class Descriptor>
00257 GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>* GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::clone() const {
00258     return new GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>(*this);
00259 }
00260 
00261 template<typename T, template<typename U> class Descriptor>
00262 int GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::getId() const {
00263     return id;
00264 }
00265 
00266 template<typename T, template<typename U> class Descriptor>
00267 void GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00268 {
00269     serializer.addValue(omega0);
00270     serializer.addValue(cSmago);
00271     IsoThermalBulkDynamics<T,Descriptor>::serialize(serializer);
00272 }
00273 
00274 template<typename T, template<typename U> class Descriptor>
00275 void GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00276 {
00277     unserializer.readValue(omega0);
00278     unserializer.readValue(cSmago);
00279     IsoThermalBulkDynamics<T,Descriptor>::unserialize(unserializer);
00280 }
00281 
00282 template<typename T, template<typename U> class Descriptor>
00283 void GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::setOmega(T omega0_)
00284 {
00285     omega0 = omega0_;
00286 }
00287 
00288 template<typename T, template<typename U> class Descriptor>
00289 T GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::getOmega() const {
00290     return IsoThermalBulkDynamics<T,Descriptor>::getOmega();
00291 }
00292 
00293 template<typename T, template<typename U> class Descriptor>
00294 void GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::computeVelocity (
00295                                         Cell<T,Descriptor> const& cell,
00296                                         Array<T,Descriptor<T>::d>& u ) const
00297 {
00298     T rhoBar;
00299     Array<T,Descriptor<T>::d> force, j;
00300     momentTemplates<T,Descriptor>::get_rhoBar_j(cell,rhoBar, j);
00301     force.from_cArray(cell.getExternal(Descriptor<T>::ExternalField::forceBeginsAt));
00302 
00303     T invRho = Descriptor<T>::invRho(rhoBar);
00304     for (plint iD = 0; iD < Descriptor<T>::d; ++iD)
00305     {
00306         u[iD] = j[iD]*invRho + force[iD]/(T)2;
00307     }
00308 
00309 }
00310 
00311 
00312 template<typename T, template<typename U> class Descriptor>
00313 void GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::collide (
00314         Cell<T,Descriptor>& cell,
00315         BlockStatistics& statistics )
00316 {
00317     T rhoBar = momentTemplates<T,Descriptor>::get_rhoBar(cell);
00318     Array<T,Descriptor<T>::d> u,j;
00319     computeVelocity(cell, u);
00320     T rho = Descriptor<T>::fullRho(rhoBar);
00321     j = rho * u;
00322     Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
00323     momentTemplates<T,Descriptor>::compute_PiNeq(cell, rhoBar, j, PiNeq);
00324     
00325     T tau = (T)1/omega0;
00326     
00327     Array<T,Descriptor<T>::d> force; force.from_cArray(cell.getExternal(Descriptor<T>::ExternalField::forceBeginsAt));
00328     Array<T,SymmetricTensor<T,Descriptor>::n> forceCorr; forceCorr.resetToZero();
00329     plint iPi = 0;
00330     for (plint iA = 0; iA < Descriptor<T>::d; ++iA) {
00331         for (plint iB = iA; iB < Descriptor<T>::d; ++iB) {
00332             forceCorr[iPi] = (u[iA]*force[iB]+u[iB]*force[iA]);
00333             ++iPi;
00334         }
00335     }
00336     PiNeq += forceCorr * (T)0.5;
00337     
00338     T normPiNeq = sqrt((T)2*SymmetricTensor<T,Descriptor>::tensorNormSqr(PiNeq));
00339     normPiNeq *= (T)2*rho*util::sqr(Descriptor<T>::cs2*cSmago*cSmago);
00340 
00341     T tauSGS = T();
00342     if (normPiNeq != T()) { // test to avoid division per 0
00343         Array<T,SymmetricTensor<T,Descriptor>::n> S =
00344                 (-rho*tau*Descriptor<T>::cs2+sqrt(util::sqr(rho*tau*Descriptor<T>::cs2)+normPiNeq)) / normPiNeq * PiNeq;
00345 
00346         T sNorm = sqrt((T)2*SymmetricTensor<T,Descriptor>::tensorNormSqr(S));
00347     
00348         tauSGS = cSmago*cSmago*sNorm*Descriptor<T>::invCs2;
00349     }
00350     T omega = (T)1/(tau+tauSGS);
00351 
00352     T uSqr = dynamicsTemplates<T,Descriptor>::bgk_ma2_collision(cell, rhoBar, j, omega);
00353     externalForceTemplates<T,Descriptor>::addGuoForce(cell, u, omega, (T)1);
00354     
00355     if (cell.takesStatistics()) {
00356         gatherStatistics(statistics, rhoBar, uSqr);
00357     }
00358 }
00359 
00360 template<typename T, template<typename U> class Descriptor>
00361 T GuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::computeEquilibrium (
00362         plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
00363         T jSqr, T thetaBar ) const
00364 {
00365     T invRho = Descriptor<T>::invRho(rhoBar);
00366     return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
00367 }
00368 
00369 
00370 /* *************** Class IncGuoExternalForceSmagorinskyBGKdynamics ************************************** */
00371 
00372 
00373 template<typename T, template<typename U> class Descriptor>
00374 int IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::id =
00375     meta::registerGeneralDynamics<T,Descriptor,IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor> >("IncBGK_Smagorinsky_Guo");
00376 
00379 template<typename T, template<typename U> class Descriptor>
00380 IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::IncGuoExternalForceSmagorinskyBGKdynamics (
00381         T omega0_, T cSmago_ )
00382     : IsoThermalBulkDynamics<T,Descriptor>(omega0_),
00383       omega0(omega0_),
00384       cSmago(cSmago_)
00385 { }
00386 
00387 template<typename T, template<typename U> class Descriptor>
00388 IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::IncGuoExternalForceSmagorinskyBGKdynamics(HierarchicUnserializer& unserializer)
00389     : IsoThermalBulkDynamics<T,Descriptor>(T()),
00390       omega0(T()),
00391       cSmago(T())
00392 {
00393     unserialize(unserializer);
00394 }
00395 
00396 template<typename T, template<typename U> class Descriptor>
00397 IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>* IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::clone() const {
00398     return new IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>(*this);
00399 }
00400 
00401 template<typename T, template<typename U> class Descriptor>
00402 int IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::getId() const {
00403     return id;
00404 }
00405 
00406 template<typename T, template<typename U> class Descriptor>
00407 void IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00408 {
00409     serializer.addValue(omega0);
00410     serializer.addValue(cSmago);
00411     IsoThermalBulkDynamics<T,Descriptor>::serialize(serializer);
00412 }
00413 
00414 template<typename T, template<typename U> class Descriptor>
00415 void IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00416 {
00417     unserializer.readValue(omega0);
00418     unserializer.readValue(cSmago);
00419     IsoThermalBulkDynamics<T,Descriptor>::unserialize(unserializer);
00420 }
00421 
00422 template<typename T, template<typename U> class Descriptor>
00423 void IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::setOmega(T omega0_)
00424 {
00425     omega0 = omega0_;
00426 }
00427 
00428 template<typename T, template<typename U> class Descriptor>
00429 T IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::getOmega() const {
00430     return IsoThermalBulkDynamics<T,Descriptor>::getOmega();
00431 }
00432 
00433 template<typename T, template<typename U> class Descriptor>
00434 void IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::computeVelocity (
00435                                         Cell<T,Descriptor> const& cell,
00436                                         Array<T,Descriptor<T>::d>& u ) const
00437 {
00438     T rhoBar;
00439     Array<T,Descriptor<T>::d> force, j;
00440     momentTemplates<T,Descriptor>::get_rhoBar_j(cell,rhoBar, j);
00441     force.from_cArray(cell.getExternal(Descriptor<T>::ExternalField::forceBeginsAt));
00442     
00443     T invRho = Descriptor<T>::invRho(rhoBar);
00444     for (plint iD = 0; iD < Descriptor<T>::d; ++iD)
00445     {
00446         u[iD] = j[iD]*invRho + force[iD]/(T)2;
00447     }
00448 }
00449 
00450 
00451 template<typename T, template<typename U> class Descriptor>
00452 void IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::collide (
00453         Cell<T,Descriptor>& cell,
00454         BlockStatistics& statistics )
00455 {
00456     T rhoBar = momentTemplates<T,Descriptor>::get_rhoBar(cell);
00457     Array<T,Descriptor<T>::d> u,j;
00458     computeVelocity(cell, u);
00459     T rho = Descriptor<T>::fullRho(rhoBar);
00460     j = rho * u;
00461     Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
00462     momentTemplates<T,Descriptor>::compute_PiNeq(cell, rhoBar, j, PiNeq);
00463 
00464     T tau = (T)1/omega0;
00465 
00466     Array<T,Descriptor<T>::d> force; force.from_cArray(cell.getExternal(Descriptor<T>::ExternalField::forceBeginsAt));
00467     Array<T,SymmetricTensor<T,Descriptor>::n> forceCorr; forceCorr.resetToZero();
00468     plint iPi = 0;
00469     for (plint iA = 0; iA < Descriptor<T>::d; ++iA) {
00470         for (plint iB = iA; iB < Descriptor<T>::d; ++iB) {
00471             forceCorr[iPi] = (u[iA]*force[iB]+u[iB]*force[iA]);
00472             ++iPi;
00473         }
00474     }
00475     PiNeq += forceCorr * (T)0.5;
00476 
00477     T normPiNeq = sqrt((T)2*SymmetricTensor<T,Descriptor>::tensorNormSqr(PiNeq));
00478     normPiNeq *= (T)2*rho*util::sqr(Descriptor<T>::cs2*cSmago*cSmago);
00479 
00480     T tauSGS = T();
00481     if (normPiNeq != T()) { // test to avoid division per 0
00482         Array<T,SymmetricTensor<T,Descriptor>::n> S =
00483                 (-rho*tau*Descriptor<T>::cs2+sqrt(util::sqr(rho*tau*Descriptor<T>::cs2)+normPiNeq)) / normPiNeq * PiNeq;
00484 
00485         T sNorm = sqrt((T)2*SymmetricTensor<T,Descriptor>::tensorNormSqr(S));
00486 
00487         tauSGS = cSmago*cSmago*sNorm*Descriptor<T>::invCs2;
00488     }
00489     
00490     T omega = (T)1/(tau+tauSGS);
00491 
00492     T uSqr = dynamicsTemplates<T,Descriptor>::bgk_inc_collision(cell, rhoBar, j, omega);
00493     
00494     externalForceTemplates<T,Descriptor>::addGuoForce(cell, u, omega, (T)1);
00495     if (cell.takesStatistics()) {
00496         gatherStatistics(statistics, rhoBar, uSqr);
00497     }
00498 }
00499 
00500 template<typename T, template<typename U> class Descriptor>
00501 T IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::computeEquilibrium (
00502         plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j,
00503         T jSqr, T thetaBar ) const
00504 {
00505     T invRho = (T)1;
00506     return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
00507 }
00508 
00509 template<typename T, template<typename U> class Descriptor>
00510 bool IncGuoExternalForceSmagorinskyBGKdynamics<T,Descriptor>::velIsJ() const {
00511     return true;
00512 }
00513 
00514 
00515 /* *************** Class SmagorinskyRegularizedDynamics ************************************ */
00516 
00517 template<typename T, template<typename U> class Descriptor>
00518 int SmagorinskyRegularizedDynamics<T,Descriptor>::id =
00519     meta::registerTwoParamDynamics<T,Descriptor,SmagorinskyRegularizedDynamics<T,Descriptor> >("Smagorinsky_Regularized");
00520 
00523 template<typename T, template<typename U> class Descriptor>
00524 SmagorinskyRegularizedDynamics<T,Descriptor>::SmagorinskyRegularizedDynamics (
00525         T omega0_, T cSmago_ )
00526     : IsoThermalBulkDynamics<T,Descriptor>(omega0_),
00527       omega0(omega0_),
00528       cSmago(cSmago_),
00529       preFactor(SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago))
00530 { }
00531 
00532 template<typename T, template<typename U> class Descriptor>
00533 void SmagorinskyRegularizedDynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
00534 {
00535     serializer.addValue(omega0);
00536     serializer.addValue(cSmago);
00537     IsoThermalBulkDynamics<T,Descriptor>::serialize(serializer);
00538 }
00539 
00540 
00541 template<typename T, template<typename U> class Descriptor>
00542 void SmagorinskyRegularizedDynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
00543 {
00544     unserializer.readValue(omega0);
00545     unserializer.readValue(cSmago);
00546     preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago);
00547     IsoThermalBulkDynamics<T,Descriptor>::unserialize(unserializer);
00548 }
00549 
00550 template<typename T, template<typename U> class Descriptor>
00551 void SmagorinskyRegularizedDynamics<T,Descriptor>::setOmega(T omega0_)
00552 {
00553     omega0 = omega0_;
00554     preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago);
00555 }
00556 
00557 template<typename T, template<typename U> class Descriptor>
00558 T SmagorinskyRegularizedDynamics<T,Descriptor>::getOmega() const {
00559     return omega0;
00560 }
00561 
00562 template<typename T, template<typename U> class Descriptor>
00563 SmagorinskyRegularizedDynamics<T,Descriptor>* SmagorinskyRegularizedDynamics<T,Descriptor>::clone() const
00564 {
00565     return new SmagorinskyRegularizedDynamics<T,Descriptor>(*this);
00566 }
00567  
00568 template<typename T, template<typename U> class Descriptor>
00569 int SmagorinskyRegularizedDynamics<T,Descriptor>::getId() const {
00570     return id;
00571 }
00572 
00573 template<typename T, template<typename U> class Descriptor>
00574 void SmagorinskyRegularizedDynamics<T,Descriptor>::collide (
00575         Cell<T,Descriptor>& cell,
00576         BlockStatistics& statistics )
00577 {
00578     T rhoBar;
00579     Array<T,Descriptor<T>::d> j;
00580     Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
00581     momentTemplates<T,Descriptor>::compute_rhoBar_j_PiNeq(cell, rhoBar, j, PiNeq);
00582     T omega = SmagoOperations<T,Descriptor>::computeOmega (
00583             omega0, preFactor, rhoBar, PiNeq );
00584     IsoThermalBulkDynamics<T,Descriptor>::setOmega(omega);
00585     T invRho = Descriptor<T>::invRho(rhoBar);
00586     T uSqr = dynamicsTemplates<T,Descriptor>::rlb_collision (
00587                  cell, rhoBar, invRho, j, PiNeq, omega );
00588     if (cell.takesStatistics()) {
00589         gatherStatistics(statistics, rhoBar, uSqr);
00590     }
00591 }
00592 
00593 template<typename T, template<typename U> class Descriptor>
00594 T SmagorinskyRegularizedDynamics<T,Descriptor>::computeEquilibrium (
00595         plint iPop, T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr, T thetaBar ) const
00596 {
00597     T invRho = Descriptor<T>::invRho(rhoBar);
00598     return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
00599 }
00600 
00601 } // namespace plb
00602 
00603 #endif  // SMAGORINSKY_DYNAMICS_HH