[Vapor-cvs] vapor/lib/flow VTFieldLine.cpp, 1.39, 1.40 VTFieldLine.h, 1.28, 1.29 VTStreamLine.cpp, 1.29, 1.30 VTTimeVaryingFieldLine.cpp, 1.19, 1.20 VectorMatrix.h, 1.3, 1.4

[Alan N. <ala...@us...>]

Update of /cvsroot/vapor/vapor/lib/flow
In directory sc8-pr-cvs9.sourceforge.net:/tmp/cvs-serv24678/lib/flow

Modified Files:
	VTFieldLine.cpp VTFieldLine.h VTStreamLine.cpp 
	VTTimeVaryingFieldLine.cpp VectorMatrix.h 
Log Message:
Fix for 1776639.
Also fixed a problem with integrating flow with very asymmetric domains,
by allowing very large changes in step size.


Index: VTTimeVaryingFieldLine.cpp
===================================================================
RCS file: /cvsroot/vapor/vapor/lib/flow/VTTimeVaryingFieldLine.cpp,v
retrieving revision 1.19
retrieving revision 1.20
diff -C2 -d -r1.19 -r1.20
*** VTTimeVaryingFieldLine.cpp	14 Aug 2007 17:44:26 -0000	1.19
--- VTTimeVaryingFieldLine.cpp	28 Aug 2007 21:05:43 -0000	1.20
***************
*** 116,120 ****
  	PointInfo thisParticle;
  	VECTOR3 thisInterpolant, prevInterpolant, second_prevInterpolant;
! 	float dt, cell_volume, mag; 
  	double curTime;
  	VECTOR3 vel;
--- 116,120 ----
  	PointInfo thisParticle;
  	VECTOR3 thisInterpolant, prevInterpolant, second_prevInterpolant;
! 	double dt, cell_volume, mag; 
  	double curTime;
  	VECTOR3 vel;
***************
*** 133,138 ****
  	// get the initial stepsize
  	cell_volume = m_pField->volume_of_cell(seedInfo.inCell);
! 	mag = vel.GetMag();
! 	dt = pow(cell_volume, (float)0.3333333f) / mag;
  	//Allow dt*mag to be 10 times the initial setting:
  	float maxDtMag = dt*mag*10.f;
--- 133,138 ----
  	// get the initial stepsize
  	cell_volume = m_pField->volume_of_cell(seedInfo.inCell);
! 	mag = vel.GetDMag();
! 	dt = pow(cell_volume, 0.3333333) / mag;
  	//Allow dt*mag to be 10 times the initial setting:
  	float maxDtMag = dt*mag*10.f;
***************
*** 146,154 ****
  	{
  		// how much advection time left
! 		float timeLeft;
  		if (m_timeDir == FORWARD) {
! 			timeLeft = (finalTime - (float)curTime);
  		} else {
! 			timeLeft = (float)curTime - finalTime;
  		}
  		if (dt > timeLeft) dt = timeLeft;
--- 146,154 ----
  	{
  		// how much advection time left
! 		double timeLeft;
  		if (m_timeDir == FORWARD) {
! 			timeLeft = (finalTime - curTime);
  		} else {
! 			timeLeft = curTime - finalTime;
  		}
  		if (dt > timeLeft) dt = timeLeft;
***************
*** 166,170 ****
  			retrace = false;
  			//Loop until dt*mag is small enough, dividing dt by 10 if necessary
! 			for (int tries = 0; tries < 20; tries++){
  				if(int_order == SECOND)
  					istat = runge_kutta2(m_timeDir, UNSTEADY, thisParticle, &curTime, dt,maxDtMag);
--- 166,170 ----
  			retrace = false;
  			//Loop until dt*mag is small enough, dividing dt by 10 if necessary
! 			for (int tries = 0; tries < 40; tries++){
  				if(int_order == SECOND)
  					istat = runge_kutta2(m_timeDir, UNSTEADY, thisParticle, &curTime, dt,maxDtMag);
***************
*** 205,209 ****
  			if(((int)seedTrace.size() > 2) && (bAdaptive))
  			{
! 				float minStepsize, maxStepsize;
  				VECTOR3 thisPhy, prevPhy, second_prevPhy;
  				list<VECTOR3*>::iterator pIter = seedTrace.end();
--- 205,209 ----
  			if(((int)seedTrace.size() > 2) && (bAdaptive))
  			{
! 				double minStepsize, maxStepsize;
  				VECTOR3 thisPhy, prevPhy, second_prevPhy;
  				list<VECTOR3*>::iterator pIter = seedTrace.end();
***************
*** 216,222 ****
  
  				cell_volume = m_pField->volume_of_cell(thisParticle.inCell);
! 				mag = vel.GetMag();
! 				minStepsize = m_fInitStepSize * pow(cell_volume, (float)0.3333333f) / mag;
! 				maxStepsize = m_fMaxStepSize * pow(cell_volume, (float)0.3333333f) / mag;
  				retrace = adapt_step(second_prevPhy, prevPhy, thisPhy, minStepsize, maxStepsize, &dt, bAdaptive);
  				
--- 216,222 ----
  
  				cell_volume = m_pField->volume_of_cell(thisParticle.inCell);
! 				mag = vel.GetDMag();
! 				minStepsize = m_fInitStepSize * pow(cell_volume, 0.3333333) / mag;
! 				maxStepsize = m_fMaxStepSize * pow(cell_volume, 0.3333333) / mag;
  				retrace = adapt_step(second_prevPhy, prevPhy, thisPhy, minStepsize, maxStepsize, &dt, bAdaptive);
  				

Index: VTStreamLine.cpp
===================================================================
RCS file: /cvsroot/vapor/vapor/lib/flow/VTStreamLine.cpp,v
retrieving revision 1.29
retrieving revision 1.30
diff -C2 -d -r1.29 -r1.30
*** VTStreamLine.cpp	16 Aug 2007 16:47:51 -0000	1.29
--- VTStreamLine.cpp	28 Aug 2007 21:05:43 -0000	1.30
***************
*** 188,192 ****
  	PointInfo thisParticle;
  	VECTOR3 thisInterpolant, prevInterpolant, second_prevInterpolant;
! 	float dt, cell_volume, mag; 
  	double curTime;
  	VECTOR3 vel;
--- 188,192 ----
  	PointInfo thisParticle;
  	VECTOR3 thisInterpolant, prevInterpolant, second_prevInterpolant;
! 	double dt, cell_volume, mag; 
  	double curTime;
  	VECTOR3 vel;
***************
*** 209,218 ****
  	// get the initial step size
  	cell_volume = m_pField->volume_of_cell(seedInfo.inCell);
! 	mag = vel.GetMag();
  	
! 	dt = m_fInitStepSize * pow(cell_volume, (float)0.3333333f) / mag;
  	
  	//Determine the value of mag*dt to project 10 times init size.  
! 	float maxMagDt = 10.*dt*mag;
  
  	int rollbackCount = 0;
--- 209,218 ----
  	// get the initial step size
  	cell_volume = m_pField->volume_of_cell(seedInfo.inCell);
! 	mag = vel.GetDMag();
  	
! 	dt = m_fInitStepSize * pow(cell_volume, 0.3333333) / mag;
  	
  	//Determine the value of mag*dt to project 10 times init size.  
! 	double maxMagDt = 10.*dt*mag;
  
  	int rollbackCount = 0;
***************
*** 231,235 ****
  			retrace = false;
  			
! 			for(int magTry = 0; magTry < 20; magTry++) {
  				if(integ_ord == SECOND)
  					istat = runge_kutta2(time_dir, time_dep, thisParticle, &curTime, dt, maxMagDt);
--- 231,235 ----
  			retrace = false;
  			
! 			for(int magTry = 0; magTry < 40; magTry++) {
  				if(integ_ord == SECOND)
  					istat = runge_kutta2(time_dir, time_dep, thisParticle, &curTime, dt, maxMagDt);
***************
*** 264,268 ****
  			if(((int)seedTrace.size() > 2)&&(onAdaptive))
  			{
! 				float minStepsize, maxStepsize;
  				VECTOR3 thisPhy, prevPhy, second_prevPhy;
  				list<VECTOR3*>::iterator pIter = seedTrace.end();
--- 264,268 ----
  			if(((int)seedTrace.size() > 2)&&(onAdaptive))
  			{
! 				double minStepsize, maxStepsize;
  				VECTOR3 thisPhy, prevPhy, second_prevPhy;
  				list<VECTOR3*>::iterator pIter = seedTrace.end();
***************
*** 275,281 ****
  
  				cell_volume = m_pField->volume_of_cell(thisParticle.inCell);
! 				mag = vel.GetMag();
! 				minStepsize = m_fInitStepSize * pow(cell_volume, (float)0.3333333f) / mag;
! 				maxStepsize = m_fMaxStepSize * pow(cell_volume, (float)0.3333333f) / mag;
  				retrace = adapt_step(second_prevPhy, prevPhy, thisPhy, minStepsize, maxStepsize, &dt, onAdaptive);
  				if(onAdaptive == false)
--- 275,281 ----
  
  				cell_volume = m_pField->volume_of_cell(thisParticle.inCell);
! 				mag = vel.GetDMag();
! 				minStepsize = m_fInitStepSize * pow(cell_volume, 0.3333333) / mag;
! 				maxStepsize = m_fMaxStepSize * pow(cell_volume, 0.3333333) / mag;
  				retrace = adapt_step(second_prevPhy, prevPhy, thisPhy, minStepsize, maxStepsize, &dt, onAdaptive);
  				if(onAdaptive == false)

Index: VectorMatrix.h
===================================================================
RCS file: /cvsroot/vapor/vapor/lib/flow/VectorMatrix.h,v
retrieving revision 1.3
retrieving revision 1.4
diff -C2 -d -r1.3 -r1.4
*** VectorMatrix.h	22 Sep 2005 18:21:29 -0000	1.3
--- VectorMatrix.h	28 Aug 2007 21:05:43 -0000	1.4
***************
*** 78,81 ****
--- 78,82 ----
  	void Zero() {vec[0] = vec[1] = vec[2] = 0.0;}					// make zero vector
  	float GetMag() {return (float)(sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]));}	// get magnitude
+ 	double GetDMag() {return (sqrt((double)vec[0]*(double)vec[0] + (double)vec[1]*(double)vec[1] + (double)vec[2]*(double)vec[2]));}	// get magnitude as double
  	const VECTOR3 & operator =(const VECTOR3 & v0)					// copy vector v0
  	{vec[0] = v0(0); vec[1] = v0(1); vec[2] = v0(2); return(*this);}

Index: VTFieldLine.h
===================================================================
RCS file: /cvsroot/vapor/vapor/lib/flow/VTFieldLine.h,v
retrieving revision 1.28
retrieving revision 1.29
diff -C2 -d -r1.28 -r1.29
*** VTFieldLine.h	14 Aug 2007 17:44:25 -0000	1.28
--- VTFieldLine.h	28 Aug 2007 21:05:43 -0000	1.29
***************
*** 135,146 ****
  	void releaseSeedMemory(void);
  	int euler_cauchy(TIME_DIR, TIME_DEP,float*, float);
! 	int runge_kutta4(TIME_DIR, TIME_DEP, PointInfo&, double*, float, float);
! 	int runge_kutta2(TIME_DIR, TIME_DEP, PointInfo&, double*, float, float);
  	int adapt_step( const VECTOR3& p2, 
  					const VECTOR3& p1, 
  					const VECTOR3& p0, 
! 					const float& minStepsize, 
! 					const float& maxStepsize, 
! 					float* dt, 
  					bool& bAdaptive);
  	
--- 135,146 ----
  	void releaseSeedMemory(void);
  	int euler_cauchy(TIME_DIR, TIME_DEP,float*, float);
! 	int runge_kutta4(TIME_DIR, TIME_DEP, PointInfo&, double*, double, double);
! 	int runge_kutta2(TIME_DIR, TIME_DEP, PointInfo&, double*, double, double);
  	int adapt_step( const VECTOR3& p2, 
  					const VECTOR3& p1, 
  					const VECTOR3& p0, 
! 					const double& minStepsize, 
! 					const double& maxStepsize, 
! 					double* dt, 
  					bool& bAdaptive);
  	

Index: VTFieldLine.cpp
===================================================================
RCS file: /cvsroot/vapor/vapor/lib/flow/VTFieldLine.cpp,v
retrieving revision 1.39
retrieving revision 1.40
diff -C2 -d -r1.39 -r1.40
*** VTFieldLine.cpp	14 Aug 2007 20:02:08 -0000	1.39
--- VTFieldLine.cpp	28 Aug 2007 21:05:43 -0000	1.40
***************
*** 77,81 ****
  int vtCFieldLine::runge_kutta2(TIME_DIR time_dir, TIME_DEP time_dep, 
  							   PointInfo& ci, 
! 							   double* t, float dt, float prevMag)
  {
  	int istat;
--- 77,81 ----
  int vtCFieldLine::runge_kutta2(TIME_DIR time_dir, TIME_DEP time_dep, 
  							   PointInfo& ci, 
! 							   double* t, double dt, double prevMag)
  {
  	int istat;
***************
*** 88,93 ****
  							   PointInfo& ci, 
  							   double* t,			// initial time
! 							   float dt,			//stepsize
! 							   float maxMagDt)		//Largest value of dt*mag(vec)
  {
  	int i, istat;
--- 88,93 ----
  							   PointInfo& ci, 
  							   double* t,			// initial time
! 							   double dt,			//stepsize
! 							   double maxMagDt)		//Largest value of dt*mag(vec)
  {
  	int i, istat;
***************
*** 106,110 ****
  		return OUT_OF_BOUND;
  
! 	if (vel.GetMag()*dt > maxMagDt) return FIELD_TOO_BIG;
  
  	for( i=0; i<3; i++ )
--- 106,112 ----
  		return OUT_OF_BOUND;
  
! 	if (vel.GetDMag()*dt > maxMagDt) {
! 		return FIELD_TOO_BIG;
! 	}
  
  	for( i=0; i<3; i++ )
***************
*** 121,126 ****
  	istat=m_pField->at_phys(fromCell, pt, ci, *t, vel);
  
! 	if (vel.GetMag()*dt > maxMagDt) return FIELD_TOO_BIG;
! 
  	if ( istat!= 1 )
  	{
--- 123,129 ----
  	istat=m_pField->at_phys(fromCell, pt, ci, *t, vel);
  
! 	if (vel.GetDMag()*dt > maxMagDt){ 
! 		return FIELD_TOO_BIG;
! 	}
  	if ( istat!= 1 )
  	{
***************
*** 138,142 ****
  	istat=m_pField->at_phys(fromCell, pt, ci, *t, vel);
  
! 	if (vel.GetMag()*dt > maxMagDt) return FIELD_TOO_BIG;
  
  	if ( istat != 1 )
--- 141,147 ----
  	istat=m_pField->at_phys(fromCell, pt, ci, *t, vel);
  
! 	if (vel.GetDMag()*dt > maxMagDt){
! 		return FIELD_TOO_BIG;
! 	}
  
  	if ( istat != 1 )
***************
*** 156,160 ****
  	fromCell = ci.inCell;
  	istat=m_pField->at_phys(fromCell, pt, ci, *t, vel);
! 	if (vel.GetMag()*dt > maxMagDt) return FIELD_TOO_BIG;
  	if ( istat != 1 )
  	{
--- 161,167 ----
  	fromCell = ci.inCell;
  	istat=m_pField->at_phys(fromCell, pt, ci, *t, vel);
! 	if (vel.GetDMag()*dt > maxMagDt){
! 		return FIELD_TOO_BIG;
! 	}
  	if ( istat != 1 )
  	{
***************
*** 180,186 ****
  							 const VECTOR3& p1, 
  							 const VECTOR3& p0,
! 							 const float& minStepsize, 
! 							 const float& maxStepsize,
! 							 float* dt,
  							 bool& bAdaptive)
  {
--- 187,193 ----
  							 const VECTOR3& p1, 
  							 const VECTOR3& p0,
! 							 const double& minStepsize, 
! 							 const double& maxStepsize,
! 							 double* dt,
  							 bool& bAdaptive)
  {
***************
*** 198,202 ****
  	if(cos_p2p1p0 < m_fLowerAngleAccuracy)
  	{
! 		*dt = (*dt) * (float)0.5;
  		retrace = true;
  		if(*dt < minStepsize)
--- 205,209 ----
  	if(cos_p2p1p0 < m_fLowerAngleAccuracy)
  	{
! 		*dt = (*dt) * 0.5;
  		retrace = true;
  		if(*dt < minStepsize)