← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/blender/intrepid-updates

« back to all changes in this revision

Viewing changes to intern/elbeem/intern/solver_adap.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Cyril Brulebois
  • Date: 2008-08-08 02:45:40 UTC
  • mfrom: (12.1.14 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080808024540-kkjp7ekfivzhuw3l
Tags: 2.46+dfsg-4
http://bugs.debian.org/492280
* Fix python syntax warning in import_dxf.py, which led to nasty output
  in installation/upgrade logs during byte-compilation, using a patch
  provided by the script author (Closes: #492280):
   - debian/patches/45_fix_python_syntax_warning

Show diffs side-by-side

added added

removed removed

Lines of Context:
intern/elbeem/intern/solver_adap.cpp
1
1
/******************************************************************************
2
2
 *
3
3
 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4
 
 * Copyright 2003,2004,2005 Nils Thuerey
 
4
 * Copyright 2003-2006 Nils Thuerey
5
5
 *
6
6
 * Adaptivity functions
7
7
 *
11
11
#include "solver_relax.h"
12
12
#include "particletracer.h"
13
13
 
14
 
 
 
14
#if (defined (__sun__) || defined (__sun)) || (!defined(linux) && (defined (__sparc) || defined (__sparc__)))
 
15
#include <ieeefp.h>
 
16
#endif
15
17
 
16
18
/*****************************************************************************/
17
19
//! coarse step functions
21
23
 
22
24
void LbmFsgrSolver::coarseCalculateFluxareas(int lev)
23
25
{
24
 
#if LBM_NOADCOARSENING==1
25
 
        if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
26
 
        lev =0; // get rid of warnings...
27
 
#else
28
26
        FSGR_FORIJK_BOUNDS(lev) {
29
27
                if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
30
28
                        if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
50
48
                }
51
49
        } // } TEST DEBUG
52
50
        if(!this->mSilent){ debMsgStd("coarseCalculateFluxareas",DM_MSG,"level "<<lev<<" calculated", 7); }
53
 
#endif //! LBM_NOADCOARSENING==1
54
51
}
55
52
        
56
53
void LbmFsgrSolver::coarseAdvance(int lev)
57
54
{
58
 
#if LBM_NOADCOARSENING==1
59
 
        if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
60
 
        lev =0; // get rid of warnings...
61
 
#else
62
55
        LbmFloat calcCurrentMass = 0.0;
63
56
        LbmFloat calcCurrentVolume = 0.0;
64
57
 
177
170
  mLevel[lev].lmass   = calcCurrentMass   * mLevel[lev].lcellfactor;
178
171
  mLevel[lev].lvolume = calcCurrentVolume * mLevel[lev].lcellfactor;
179
172
#if ELBEEM_PLUGIN!=1
180
 
  errMsg("DFINI", " m l"<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<"  lcf="<< mLevel[lev].lcellfactor );
181
 
  errMsg("DFINI", " v l"<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<"  lcf="<< mLevel[lev].lcellfactor );
 
173
  debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "mass: lev="<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<"  lcf="<< mLevel[lev].lcellfactor, 8 );
 
174
  debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "volume: lev="<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<"  lcf="<< mLevel[lev].lcellfactor, 8 );
182
175
#endif // ELBEEM_PLUGIN
183
 
#endif //! LBM_NOADCOARSENING==1
184
176
}
185
177
 
186
178
/*****************************************************************************/
191
183
// get dfs from level (lev+1) to (lev) coarse border nodes
192
184
void LbmFsgrSolver::coarseRestrictFromFine(int lev)
193
185
{
194
 
#if LBM_NOADCOARSENING==1
195
 
        if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
196
 
        lev =0; // get rid of warnings...
197
 
#else
198
186
        if((lev<0) || ((lev+1)>mMaxRefine)) return;
199
187
#if FSGR_STRICT_DEBUG==1
200
188
        // reset all unused cell values to invalid
245
233
                } // & fluid
246
234
        }}}
247
235
        if(!this->mSilent){ errMsg("coarseRestrictFromFine"," from l"<<(lev+1)<<",s"<<mLevel[lev+1].setCurr<<" to l"<<lev<<",s"<<mLevel[lev].setCurr); }
248
 
#endif //! LBM_NOADCOARSENING==1
249
236
}
250
237
 
251
238
bool LbmFsgrSolver::adaptGrid(int lev) {
252
 
#if LBM_NOADCOARSENING==1
253
 
        if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
254
 
        lev =0; // get rid of warnings...
255
 
        return false;
256
 
#else
257
239
        if((lev<0) || ((lev+1)>mMaxRefine)) return false;
258
240
        bool change = false;
259
241
        { // refinement, PASS 1-3
706
688
 
707
689
        if(!this->mSilent){ errMsg("adaptGrid"," for l"<<lev<<" done " ); }
708
690
        return change;
709
 
#endif //! LBM_NOADCOARSENING==1
710
691
}
711
692
 
712
693
/*****************************************************************************/
715
696
 
716
697
void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet)
717
698
{
718
 
#if LBM_NOADCOARSENING==1
719
 
        if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
720
 
        i=j=k=srcSet=dstSet=lev =0; // get rid of warnings...
721
 
#else
722
699
        LbmFloat *ccel = RACPNT(lev+1, 2*i,2*j,2*k,srcSet);
723
700
        LbmFloat *tcel = RACPNT(lev  , i,j,k      ,dstSet);
724
701
 
862
839
        RAC(tcel, dWT) = (lcsmeq[dWT] + (MSRC_WT-lcsmeq[dWT] )*lcsmdfscale);
863
840
        RAC(tcel, dWB) = (lcsmeq[dWB] + (MSRC_WB-lcsmeq[dWB] )*lcsmdfscale);
864
841
#                               endif // OPT3D==0
865
 
#endif //! LBM_NOADCOARSENING==1
866
842
}
867
843
 
868
844
void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet, bool markNbs) {
869
 
#if LBM_NOADCOARSENING==1
870
 
        if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
871
 
        i=j=k=dstSet=lev =0; // get rid of warnings...
872
 
        t=0.0; flagSet=0; markNbs=false;
873
 
#else
874
845
        LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;
875
846
        LbmFloat intDf[19] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
876
847
 
952
923
 
953
924
        IDF_WRITEBACK;
954
925
        return;
955
 
#endif //! LBM_NOADCOARSENING==1
956
926
}
957
927
 
958
928
 
959
929
 
 
930
// required globals
 
931
extern bool glob_mpactive;
 
932
extern int glob_mpnum, glob_mpindex;
 
933
#define MPTADAP_INTERV 4
 
934
 
960
935
/*****************************************************************************/
961
936
/*! change the  size of the LBM time step */
962
937
/*****************************************************************************/
966
941
 
967
942
        bool rescale = false;  // do any rescale at all?
968
943
        LbmFloat scaleFac = -1.0; // timestep scaling
969
 
        if(this->mPanic) return;
 
944
        if(mPanic) return;
970
945
 
971
946
        LbmFloat levOldOmega[FSGR_MAXNOOFLEVELS];
972
947
        LbmFloat levOldStepsize[FSGR_MAXNOOFLEVELS];
974
949
                levOldOmega[lev] = mLevel[lev].omega;
975
950
                levOldStepsize[lev] = mLevel[lev].timestep;
976
951
        }
977
 
        //if(mTimeSwitchCounts>0){ errMsg("DEB CSKIP",""); return; } // DEBUG
978
952
 
979
953
        const LbmFloat reduceFac = 0.8;          // modify time step by 20%, TODO? do multiple times for large changes?
980
954
        LbmFloat diffPercent = 0.05; // dont scale if less than 5%
981
 
        LbmFloat allowMax = this->mpParam->getTadapMaxSpeed();  // maximum allowed velocity
982
 
        LbmFloat nextmax = this->mpParam->getSimulationMaxSpeed() + norm(mLevel[mMaxRefine].gravity);
983
 
 
984
 
        //newdt = this->mpParam->getTimestep() * (allowMax/nextmax);
985
 
        LbmFloat newdt = this->mpParam->getTimestep(); // newtr
 
955
        LbmFloat allowMax = mpParam->getTadapMaxSpeed();  // maximum allowed velocity
 
956
        LbmFloat nextmax = mpParam->getSimulationMaxSpeed() + norm(mLevel[mMaxRefine].gravity);
 
957
 
 
958
        // sync nextmax
 
959
#if LBM_INCLUDE_TESTSOLVERS==1
 
960
        if(glob_mpactive) {
 
961
                if(mLevel[mMaxRefine].lsteps % MPTADAP_INTERV != MPTADAP_INTERV-1) {
 
962
                        debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "mpact:"<<glob_mpactive<<","<<glob_mpindex<<"/"<<glob_mpnum<<" step:"<<mLevel[mMaxRefine].lsteps<<" skipping tadap...",8);
 
963
                        return;
 
964
                }
 
965
                nextmax = mrInitTadap(nextmax);
 
966
        }
 
967
#endif // LBM_INCLUDE_TESTSOLVERS==1
 
968
 
 
969
        LbmFloat newdt = mpParam->getTimestep(); // newtr
986
970
        if(nextmax > allowMax/reduceFac) {
987
971
                mTimeMaxvelStepCnt++; }
988
972
        else { mTimeMaxvelStepCnt=0; }
989
973
        
990
974
        // emergency, or 10 steps with high vel
991
975
        if((mTimeMaxvelStepCnt>5) || (nextmax> (1.0/3.0)) || (mForceTimeStepReduce) ) {
992
 
        //if(nextmax > allowMax/reduceFac) {
993
 
                newdt = this->mpParam->getTimestep() * reduceFac;
 
976
                newdt = mpParam->getTimestep() * reduceFac;
994
977
        } else {
995
978
                if(nextmax<allowMax*reduceFac) {
996
 
                        newdt = this->mpParam->getTimestep() / reduceFac;
 
979
                        newdt = mpParam->getTimestep() / reduceFac;
997
980
                }
998
981
        } // newtr
999
 
        //errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<reduceFac<<" simmaxv="<< this->mpParam->getSimulationMaxSpeed() );
 
982
        //errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<reduceFac<<" simmaxv="<< mpParam->getSimulationMaxSpeed() );
1000
983
 
1001
984
        bool minCutoff = false;
1002
985
        LbmFloat desireddt = newdt;
1003
 
        if(newdt>this->mpParam->getMaxTimestep()){ newdt = this->mpParam->getMaxTimestep(); }
1004
 
        if(newdt<this->mpParam->getMinTimestep()){ 
1005
 
                newdt = this->mpParam->getMinTimestep(); 
 
986
        if(newdt>mpParam->getMaxTimestep()){ newdt = mpParam->getMaxTimestep(); }
 
987
        if(newdt<mpParam->getMinTimestep()){ 
 
988
                newdt = mpParam->getMinTimestep(); 
1006
989
                if(nextmax>allowMax/reduceFac){ minCutoff=true; } // only if really large vels...
1007
990
        }
1008
991
 
1009
 
        LbmFloat dtdiff = fabs(newdt - this->mpParam->getTimestep());
1010
 
        if(!this->mSilent) {
1011
 
                debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt<<" max"<<this->mpParam->getMaxTimestep()<<" min"<<this->mpParam->getMinTimestep()<<" diff"<<dtdiff<<
1012
 
                        " simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep) , 10); }
 
992
        LbmFloat dtdiff = fabs(newdt - mpParam->getTimestep());
 
993
        if(!mSilent) {
 
994
                debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt
 
995
                        <<" max"<<mpParam->getMaxTimestep()<<" min"<<mpParam->getMinTimestep()<<" diff"<<dtdiff
 
996
                        <<" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep)<<
 
997
                        " olddt"<<levOldStepsize[mMaxRefine]<<" redlock"<<mTimestepReduceLock 
 
998
                        , 10); }
1013
999
 
1014
1000
        // in range, and more than X% change?
1015
 
        //if( newdt <  this->mpParam->getTimestep() ) // DEBUG
 
1001
        //if( newdt <  mpParam->getTimestep() ) // DEBUG
1016
1002
        LbmFloat rhoAvg = mCurrentMass/mCurrentVolume;
1017
 
        if( (newdt<=this->mpParam->getMaxTimestep()) && (newdt>=this->mpParam->getMinTimestep()) 
1018
 
                        && (dtdiff>(this->mpParam->getTimestep()*diffPercent)) ) {
 
1003
        if( (newdt<=mpParam->getMaxTimestep()) && (newdt>=mpParam->getMinTimestep()) 
 
1004
                        && (dtdiff>(mpParam->getTimestep()*diffPercent)) ) {
1019
1005
                if((newdt>levOldStepsize[mMaxRefine])&&(mTimestepReduceLock)) {
1020
1006
                        // wait some more...
1021
1007
                        //debMsgNnl("LbmFsgrSolver::TAdp",DM_NOTIFY," Delayed... "<<mTimestepReduceLock<<" ",10);
1022
1008
                        debMsgDirect("D");
1023
1009
                } else {
1024
 
                        this->mpParam->setDesiredTimestep( newdt );
 
1010
                        mpParam->setDesiredTimestep( newdt );
1025
1011
                        rescale = true;
1026
 
                        if(!this->mSilent) {
 
1012
                        if(!mSilent) {
1027
1013
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"\n\n\n\n",10);
1028
 
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: new="<<newdt<<" old="<<this->mpParam->getTimestep()<<" maxSpeed:"<<this->mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<this->mStepCnt, 10 );
1029
 
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: "<<
1030
 
                                                "rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
 
1014
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: new="<<newdt<<" old="<<mpParam->getTimestep()
 
1015
                                                <<" maxSpeed:"<<mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<mStepCnt, 10 );
 
1016
                                //debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: "<< "rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
1031
1017
                        }
1032
1018
                } // really change dt
1033
1019
        }
1039
1025
        /*const int tadtogInter = 100;
1040
1026
        const double tadtogSwitch = 0.66;
1041
1027
        errMsg("TIMESWITCHTOGGLETEST","warning enabled "<< tadtogSwitch<<","<<tadtogSwitch<<" !!!!!!!!!!!!!!!!!!!");
1042
 
        if( ((this->mStepCnt% tadtogInter)== (tadtogInter/4*1)-1) ||
1043
 
            ((this->mStepCnt% tadtogInter)== (tadtogInter/4*2)-1) ){
 
1028
        if( ((mStepCnt% tadtogInter)== (tadtogInter/4*1)-1) ||
 
1029
            ((mStepCnt% tadtogInter)== (tadtogInter/4*2)-1) ){
1044
1030
                rescale = true; minCutoff = false;
1045
 
                newdt = tadtogSwitch * this->mpParam->getTimestep();
1046
 
                this->mpParam->setDesiredTimestep( newdt );
 
1031
                newdt = tadtogSwitch * mpParam->getTimestep();
 
1032
                mpParam->setDesiredTimestep( newdt );
1047
1033
        } else 
1048
 
        if( ((this->mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) ||
1049
 
            ((this->mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){
 
1034
        if( ((mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) ||
 
1035
            ((mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){
1050
1036
                rescale = true; minCutoff = false;
1051
 
                newdt = this->mpParam->getTimestep()/tadtogSwitch ;
1052
 
                this->mpParam->setDesiredTimestep( newdt );
 
1037
                newdt = mpParam->getTimestep()/tadtogSwitch ;
 
1038
                mpParam->setDesiredTimestep( newdt );
1053
1039
        } else {
1054
1040
                rescale = false; minCutoff = false;
1055
1041
        }
1057
1043
 
1058
1044
        // test mass rescale
1059
1045
 
1060
 
        scaleFac = newdt/this->mpParam->getTimestep();
 
1046
        scaleFac = newdt/mpParam->getTimestep();
1061
1047
        if(rescale) {
1062
1048
                // perform acutal rescaling...
1063
1049
                mTimeMaxvelStepCnt=0; 
1064
1050
                mForceTimeStepReduce = false;
1065
1051
 
1066
1052
                // FIXME - approximate by averaging, take gravity direction here?
1067
 
                mTimestepReduceLock = 4*(mLevel[mMaxRefine].lSizey+mLevel[mMaxRefine].lSizez+mLevel[mMaxRefine].lSizex)/3;
 
1053
                //mTimestepReduceLock = 4*(mLevel[mMaxRefine].lSizey+mLevel[mMaxRefine].lSizez+mLevel[mMaxRefine].lSizex)/3;
 
1054
                // use z as gravity direction
 
1055
                mTimestepReduceLock = 4*mLevel[mMaxRefine].lSizez;
1068
1056
 
1069
1057
                mTimeSwitchCounts++;
1070
 
                this->mpParam->calculateAllMissingValues( mSimulationTime, this->mSilent );
 
1058
                mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
1071
1059
                recalculateObjectSpeeds();
1072
1060
                // calc omega, force for all levels
1073
1061
                mLastOmega=1e10; mLastGravity=1e10;
1074
1062
                initLevelOmegas();
1075
 
                if(this->mpParam->getTimestep()<mMinTimestep) mMinTimestep = this->mpParam->getTimestep();
1076
 
                if(this->mpParam->getTimestep()>mMaxTimestep) mMaxTimestep = this->mpParam->getTimestep();
 
1063
                if(mpParam->getTimestep()<mMinTimestep) mMinTimestep = mpParam->getTimestep();
 
1064
                if(mpParam->getTimestep()>mMaxTimestep) mMaxTimestep = mpParam->getTimestep();
1077
1065
 
1078
1066
                // this might be called during init, before we have any particles
1079
1067
                if(mpParticles) { mpParticles->adaptPartTimestep(scaleFac); }
1088
1076
                        LbmFloat newSteptime = mLevel[lev].timestep;
1089
1077
                        LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]);
1090
1078
 
1091
 
                        if(!this->mSilent) {
1092
 
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Level: "<<lev<<" Timestep change: "<<
 
1079
                        if(!mSilent) {
 
1080
                                debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Level: "<<lev<<" Timestep chlevel: "<<
1093
1081
                                                " scaleFac="<<dfScaleFac<<" newDt="<<newSteptime<<" newOmega="<<mLevel[lev].omega,10);
1094
1082
                        }
1095
1083
                        if(lev!=mMaxRefine) coarseCalculateFluxareas(lev);
1109
1097
                                                        // init empty zero vel  interface cell...
1110
1098
                                                        initVelocityCell(lev, i,j,k, CFInter, 1.0, 0.01, LbmVec(0.) );
1111
1099
                                                } else {// */
1112
 
                                                        for(int l=0; l<this->cDfNum; l++) {
 
1100
                                                        for(int l=0; l<cDfNum; l++) {
1113
1101
                                                                QCELL(lev, i, j, k, workSet, l) = QCELL(lev, i, j, k, workSet, l)* scaleFac; 
1114
1102
                                                        }
1115
1103
                                                //} //  ok
1132
1120
                                        LbmVec velOld;
1133
1121
                                        LbmFloat rho, ux,uy,uz;
1134
1122
                                        rho=0.0; ux =  uy = uz = 0.0;
1135
 
                                        for(int l=0; l<this->cDfNum; l++) {
 
1123
                                        for(int l=0; l<cDfNum; l++) {
1136
1124
                                                LbmFloat m = QCELL(lev, i, j, k, workSet, l); 
1137
1125
                                                rho += m;
1138
 
                                                ux  += (this->dfDvecX[l]*m);
1139
 
                                                uy  += (this->dfDvecY[l]*m); 
1140
 
                                                uz  += (this->dfDvecZ[l]*m); 
 
1126
                                                ux  += (dfDvecX[l]*m);
 
1127
                                                uy  += (dfDvecY[l]*m); 
 
1128
                                                uz  += (dfDvecZ[l]*m); 
1141
1129
                                        } 
1142
1130
                                        rhoOld = rho;
1143
1131
                                        velOld = LbmVec(ux,uy,uz);
1148
1136
                                        LbmFloat df[LBM_DFNUM];
1149
1137
                                        LbmFloat feqOld[LBM_DFNUM];
1150
1138
                                        LbmFloat feqNew[LBM_DFNUM];
1151
 
                                        for(int l=0; l<this->cDfNum; l++) {
1152
 
                                                feqOld[l] = this->getCollideEq(l,rhoOld, velOld[0],velOld[1],velOld[2] );
1153
 
                                                feqNew[l] = this->getCollideEq(l,rhoNew, velNew[0],velNew[1],velNew[2] );
 
1139
                                        for(int l=0; l<cDfNum; l++) {
 
1140
                                                feqOld[l] = getCollideEq(l,rhoOld, velOld[0],velOld[1],velOld[2] );
 
1141
                                                feqNew[l] = getCollideEq(l,rhoNew, velNew[0],velNew[1],velNew[2] );
1154
1142
                                                df[l] = QCELL(lev, i,j,k,workSet, l);
1155
1143
                                        }
1156
 
                                        const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feqOld);
1157
 
                                        const LbmFloat oldOmega = this->getLesOmega(levOldOmega[lev], mLevel[lev].lcsmago,Qo);
1158
 
                                        const LbmFloat newOmega = this->getLesOmega(mLevel[lev].omega,mLevel[lev].lcsmago,Qo);
 
1144
                                        const LbmFloat Qo = getLesNoneqTensorCoeff(df,feqOld);
 
1145
                                        const LbmFloat oldOmega = getLesOmega(levOldOmega[lev], mLevel[lev].lcsmago,Qo);
 
1146
                                        const LbmFloat newOmega = getLesOmega(mLevel[lev].omega,mLevel[lev].lcsmago,Qo);
1159
1147
                                        //newOmega = mLevel[lev].omega; // FIXME debug test
1160
1148
 
1161
1149
                                        //LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]);
1162
1150
                                        const LbmFloat dfScale = (newSteptime/newOmega)/(levOldStepsize[lev]/oldOmega);
1163
1151
                                        //dfScale = dfScaleFac/newOmega;
1164
1152
                                        
1165
 
                                        for(int l=0; l<this->cDfNum; l++) {
 
1153
                                        for(int l=0; l<cDfNum; l++) {
1166
1154
                                                // org scaling
1167
1155
                                                //df = eqOld + (df-eqOld)*dfScale; df *= (eqNew/eqOld); // non-eq. scaling, important
1168
1156
                                                // new scaling
1206
1194
 
1207
1195
                } // lev
1208
1196
 
1209
 
                if(!this->mSilent) {
1210
 
                        debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<this->mStepCnt<<
 
1197
                if(!mSilent) {
 
1198
                        debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<mStepCnt<<
1211
1199
                                        " no"<<mTimeSwitchCounts<<" maxdt"<<mMaxTimestep<<
1212
1200
                                        " mindt"<<mMinTimestep<<" currdt"<<mLevel[mMaxRefine].timestep, 10);
1213
1201
                        debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE  masst:"<<massTNew<<","<<massTOld<<" org:"<<mCurrentMass<<"; "<<
1214
1202
                                        " volt:"<<volTNew<<","<<volTOld<<" org:"<<mCurrentVolume, 10);
1215
1203
                } else {
1216
 
                        debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep change by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<<newdt
1217
 
                                        <<", oldDt:"<<levOldStepsize[mMaxRefine]<<" newOmega:"<<this->mOmega<<" gStar:"<<this->mpParam->getCurrentGStar() , 10);
 
1204
                        debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep changed by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<<newdt
 
1205
                                        <<", oldDt:"<<levOldStepsize[mMaxRefine]<<" newOmega:"<<mOmega<<" gStar:"<<mpParam->getCurrentGStar()<<", step:"<<mStepCnt , 10);
1218
1206
                }
1219
1207
        } // rescale?
1220
1208
        //NEWDEBCHECK("tt2");
1221
1209
        
1222
1210
        //errMsg("adaptTimestep","Warning - brute force rescale off!"); minCutoff = false; // DEBUG
1223
1211
        if(minCutoff) {
1224
 
                errMsg("adaptTimestep","Warning - performing Brute-Force rescale... (sim:"<<this->mName<<" step:"<<this->mStepCnt<<" newdt="<<desireddt<<" mindt="<<this->mpParam->getMinTimestep()<<") " );
 
1212
                errMsg("adaptTimestep","Warning - performing Brute-Force rescale... (sim:"<<mName<<" step:"<<mStepCnt<<" newdt="<<desireddt<<" mindt="<<mpParam->getMinTimestep()<<") " );
1225
1213
                //brute force resacle all the time?
1226
1214
 
1227
1215
                for(int lev=mMaxRefine; lev>=0 ; lev--) {
1250
1238
                        // collide on current set
1251
1239
                        LbmFloat rho, ux,uy,uz;
1252
1240
                        rho=0.0; ux =  uy = uz = 0.0;
1253
 
                        for(int l=0; l<this->cDfNum; l++) {
 
1241
                        for(int l=0; l<cDfNum; l++) {
1254
1242
                                LbmFloat m = QCELL(lev, i, j, k, workSet, l); 
1255
1243
                                rho += m;
1256
 
                                ux  += (this->dfDvecX[l]*m);
1257
 
                                uy  += (this->dfDvecY[l]*m); 
1258
 
                                uz  += (this->dfDvecZ[l]*m); 
 
1244
                                ux  += (dfDvecX[l]*m);
 
1245
                                uy  += (dfDvecY[l]*m); 
 
1246
                                uz  += (dfDvecZ[l]*m); 
1259
1247
                        } 
1260
1248
#ifndef WIN32
1261
1249
                        if (!finite(rho)) {
1272
1260
                                ux *= cfac;
1273
1261
                                uy *= cfac;
1274
1262
                                uz *= cfac;
1275
 
                                for(int l=0; l<this->cDfNum; l++) {
1276
 
                                        QCELL(lev, i, j, k, workSet, l) = this->getCollideEq(l, rho, ux,uy,uz); }
 
1263
                                for(int l=0; l<cDfNum; l++) {
 
1264
                                        QCELL(lev, i, j, k, workSet, l) = getCollideEq(l, rho, ux,uy,uz); }
1277
1265
                                rescs++;
1278
1266
                                debMsgDirect("B");
1279
1267
                        }