Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
CLHEP
CLHEP
Commits
c4358e8a
Commit
c4358e8a
authored
Aug 06, 2013
by
Joseph Boudreau
Browse files
bugfix. adaptive stepper was not adaptive.
parent
96cce689
Changes
1
Hide whitespace changes
Inline
Side-by-side
GenericFunctions/src/AdaptiveRKStepper.cc
View file @
c4358e8a
...
...
@@ -38,9 +38,9 @@ namespace Genfun {
if
(
s
.
time
==
0.0
)
{
stepsize
=
sStepsize
;
}
const
unsigned
int
p
=
eeStepper
->
order
();
// Order of the stepper
const
unsigned
int
p
=
eeStepper
->
order
();
// Order of the stepper
const
double
deltaMax
=
T
*
std
::
pow
(
S
/
Rmax
,
(
int
)(
p
+
1
));
// Maximum error 4 adjustment.
const
double
TINY
=
1.0E-30
;
// Denominator regularization
const
double
TINY
=
1.0E-30
;
// Denominator regularization
double
hnext
;
//
// Time limited step ?
...
...
@@ -66,16 +66,18 @@ namespace Genfun {
// Bail out and try a smaller step.
//
h
=
std
::
max
(
S
*
h
*
std
::
pow
(
T
/
(
delta
+
TINY
),
1.0
/
(
p
+
1
)),
Rmin
*
h
);
if
(
!
(((
float
)
s
.
time
+
h
-
(
float
)
s
.
time
)
>
0
)
)
{
std
::
cerr
<<
"Warning, RK Integrator step underflow"
<<
std
::
endl
;
if
(
!
(((
double
)
(
s
.
time
+
h
)
-
(
double
)
s
.
time
)
>
0
)
)
{
throw
std
::
runtime_error
(
"Warning, RK Integrator step underflow"
)
;
}
d
.
time
=
s
.
time
+
h
;
hnext
=
h
;
continue
;
}
else
{
if
(
delta
>
deltaMax
)
{
if
(
delta
<
deltaMax
)
{
hnext
=
S
*
h
*
std
::
pow
(
T
/
(
delta
+
TINY
),
1.0
/
(
p
+
1
));
// stepsize is supposed to increase;
if
(
hnext
<
h
)
hnext
=
h
;
}
else
{
hnext
=
Rmax
*
h
;
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment