Dear All
I wrote a Smagorinsky dynamics class for MRT which is shown below. However, I am not quite sure about the way to modify the corresponding relaxation frequencies in S vector when the new omega is updated every step. I just use the " IsoThermalBulkDynamics<T,Descriptor>::setOmega(omega);" in the member function ::collide. Please look at the codes below. Any suggestions will be appreciated. Thank you very much.
/* *************** Class SmagorinskyMRTdynamics ************************************** /
template<typename T, template class Descriptor>
int SmagorinskyMRTdynamics<T,Descriptor>::id =
meta::registerGeneralDynamics<T,Descriptor,SmagorinskyMRTdynamics<T,Descriptor> >(“MRT_Smagorinsky”);
/* \param omega_ relaxation parameter, related to the dynamic viscosity
/
template<typename T, template class Descriptor>
SmagorinskyMRTdynamics<T,Descriptor>::SmagorinskyMRTdynamics (
plint externalParam_, T omega0_, T cSmago_ )
: IsoThermalBulkDynamics<T,Descriptor>(mrtParam<T,Descriptor>().get(externalParam_).getOmega()),
param(0),
externalParam(externalParam_),
omega0(omega0_),
cSmago(cSmago_),
preFactor(SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago))
{ }
template<typename T, template class Descriptor>
SmagorinskyMRTdynamics<T,Descriptor>::SmagorinskyMRTdynamics (
HierarchicUnserializer& unserializer )
: IsoThermalBulkDynamics<T,Descriptor>(T()),
param(0),
externalParam(-1),
omega0(T()),
cSmago(T()),
preFactor(T())
{
unserialize(unserializer);
}
template<typename T, template class Descriptor>
void SmagorinskyMRTdynamics<T,Descriptor>::setOmega(T omega0_)
{
omega0 = omega0_;
preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago);
}
template<typename T, template class Descriptor>
T SmagorinskyMRTdynamics<T,Descriptor>::getOmega() const {
return omega0;
}
template<typename T, template class Descriptor>
SmagorinskyMRTdynamics<T,Descriptor> SmagorinskyMRTdynamics<T,Descriptor>::clone() const {
return new SmagorinskyMRTdynamics<T,Descriptor>(this);
}
template<typename T, template class Descriptor>
int SmagorinskyMRTdynamics<T,Descriptor>::getId() const {
return id;
}
template<typename T, template class Descriptor>
void SmagorinskyMRTdynamics<T,Descriptor>::serialize(HierarchicSerializer& serializer) const
{
int useExternalParamFlag = param ? 0:1;
serializer.addValue(useExternalParamFlag);
if (param) {
param->serialize(serializer);
}
else {
serializer.addValue(externalParam);
}
serializer.addValue(omega0);
serializer.addValue(cSmago);
IsoThermalBulkDynamics<T,Descriptor>::serialize(serializer);
}
template<typename T, template class Descriptor>
void SmagorinskyMRTdynamics<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer)
{
int useExternalParamFlag = unserializer.readValue();
if (useExternalParamFlag) {
externalParam = unserializer.readValue();
param = 0;
}
else {
delete param;
param = new MRTparam<T,Descriptor>(unserializer);
externalParam = -1;
}
unserializer.readValue(omega0);
unserializer.readValue(cSmago);
preFactor = SmagoOperations<T,Descriptor>::computePrefactor(omega0,cSmago);
IsoThermalBulkDynamics<T,Descriptor>::unserialize(unserializer);
}
template<typename T, template class Descriptor>
void SmagorinskyMRTdynamics<T,Descriptor>::collide (
Cell<T,Descriptor>& cell,
BlockStatistics& statistics )
{
typedef mrtTemplates<T,Descriptor> mrtTemp;
T rhoBar;
Array<T,Descriptor::d> j;
Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
momentTemplates<T,Descriptor>::compute_rhoBar_j_PiNeq(cell, rhoBar, j, PiNeq);
T omega = SmagoOperations<T,Descriptor>::computeOmega (
omega0, preFactor, rhoBar, PiNeq );
IsoThermalBulkDynamics<T,Descriptor>::setOmega(omega);
MRTparam<T,Descriptor> parameter = param ? param : &(mrtParam<T,Descriptor>().get(externalParam));
T jSqr = mrtTemp::mrtCollision(cell, rhoBar, j, parameter->getInvM());
if (cell.takesStatistics()) {
gatherStatistics(statistics, rhoBar, jSqr * Descriptor::invRho(rhoBar) * Descriptor::invRho(rhoBar) );
}
}
template<typename T, template class Descriptor>
T SmagorinskyMRTdynamics<T,Descriptor>::computeEquilibrium(plint iPop, T rhoBar, Array<T,Descriptor::d> const& j,
T jSqr, T thetaBar) const
{
T invRho = Descriptor::invRho(rhoBar);
return dynamicsTemplates<T,Descriptor>::bgk_ma2_equilibrium(iPop, rhoBar, invRho, j, jSqr);
}