Commit 19565270 authored by allens's avatar allens
Browse files

added some comments and moved a bit of code so to be similar across all integrators

git-svn-id: http://visit.ilight.com/svn/visit/trunk/src@23212 18c085ea-50e0-402c-830e-de6fd14e8384
parent ef488c08
......@@ -50,7 +50,7 @@
static const double epsilon = std::numeric_limits<double>::epsilon();
// coefficients for the Adams-Bashforth scheme for up to five steps.
static const double bashforth[ADAMS_BASHFORTH_NSTEPS][ADAMS_BASHFORTH_NSTEPS] =
static const double bashforth[5][5] = //[ADAMS_BASHFORTH_NSTEPS][ADAMS_BASHFORTH_NSTEPS] =
{{ 1., 0., 0., 0., 0. },
{ 3./2., -1./2., 0., 0., 0. },
{ 23./12., -4./3., 5./12., 0., 0. },
......@@ -279,7 +279,7 @@ avtIVPAdamsBashforth::Step(avtIVPField* field, double t_max,
return avtIVPSolver::STEPSIZE_UNDERFLOW;
avtIVPField::Result fieldResult;
avtVector yNew = yCur, vCur;
avtVector yNew, vCur, vNew(0,0,0);
// Get the first vector value for the history.
if ((fieldResult = (*field)(t_local, yCur, vCur)) != avtIVPField::OK)
......@@ -289,9 +289,16 @@ avtIVPAdamsBashforth::Step(avtIVPField* field, double t_max,
// value when calculating pathlines.
history[0] = vCur;
// Calculate the new postion using the Adams-Bashforth algorithm
// Calculate the new velocity using the Adams-Bashforth algorithm
for (size_t i = 0; i < abNSteps; ++i)
yNew += h * bashforth[abStep][i] * history[i];
vNew += bashforth[abStep][i] * history[i];
// Calculate the new position.
yNew = yCur + h * vNew;
// Update the history to save the last N vector values.
for (size_t i = abStep; i>0; --i)
history[i] = history[i-1];
// Increment the number of steps to be taken.
if( abNSteps < ADAMS_BASHFORTH_NSTEPS )
......@@ -300,6 +307,7 @@ avtIVPAdamsBashforth::Step(avtIVPField* field, double t_max,
++abNSteps; // Number of steps to be taken
}
// Convert and save the position.
ivpstep->resize(2);
if( convertToCartesian )
......@@ -315,6 +323,8 @@ avtIVPAdamsBashforth::Step(avtIVPField* field, double t_max,
ivpstep->t0 = t;
ivpstep->t1 = t + h;
// Update for the next step.
numStep++;
yCur = yNew;
......@@ -323,10 +333,6 @@ avtIVPAdamsBashforth::Step(avtIVPField* field, double t_max,
// Reset the step size on sucessful step.
h = h_max;
// Update the history to save the last N vector values.
for (size_t i = abStep; i>0; --i)
history[i] = history[i-1];
return (last ? avtIVPSolver::TERMINATE : avtIVPSolver::OK);
}
......
......@@ -119,6 +119,6 @@ class IVP_API avtIVPAdamsBashforth: public avtIVPSolver
int degenerate_iterations;
double stiffness_eps;
avtVector history[ADAMS_BASHFORTH_NSTEPS];
avtVector dhistory[ADAMS_BASHFORTH_NSTEPS];
// avtVector dhistory[ADAMS_BASHFORTH_NSTEPS];
};
#endif
......@@ -568,6 +568,7 @@ avtIVPDopri5::Step(avtIVPField* field, double t_max,
// make interpolation polynomial
if( ivpstep )
{
// Convert and save the position.
ivpstep->resize(5);
if( convertToCartesian )
......@@ -599,11 +600,15 @@ avtIVPDopri5::Step(avtIVPField* field, double t_max,
// update internal state
// first-same-as-last for k1
k1 = k7;
yCur = y_new;
// Update for the next step.
numStep++;
yCur = y_new;
t = t+h;
// Set the step size on sucessful step.
h = h_new;
numStep++;
return last ? avtIVPSolver::TERMINATE : avtIVPSolver::OK;
}
......@@ -616,6 +621,7 @@ avtIVPDopri5::Step(avtIVPField* field, double t_max,
if( n_accepted >= 1 )
n_rejected++;
// Update the step size to the new step size.
h = h_new;
if (DebugStream::Level5())
......
......@@ -264,9 +264,11 @@ avtIVPEuler::Step(avtIVPField* field, double t_max, avtIVPStep* ivpstep)
{
if ((fieldResult = (*field)(t_local, yCur, vCur)) != avtIVPField::OK)
return ConvertResult(fieldResult);
yNew = yCur + h * vCur; // New position
}
// Convert and save the position.
ivpstep->resize(2);
if( convertToCartesian )
......@@ -282,6 +284,8 @@ avtIVPEuler::Step(avtIVPField* field, double t_max, avtIVPStep* ivpstep)
ivpstep->t0 = t;
ivpstep->t1 = t + h;
// Update for the next step.
numStep++;
yCur = yNew;
......
......@@ -218,20 +218,23 @@ avtIVPLeapfrog::Step(avtIVPField* field, double t_max, avtIVPStep* ivpstep)
avtVector aCur;
if ((fieldResult = (*field)(t_local, yCur, vCur, aCur)) != avtIVPField::OK)
return ConvertResult(fieldResult);
if( numStep )
vNew = vCur + aCur * h; // New velocity
else
vNew = vCur + aCur * h/2.0; // Initial velocity at half step
yNew = yCur + vNew * h; // New position
yNew = yCur + vNew * h; // New position
}
else //if( field->GetOrder() == 1 )
{
if ((fieldResult = (*field)(t_local, yCur, vCur)) != avtIVPField::OK)
return ConvertResult(fieldResult);
yNew = yCur + vCur * h; // New position
}
// Convert and save the position.
ivpstep->resize(2);
if( convertToCartesian )
......@@ -247,6 +250,8 @@ avtIVPLeapfrog::Step(avtIVPField* field, double t_max, avtIVPStep* ivpstep)
ivpstep->t0 = t;
ivpstep->t1 = t + h;
// Update for the next step.
numStep++;
yCur = yNew;
......
......@@ -210,6 +210,7 @@ avtIVPM3DC1Integrator::Step(avtIVPField* field, double t_max,
if( res == avtIVPSolver::OK )
{
// Convert and save the position.
ivpstep->resize( 2 );
// The integration is done in cylindrical coordinates so
......
......@@ -209,9 +209,9 @@ avtIVPRK4::Step(avtIVPField* field, double t_max, avtIVPStep* ivpstep)
if( 0.1*std::abs(h) <= std::abs(t_local)*epsilon )
return avtIVPSolver::STEPSIZE_UNDERFLOW;
//avtIVPSolver::Result res = avtIVPSolverResult::OK;
avtIVPField::Result fieldResult;
// Compute the RK4 values.
avtVector k1;
if ((fieldResult = (*field)(t_local, yCur, k1)) != avtIVPField::OK )
return ConvertResult(fieldResult);
......@@ -228,8 +228,10 @@ avtIVPRK4::Step(avtIVPField* field, double t_max, avtIVPStep* ivpstep)
if ((fieldResult = (*field)(t_local+h, yCur+h*k3, k4)) != avtIVPField::OK)
return ConvertResult(fieldResult);
// Calculate the new position.
avtVector yNew = yCur + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
// Convert and save the position.
ivpstep->resize(2);
if( convertToCartesian )
......@@ -245,6 +247,8 @@ avtIVPRK4::Step(avtIVPField* field, double t_max, avtIVPStep* ivpstep)
ivpstep->t0 = t;
ivpstep->t1 = t + h;
// Update for the next step.
numStep++;
yCur = yNew;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment