apply first part of patch #6994 - elbeem_warning_patch.diff
[blender.git] / intern / elbeem / intern / parametrizer.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2003-2006 Nils Thuerey
5  *
6  * Parameter calculator for the LBM Solver class
7  *
8  *****************************************************************************/
9
10 #include <sstream>
11 #include "parametrizer.h"
12
13 // debug output flag, has to be off for win32 for some reason...
14 #define DEBUG_PARAMCHANNELS 0
15
16 /*! param seen debug string array */
17 const char *ParamStrings[] = {
18         "RelaxTime",
19         "Reynolds",
20         "Viscosity",
21         "SoundSpeed",
22         "DomainSize",
23         "GravityForce",
24         "TimeLength",
25         "Timestep",
26         "Size",
27         "TimeFactor",
28         "AniFrames",
29         "AniFrameTime",
30         "AniStart",
31         "SurfaceTension",
32         "Density",
33         "CellSize",
34         "GStar",
35         "MaxSpeed",
36         "SimMaxSpeed",
37         "FluidVolHeight",
38         "NormalizedGStar",
39         "PSERR", "PSERR", "PSERR", "PSERR"
40 };
41
42
43
44 /******************************************************************************
45  * Default constructor
46  *****************************************************************************/
47 Parametrizer::Parametrizer( void ) :
48         mcViscosity( 8.94e-7 ), 
49         mSoundSpeed( 1500 ),
50         mDomainSize( 0.1 ), mCellSize( 0.01 ),
51         mcGravity( ParamVec(0.0) ),
52         mTimestep(0.0001), mDesiredTimestep(-1.0),
53         mMaxTimestep(-1.0),
54         mMinTimestep(-1.0), 
55         mSizex(50), mSizey(50), mSizez(50),
56         mTimeFactor( 1.0 ),
57         mcAniFrameTime(0.0001),
58         mTimeStepScale(1.0),
59         mAniStart(0.0),
60         //mExtent(1.0, 1.0, 1.0), //mSurfaceTension( 0.0 ),
61         mDensity(1000.0), mGStar(0.0001), mFluidVolumeHeight(0.0),
62         mSimulationMaxSpeed(0.0),
63         mTadapMaxOmega(2.0), mTadapMaxSpeed(0.1), mTadapLevels(1),
64         mFrameNum(0),
65         mSeenValues( 0 ), mCalculatedValues( 0 )
66 {
67 }
68
69
70 /******************************************************************************
71  * Destructor
72  *****************************************************************************/
73 Parametrizer::~Parametrizer() 
74 {
75         /* not much to do... */
76 }
77
78 /******************************************************************************
79  * Init from attr list
80  *****************************************************************************/
81 void Parametrizer::parseAttrList() 
82 {
83         if(!mpAttrs) {
84                 errFatal("Parametrizer::parseAttrList", "mpAttrs pointer not initialized!", SIMWORLD_INITERROR);
85                 return;
86         }
87
88         // unused
89         string  mSetupType = "";
90         mSetupType = mpAttrs->readString("p_setup",mSetupType, "Parametrizer","mSetupType", false); 
91
92         // real params
93         if(getAttributeList()->exists("p_viscosity")) {
94                         mcViscosity = mpAttrs->readChannelFloat("p_viscosity"); seenThis( PARAM_VISCOSITY ); }
95
96         mSoundSpeed = mpAttrs->readFloat("p_soundspeed",mSoundSpeed, "Parametrizer","mSoundSpeed", false); 
97         if(getAttributeList()->exists("p_soundspeed")) seenThis( PARAM_SOUNDSPEED );
98
99         mDomainSize = mpAttrs->readFloat("p_domainsize",mDomainSize, "Parametrizer","mDomainSize", false); 
100         if(getAttributeList()->exists("p_domainsize")) seenThis( PARAM_DOMAINSIZE );
101         if(mDomainSize<=0.0) {
102                 errMsg("Parametrizer::parseAttrList","Invalid real world domain size:"<<mDomainSize<<", resetting to 0.1");
103                 mDomainSize = 0.1;
104         }
105
106         if(getAttributeList()->exists("p_gravity")) { // || (!mcGravity.isInited()) ) {
107                 mcGravity = mpAttrs->readChannelVec3d("p_gravity"); seenThis( PARAM_GRAVITY );
108         }
109
110         mTimestep = mpAttrs->readFloat("p_steptime",mTimestep, "Parametrizer","mTimestep", false); 
111         if(getAttributeList()->exists("p_steptime")) seenThis( PARAM_STEPTIME );
112
113         mTimeFactor = mpAttrs->readFloat("p_timefactor",mTimeFactor, "Parametrizer","mTimeFactor", false); 
114         if(getAttributeList()->exists("p_timefactor")) seenThis( PARAM_TIMEFACTOR );
115
116         if(getAttributeList()->exists("p_aniframetime")) { //|| (!mcAniFrameTime.isInited()) ) {
117                 mcAniFrameTime = mpAttrs->readChannelFloat("p_aniframetime");seenThis( PARAM_ANIFRAMETIME ); 
118         }
119         mTimeStepScale = mpAttrs->readFloat("p_timestepscale",mTimeStepScale, "Parametrizer","mTimeStepScale", false); 
120
121         mAniStart = mpAttrs->readFloat("p_anistart",mAniStart, "Parametrizer","mAniStart", false); 
122         if(getAttributeList()->exists("p_anistart")) seenThis( PARAM_ANISTART );
123         if(mAniStart<0.0) {
124                 errMsg("Parametrizer::parseAttrList","Invalid start time:"<<mAniStart<<", resetting to 0.0");
125                 mAniStart = 0.0;
126         }
127
128         //mSurfaceTension = mpAttrs->readFloat("p_surfacetension",mSurfaceTension, "Parametrizer","mSurfaceTension", false); 
129         //if(getAttributeList()->exists("p_surfacetension")) seenThis( PARAM_SURFACETENSION );
130
131         mDensity = mpAttrs->readFloat("p_density",mDensity, "Parametrizer","mDensity", false); 
132         if(getAttributeList()->exists("p_density")) seenThis( PARAM_DENSITY );
133
134         ParamFloat cellSize = 0.0; // unused, deprecated
135         cellSize = mpAttrs->readFloat("p_cellsize",cellSize, "Parametrizer","cellSize", false); 
136
137         mGStar = mpAttrs->readFloat("p_gstar",mGStar, "Parametrizer","mGStar", false); 
138         if(getAttributeList()->exists("p_gstar")) seenThis( PARAM_GSTAR );
139
140         mNormalizedGStar = mpAttrs->readFloat("p_normgstar",mNormalizedGStar, "Parametrizer","mNormalizedGStar", false); 
141         if(getAttributeList()->exists("p_normgstar")) seenThis( PARAM_NORMALIZEDGSTAR );
142
143         mTadapMaxOmega = mpAttrs->readFloat("p_tadapmaxomega",mTadapMaxOmega, "Parametrizer","mTadapMaxOmega", false); 
144         mTadapMaxSpeed = mpAttrs->readFloat("p_tadapmaxspeed",mTadapMaxSpeed, "Parametrizer","mTadapMaxSpeed", false); 
145 }
146
147 /******************************************************************************
148  *! advance to next render/output frame 
149  *****************************************************************************/
150 void Parametrizer::setFrameNum(int frame) {
151         mFrameNum = frame;
152 #if DEBUG_PARAMCHANNELS>0
153         errMsg("DEBUG_PARAMCHANNELS","setFrameNum frame-num="<<mFrameNum);
154 #endif // DEBUG_PARAMCHANNELS>0
155 }
156 /*! get time of an animation frame (renderer)  */
157 // also used by: mpParam->getCurrentAniFrameTime() , e.g. for velocity dump
158 ParamFloat Parametrizer::getAniFrameTime( int frame )   { 
159         double frametime = (double)frame;
160         ParamFloat anift = mcAniFrameTime.get(frametime);
161         if(anift<0.0) {
162                 ParamFloat resetv = 0.;
163                 errMsg("Parametrizer::setFrameNum","Invalid frame time:"<<anift<<" at frame "<<frame<<", resetting to "<<resetv);
164                 anift = resetv; 
165         }
166 #if DEBUG_PARAMCHANNELS>0
167         if((0)|| (DEBUG_PARAMCHANNELS)) errMsg("DEBUG_PARAMCHANNELS","getAniFrameTime frame="<<frame<<", frametime="<<anift<<" ");
168 #endif // DEBUG_PARAMCHANNELS>0
169         return anift; 
170 }
171
172 /******************************************************************************
173  * scale a given speed vector in m/s to lattice values 
174  *****************************************************************************/
175 ParamVec Parametrizer::calculateAddForce(ParamVec vec, string usage)
176 {
177         ParamVec ret = vec * (mTimestep*mTimestep) /mCellSize;
178         debMsgStd("Parametrizer::calculateVector", DM_MSG, "scaled vector = "<<ret<<" for '"<<usage<<"', org = "<<vec<<" dt="<<mTimestep ,10);
179         return ret;
180 }
181
182
183 /******************************************************************************
184  * calculate size of a single cell
185  *****************************************************************************/
186 ParamFloat Parametrizer::calculateCellSize(void)
187 {
188         int maxsize = mSizex; // get max size
189         if(mSizey>maxsize) maxsize = mSizey;
190         if(mSizez>maxsize) maxsize = mSizez;
191         maxsize = mSizez; // take along gravity dir for now!
192         ParamFloat cellSize = 1.0 / (ParamFloat)maxsize;
193         return cellSize;
194 }
195
196
197 /*****************************************************************************/
198 /* simple calulation functions */
199 /*****************************************************************************/
200
201 /*! get omega for LBM from channel */
202 ParamFloat Parametrizer::calculateOmega( double time ) { 
203         ParamFloat viscStar = calculateLatticeViscosity(time);
204         ParamFloat relaxTime = (6.0 * viscStar + 1) * 0.5;
205 #if DEBUG_PARAMCHANNELS>0
206         errMsg("DEBUG_PARAMCHANNELS","calculateOmega viscStar="<<viscStar<<" relaxtime="<<relaxTime);
207 #endif // DEBUG_PARAMCHANNELS>0
208         return (1.0/relaxTime); 
209 }
210
211 /*! get external force x component */
212 ParamVec Parametrizer::calculateGravity( double time ) { 
213         ParamVec grav = mcGravity.get(time);
214         ParamFloat forceFactor = (mTimestep *mTimestep)/mCellSize;
215         ParamVec latticeGravity = grav * forceFactor;
216 #if DEBUG_PARAMCHANNELS>0
217         errMsg("DEBUG_PARAMCHANNELS","calculateGravity grav="<<grav<<" ff"<<forceFactor<<" lattGrav="<<latticeGravity);
218 #endif // DEBUG_PARAMCHANNELS>0
219         return latticeGravity; 
220 }
221
222 /*! calculate the lattice viscosity */
223 ParamFloat Parametrizer::calculateLatticeViscosity( double time ) { 
224         // check seen values
225         int reqValues = PARAM_VISCOSITY | PARAM_STEPTIME;
226         if(!checkSeenValues( reqValues ) ){
227                 errMsg("Parametrizer::calculateLatticeViscosity"," Missing arguments!");
228         }
229         ParamFloat viscStar = mcViscosity.get(time) * mTimestep / (mCellSize*mCellSize);
230 #if DEBUG_PARAMCHANNELS>0
231         errMsg("DEBUG_PARAMCHANNELS","calculateLatticeViscosity viscStar="<<viscStar);
232 #endif // DEBUG_PARAMCHANNELS>0
233         return viscStar; 
234 }
235
236 /*! get no of steps for the given length in seconds */
237 int Parametrizer::calculateStepsForSecs( ParamFloat s ) { 
238         return (int)(s/mTimestep); 
239 }
240
241 /*! get start time of animation */
242 int Parametrizer::calculateAniStart( void )   { 
243         return (int)(mAniStart/mTimestep); 
244 }
245
246 /*! get no of steps for a single animation frame */
247 int Parametrizer::calculateAniStepsPerFrame(int frame)   { 
248         if(!checkSeenValues(PARAM_ANIFRAMETIME)) {
249                 errFatal("Parametrizer::calculateAniStepsPerFrame", "Missing ani frame time argument!", SIMWORLD_INITERROR);
250                 return 1;
251         }
252         int value = (int)(getAniFrameTime(frame)/mTimestep); 
253         if((value<0) || (value>1000000)) {
254                 errFatal("Parametrizer::calculateAniStepsPerFrame", "Invalid step-time (="<<value<<") <> ani-frame-time ("<<mTimestep<<") settings, aborting...", SIMWORLD_INITERROR);
255                 return 1;
256         }
257         return value;
258 }
259
260 /*! get no. of timesteps for LBM */
261 int Parametrizer::calculateNoOfSteps( ParamFloat timelen ) { 
262         return (int)(timelen/mTimestep); 
263 }
264
265 /*! get extent of the domain = (1,1,1) if parametrizer not used, (x,y,z) [m] otherwise */
266 //ParamVec Parametrizer::calculateExtent( void ) { 
267         //return mExtent; 
268 //}
269
270 /*! get (scaled) surface tension */
271 //ParamFloat Parametrizer::calculateSurfaceTension( void ) { return mSurfaceTension; }
272
273 /*! get the length of a single time step */
274 // explicity scaled by time factor for refinement 
275 ParamFloat Parametrizer::getTimestep( void ) { 
276         return mTimestep; 
277 }
278
279 /*! calculate lattice velocity from real world value [m/s] */
280 ParamVec Parametrizer::calculateLattVelocityFromRw( ParamVec ivel ) { 
281         ParamVec velvec = ivel;
282         velvec /= mCellSize;
283         velvec *= mTimestep;
284         return velvec; 
285 }
286 /*! calculate real world [m/s] velocity from lattice value */
287 ParamVec Parametrizer::calculateRwVelocityFromLatt( ParamVec ivel ) { 
288         ParamVec velvec = ivel;
289         velvec *= mCellSize;
290         velvec /= mTimestep;
291         return velvec; 
292 }
293
294
295 /*! get g star value with fhvol calculations */
296 ParamFloat Parametrizer::getCurrentGStar( void ) {
297         ParamFloat gStar = mGStar; // check? TODO get from mNormalizedGStar?
298         if(mFluidVolumeHeight>0.0) { gStar = mGStar/mFluidVolumeHeight; }
299         return gStar;
300 }
301
302 /******************************************************************************
303  * function that tries to calculate all the missing values from the given ones
304  * prints errors and returns false if thats not possible 
305  *****************************************************************************/
306 bool Parametrizer::calculateAllMissingValues( double time, bool silent )
307 {
308         bool init = false;  // did we init correctly?
309         int  valuesChecked = 0;
310         int reqValues;
311
312         // we always need the sizes
313         reqValues = PARAM_SIZE;
314         valuesChecked |= reqValues;
315         if(!checkSeenValues(reqValues)) {
316                 errMsg("Parametrizer::calculateAllMissingValues"," Missing size argument!");
317                 return false;
318         }
319
320         if(!checkSeenValues(PARAM_DOMAINSIZE)) {
321                 errMsg("Parametrizer::calculateAllMissingValues"," Missing domain size argument!");
322                 return false;
323         }
324         int maxsize = mSizex; // get max size
325         if(mSizey>maxsize) maxsize = mSizey;
326         if(mSizez>maxsize) maxsize = mSizez;
327         maxsize = mSizez; // take along gravity dir for now!
328         mCellSize = ( mDomainSize * calculateCellSize() ); // sets mCellSize
329         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," max domain resolution="<<(maxsize)<<" cells , cellsize="<<mCellSize ,10);
330
331                         
332         /* Carolin init , see DA for details */
333         ParamFloat maxDeltaT = 0.0;
334
335         /* normalized gstar init */
336         reqValues = PARAM_NORMALIZEDGSTAR;
337         valuesChecked |= reqValues;
338         if(checkSeenValues( reqValues ) ){
339                 //if(checkSeenValues( PARAM_GSTAR ) ){ if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_WARNING," g star value override by normalizedGStar!",1); }
340                 const ParamFloat normgstarReset = 0.005;
341                 if(mNormalizedGStar<=1e-6) {
342                         errMsg("Parametrizer::calculateAllMissingValues","Invalid NormGstar: "<<mNormalizedGStar<<"... resetting to "<<normgstarReset);
343                         mNormalizedGStar = normgstarReset;
344                 }
345
346                 mGStar = mNormalizedGStar/maxsize;
347
348 // TODO FIXME add use testdata check!
349 mGStar = mNormalizedGStar/mSizez;
350 errMsg("Warning","Used z-dir for gstar!");
351
352                 if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," g star set to "<<mGStar<<" from normalizedGStar="<<mNormalizedGStar ,1);
353                 seenThis(PARAM_GSTAR);
354         }
355
356         reqValues = PARAM_GSTAR | PARAM_VISCOSITY;
357         if((checkSeenValues(PARAM_SURFACETENSION))) reqValues |= PARAM_DENSITY; // surface tension optional now...
358         valuesChecked |= reqValues;
359         if(checkSeenValues( reqValues ) ){
360                 const ParamFloat gstarReset = 0.0005;
361                 if(getCurrentGStar()<=1e-6) {
362                         errMsg("Parametrizer::calculateAllMissingValues","Invalid Gstar: "<<getCurrentGStar()<<" (set to "<<mGStar<<") ... resetting to "<<gstarReset);
363                         mGStar = gstarReset;
364                 }
365
366                 ParamFloat gStar = getCurrentGStar(); // mGStar
367                 if(mFluidVolumeHeight>0.0) {
368                         debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," height"<<mFluidVolumeHeight<<" resGStar = "<<gStar, 10);
369                 }
370                 if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," g star = "<<gStar, 10);
371
372                 ParamFloat forceStrength = 0.0;
373                 //if(checkSeenValues(PARAM_GRAVITY)) { forceStrength = norm( calculateGravity(time) ); }
374                 if(checkSeenValues(PARAM_GRAVITY)) { forceStrength = norm( mcGravity.get(time) ); }
375
376                 // determine max. delta density per timestep trough gravity force
377                 if(forceStrength>0.0) {
378                         maxDeltaT = sqrt( gStar*mCellSize *mTimeStepScale /forceStrength );
379                 } else {
380                         // use 1 lbm setp = 1 anim step as max
381                         maxDeltaT = getAniFrameTime(0);
382                 }
383
384                 if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," targeted step time = "<<maxDeltaT, 10);
385
386                 //ParamFloat viscStarFac = mViscosity/(mCellSize*mCellSize);
387                 //if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," viscStarFac = "<<viscStarFac<<" viscosity:"<<mViscosity, 10);
388
389                 // time step adaptivty, only for caro with max sim speed
390                 ParamFloat setDeltaT = maxDeltaT;
391                 if(mDesiredTimestep>0.0) {
392                         // explicitly set step time according to max velocity in sim
393                         setDeltaT = mDesiredTimestep;
394                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," desired step time = "<<setDeltaT, 10);
395                         mDesiredTimestep = -1.0;
396                 } else {
397                         // just use max delta t as current
398                 }
399                 
400                 // and once for init determine minimal delta t by omega max. 
401                 if((mMinTimestep<0.0) || (mMaxTimestep<0.0)) {
402                         ParamFloat minDeltaT; 
403                         ParamFloat maxOmega = mTadapMaxOmega;
404                         ParamFloat minRelaxTime = 1.0/maxOmega;
405                         for(int lev=1; lev<mTadapLevels; lev++) {
406                                 // make minRelaxTime larger for each level that exists...
407                                 minRelaxTime = 2.0 * (minRelaxTime-0.5) + 0.5;
408                         }
409                         maxOmega = 1.0/minRelaxTime;
410                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," maxOmega="<<maxOmega<<" minRelaxTime="<<minRelaxTime<<" levels="<<mTadapLevels, 1);
411                         // visc-star for min relax time to calculate min delta ta
412                         if(mcViscosity.get(time)>0.0) {
413                                 minDeltaT = ((2.0*minRelaxTime-1.0)/6.0) * mCellSize * mCellSize / mcViscosity.get(time);
414                         } else {
415                                 // visc=0, this is not physical, but might happen
416                                 minDeltaT = 0.0;
417                         }
418                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," min delta t = "<<minDeltaT<<" , range = " << (maxDeltaT/minDeltaT) ,1);
419
420                         // sim speed + accel shouldnt exceed 0.1?
421                         mMaxTimestep = maxDeltaT;
422                         mMinTimestep = minDeltaT;
423                         // only use once...  
424                 } 
425
426                 setTimestep( setDeltaT ); // set mTimestep to new value
427                 init = true;    
428         }
429
430         // finish init
431         if(init) {
432                 if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," omega = "<<calculateOmega(0.0)<<", delt="<<mTimestep,1);
433
434                 if(checkSeenValues(PARAM_GRAVITY)) {
435                         ParamFloat forceFactor = (mTimestep *mTimestep)/mCellSize; // only used for printing...
436                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," gravity force = "<<PRINT_NTLVEC(mcGravity.get(time))<<", scaled with "<<forceFactor<<" to "<<calculateGravity(time),1);
437                 }
438                 
439                 //mExtent = ParamVec( mCellSize*mSizex, mCellSize*mSizey, mCellSize*mSizez );
440                 //if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," domain extent = "<<PRINT_NTLVEC(mExtent)<<"m , gs:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<" cs:"<<mCellSize,1);
441                 
442                 if(!checkSeenValues(PARAM_ANIFRAMETIME)) {
443                         errFatal("Parametrizer::calculateAllMissingValues"," Warning no ani frame time given!", SIMWORLD_INITERROR);
444                         setAniFrameTimeChannel( mTimestep );
445                 } 
446
447                 if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," ani frame steps = "<<calculateAniStepsPerFrame(mFrameNum)<<" for frame "<<mFrameNum, 1);
448
449                 if((checkSeenValues(PARAM_ANISTART))&&(calculateAniStart()>0)) {
450                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," ani start steps = "<<calculateAniStart()<<" ",1); 
451                 }
452                         
453                 if(! isSimworldOk() ) return false;
454                 // everything ok
455                 return true;
456         }
457
458         // init failed ... failure:
459         errMsg("Parametrizer::calculateAllMissingValues "," invalid configuration!");
460         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues ",DM_WARNING, " values seen:", 1);
461         for(int i=0;i<PARAM_NUMIDS;i++) {
462                 if(checkSeenValues( 1<<i )) {
463                         if(!silent) debMsgStd("  ",DM_NOTIFY, ParamStrings[i], 1);
464                 }
465         }
466         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues ",DM_WARNING, "values checked but missing:", 1);
467         for(int i=0;i<PARAM_NUMIDS;i++) {
468                 if((!checkSeenValues( 1<<i ))&&
469                          ( (valuesChecked&(1<<i))==(1<<i)) ) {
470                         debMsgStd("  ",DM_IMPORTANT, ParamStrings[i], 1);
471                 }
472         }
473
474         // print values?
475         return false;
476 }
477
478
479 /******************************************************************************
480  * init debug functions
481  *****************************************************************************/
482
483
484 /*! set kinematic viscosity */
485 void Parametrizer::setViscosity(ParamFloat set) { 
486         mcViscosity = AnimChannel<ParamFloat>(set); 
487         seenThis( PARAM_VISCOSITY ); 
488 #if DEBUG_PARAMCHANNELS>0
489         { errMsg("DebugChannels","Parametrizer::mcViscosity set = "<< mcViscosity.printChannel() ); }
490 #endif // DEBUG_PARAMCHANNELS>0
491 }
492 void Parametrizer::initViscosityChannel(vector<ParamFloat> val, vector<double> time) { 
493         mcViscosity = AnimChannel<ParamFloat>(val,time); 
494         seenThis( PARAM_VISCOSITY ); 
495 #if DEBUG_PARAMCHANNELS>0
496         { errMsg("DebugChannels","Parametrizer::mcViscosity initc = "<< mcViscosity.printChannel() ); }
497 #endif // DEBUG_PARAMCHANNELS>0
498 }
499
500 /*! set the external force */
501 void Parametrizer::setGravity(ParamFloat setx, ParamFloat sety, ParamFloat setz) { 
502         mcGravity = AnimChannel<ParamVec>(ParamVec(setx,sety,setz)); 
503         seenThis( PARAM_GRAVITY ); 
504 #if DEBUG_PARAMCHANNELS>0
505         { errMsg("DebugChannels","Parametrizer::mcGravity set = "<< mcGravity.printChannel() ); }
506 #endif // DEBUG_PARAMCHANNELS>0
507 }
508 void Parametrizer::setGravity(ParamVec set) { 
509         mcGravity = AnimChannel<ParamVec>(set); 
510         seenThis( PARAM_GRAVITY ); 
511 #if DEBUG_PARAMCHANNELS>0
512         { errMsg("DebugChannels","Parametrizer::mcGravity set = "<< mcGravity.printChannel() ); }
513 #endif // DEBUG_PARAMCHANNELS>0
514 }
515 void Parametrizer::initGravityChannel(vector<ParamVec> val, vector<double> time) { 
516         mcGravity = AnimChannel<ParamVec>(val,time); 
517         seenThis( PARAM_GRAVITY ); 
518 #if DEBUG_PARAMCHANNELS>0
519         { errMsg("DebugChannels","Parametrizer::mcGravity initc = "<< mcGravity.printChannel() ); }
520 #endif // DEBUG_PARAMCHANNELS>0
521 }
522
523 /*! set time of an animation frame (renderer)  */
524 void Parametrizer::setAniFrameTimeChannel(ParamFloat set) { 
525         mcAniFrameTime = AnimChannel<ParamFloat>(set); 
526         seenThis( PARAM_ANIFRAMETIME ); 
527 #if DEBUG_PARAMCHANNELS>0
528         { errMsg("DebugChannels","Parametrizer::mcAniFrameTime set = "<< mcAniFrameTime.printChannel() ); }
529 #endif // DEBUG_PARAMCHANNELS>0
530 }
531 void Parametrizer::initAniFrameTimeChannel(vector<ParamFloat> val, vector<double> time) { 
532         mcAniFrameTime = AnimChannel<ParamFloat>(val,time); 
533         seenThis( PARAM_ANIFRAMETIME ); 
534 #if DEBUG_PARAMCHANNELS>0
535         { errMsg("DebugChannels","Parametrizer::mcAniFrameTime initc = "<< mcAniFrameTime.printChannel() ); }
536 #endif // DEBUG_PARAMCHANNELS>0
537 }
538
539 // OLD interface stuff
540 // reactivate at some point?
541
542                 /*! surface tension, [kg/s^2] */
543                 //ParamFloat mSurfaceTension;
544                 /*! set starting time of the animation (renderer) */
545                 //void setSurfaceTension(ParamFloat set) { mSurfaceTension = set; seenThis( PARAM_SURFACETENSION ); }
546                 /*! get starting time of the animation (renderer) */
547                 //ParamFloat getSurfaceTension( void )   { return mSurfaceTension; }
548                 /*if((checkSeenValues(PARAM_SURFACETENSION))&&(mSurfaceTension>0.0)) {
549                         ParamFloat massDelta = 1.0;
550                         ParamFloat densityStar = 1.0;
551                         massDelta = mDensity / densityStar *mCellSize*mCellSize*mCellSize;
552                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," massDelta = "<<massDelta, 10);
553
554                         mSurfaceTension = mSurfaceTension*mTimestep*mTimestep/massDelta;
555                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," surface tension = "<<mSurfaceTension<<" ",1);
556                 } // */
557
558 // probably just delete:
559
560                 /*! reynolds number (calculated from domain length and max. speed [dimensionless] */
561                 //ParamFloat mReynolds;
562
563                 /*! set relaxation time */
564                 //void setRelaxTime(ParamFloat set) { mRelaxTime = set; seenThis( PARAM_RELAXTIME ); }
565                 /*! get relaxation time */
566                 //ParamFloat getRelaxTime( void )   { return mRelaxTime; }
567                 /*! set reynolds number */
568                 //void setReynolds(ParamFloat set) { mReynolds = set; seenThis( PARAM_REYNOLDS ); }
569                 /*! get reynolds number */
570                 //ParamFloat getReynolds( void )   { return mReynolds; }
571
572                 // calculate reynolds number
573                 /*if(mViscosity>0.0) {
574                         ParamFloat maxSpeed  = 1.0/6.0; // for rough reynolds approx
575                         ParamFloat reynoldsApprox = -1.0;
576                         ParamFloat gridSpeed = (maxSpeed*mCellSize/mTimestep);
577                         reynoldsApprox = (mDomainSize*gridSpeed) / mViscosity;
578                         if(!silent) debMsgStd("Parametrizer::calculateAllMissingValues",DM_MSG," reynolds number (D="<<mDomainSize<<", assuming V="<<gridSpeed<<")= "<<reynoldsApprox<<" ", 1);
579                 } // */
580
581         //? mRelaxTime = mpAttrs->readFloat("p_relaxtime",mRelaxTime, "Parametrizer","mRelaxTime", false); 
582         //if(getAttributeList()->exists("p_relaxtime")) seenThis( PARAM_RELAXTIME );
583         //? mReynolds = mpAttrs->readFloat("p_reynolds",mReynolds, "Parametrizer","mReynolds", false); 
584         //if(getAttributeList()->exists("p_reynolds")) seenThis( PARAM_REYNOLDS );
585
586         //mViscosity = mpAttrs->readFloat("p_viscosity",mViscosity, "Parametrizer","mViscosity", false); 
587         //if(getAttributeList()->exists("p_viscosity") || (!mcViscosity.isInited()) ) { }
588         //if(getAttributeList()->exists("p_viscosity")) 
589
590
591                 //ParamFloat viscStar = calculateLatticeViscosity(time);
592                 //RelaxTime = (6.0 * viscStar + 1) * 0.5;
593
594