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