[Debian-astro-commits] [gyoto] 102/221: Use the integrator set in state_ when refining geodesic in Worldline::getCoord()

[Thibaut Jean-Claude Paumard]

This is an automated email from the git hooks/post-receive script.

thibaut pushed a commit to branch master
in repository gyoto.

commit 6f7f809a06a58b304d99decacd62d5c2570949c3
Author: Thibaut Paumard <paumard at users.sourceforge.net>
Date:   Fri Nov 21 14:50:11 2014 +0100

    Use the integrator set in state_ when refining geodesic in Worldline::getCoord()
---
 include/GyotoWorldline.h  | 25 ++++++++++++++++++++-----
 lib/Worldline.C           |  4 ++--
 lib/WorldlineIntegState.C | 20 ++++++++++++++++++++
 yorick/check-kerrbl.i     |  2 +-
 4 files changed, 43 insertions(+), 8 deletions(-)

diff --git a/include/GyotoWorldline.h b/include/GyotoWorldline.h
index e8a4b92..6044d9e 100644
--- a/include/GyotoWorldline.h
+++ b/include/GyotoWorldline.h
@@ -636,6 +636,19 @@ class Gyoto::Worldline::IntegState::Generic : SmartPointee {
    */
   virtual int nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0;
 
+  /// Make one step of exactly this size.
+  /**
+   * doStep() is meant to refine a computation made using
+   * nextStep(). In particular, there is no checking for norm
+   * conservation.
+   *
+   * \param[in] coordin current position-velocity;
+   * \param[in] step  exact step to use.
+   * \param[out] coordout next position-velocity;
+   */
+  virtual void doStep(double const coordin[8], 
+		      double step,
+		      double coordout[8])=0;
 };
 
 /**
@@ -664,13 +677,12 @@ class Gyoto::Worldline::IntegState::Legacy : public Generic {
   void init(Worldline * line, const double *coord, const double delta);
   virtual std::string kind();
 
-  /// Make one step.
-  /**
-   * \param[out] coord Next position-velocity;
-   * \param[in] h1max maximum step in case of adaptive integration
-   */
   virtual int nextStep(double *coord, double h1max=1e6);
 
+  virtual void doStep(double const coordin[8], 
+		      double step,
+		      double coordout[8]);
+
   virtual ~Legacy();
 };
 
@@ -726,6 +738,9 @@ class Gyoto::Worldline::IntegState::Boost : public Generic {
   virtual void init();
   virtual void init(Worldline * line, const double *coord, const double delta);
   virtual int nextStep(double *coord, double h1max=1e6);
+  virtual void doStep(double const coordin[8], 
+		      double step,
+		      double coordout[8]);
   virtual std::string kind();
   
 };
diff --git a/lib/Worldline.C b/lib/Worldline.C
index 83bbb7c..ff21efa 100644
--- a/lib/Worldline.C
+++ b/lib/Worldline.C
@@ -785,7 +785,7 @@ void Worldline::getCoord(double const * const dates, size_t const n_dates,
     bestl[5] = x1dot_[curl];
     bestl[6] = x2dot_[curl];
     bestl[7] = x3dot_[curl];
-    metric_ -> myrk4(this, bestl, dtaul, resl);
+    state_ -> doStep(bestl, dtaul, resl);
 
     // from above...
     besth[0] =    x0_[curh];
@@ -796,7 +796,7 @@ void Worldline::getCoord(double const * const dates, size_t const n_dates,
     besth[5] = x1dot_[curh];
     besth[6] = x2dot_[curh];
     besth[7] = x3dot_[curh];
-    metric_ -> myrk4(this, besth, dtauh, resh);
+    state_ -> doStep(besth, dtauh, resh);
 
 #   if GYOTO_DEBUG_ENABLED
     GYOTO_IF_DEBUG
diff --git a/lib/WorldlineIntegState.C b/lib/WorldlineIntegState.C
index 8ae18ac..624ed0c 100644
--- a/lib/WorldlineIntegState.C
+++ b/lib/WorldlineIntegState.C
@@ -146,6 +146,11 @@ int Worldline::IntegState::Legacy::nextStep(double coord[8], double h1max) {
   return 0;
 }
 
+void Worldline::IntegState::Legacy::doStep(double const coordin[8],
+					   double step, 
+					   double coordout[8]) {
+  gg_ -> myrk4(line_, coordin, step, coordout); 
+}
 
 std::string Worldline::IntegState::Legacy::kind() { return "Legacy"; } 
 
@@ -273,6 +278,21 @@ int Worldline::IntegState::Boost::nextStep(double coord[8], double h1max) {
   return line_->stopcond;
 }
 
+void Worldline::IntegState::Boost::doStep(double const coordin[8],
+					 double step, 
+					 double coordout[8]) {
+  // We first make a C++ std::array out of the bare C array:
+  std::array<double, 8> inout = {
+    coordin[0], coordin[1], coordin[2], coordin[3],
+    coordin[4], coordin[5], coordin[6], coordin[7]};
+
+  // We call the Boost stepper
+  do_step_(inout, step);
+
+  // Copy the result
+  for (size_t i=0; i<=7; ++i) coordout[i]=inout[i];
+}
+
 std::string Worldline::IntegState::Boost::kind() {
   if (kind_== Kind::runge_kutta_cash_karp54) return "runge_kutta_cash_karp54";
   if (kind_== Kind::runge_kutta_fehlberg78) return "runge_kutta_fehlberg78";
diff --git a/yorick/check-kerrbl.i b/yorick/check-kerrbl.i
index ca78e39..f92032f 100644
--- a/yorick/check-kerrbl.i
+++ b/yorick/check-kerrbl.i
@@ -124,7 +124,7 @@ if (gyoto_haveBoost()) {
   //coords7=st7(get_coord=dates);
   if (max(abs(coords-coords3))>1e-3)
     error, "two integrators don't yield the same result";
-  if (max(abs(coords-coords5))>1e-3)
+  if (max(abs(coords-coords5))>3e-3)
     error, "two integrators don't yield the same result";
   if (max(abs(coords-coords4))>1e-3)
     error, "two integrators don't yield the same result";

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-astro/packages/gyoto.git