]> https://gitweb.dealii.org/ - dealii.git/commitdiff
SUNDIALS: Use N_Vector module in ARKode wrapper 11655/head
authorSebastian Proell <proell@lnm.mw.tum.de>
Sat, 30 Jan 2021 11:23:31 +0000 (12:23 +0100)
committerSebastian Proell <proell@lnm.mw.tum.de>
Sun, 31 Jan 2021 19:21:54 +0000 (20:21 +0100)
include/deal.II/sundials/arkode.h
source/sundials/arkode.cc
tests/sundials/arkode_03.with_sundials.geq.5.4.0.output

index 5f565e36405bf9120fec60da8476d4ae057b2f67..bd5bda0dcf80850fd6d83a24b7a4f5ca55a68236 100644 (file)
@@ -41,6 +41,8 @@
 #  ifdef DEAL_II_WITH_MPI
 #    include <nvector/nvector_parallel.h>
 #  endif
+#  include <deal.II/sundials/n_vector.h>
+
 #  include <boost/signals2.hpp>
 
 #  include <sundials/sundials_linearsolver.h>
@@ -610,26 +612,7 @@ namespace SUNDIALS
      * @param[in,out] y   The new initial solution
      */
     void
-    reset(const double t, const double h, const VectorType &y);
-
-    /*!
-     * Create a new SUNDIALS vector from a given template.
-     *
-     * @note Vectors created this way should be freed with free_vector().
-     *
-     * @param template_vector The vector to use as a template for the layout of
-     *   a new vector.
-     */
-    N_Vector
-    create_vector(const VectorType &template_vector) const;
-
-    /**
-     * Free a SUNDIALS vector created with create_vector().
-     *
-     * @param vector The vector to free
-     */
-    void
-    free_vector(N_Vector vector) const;
+    reset(const double t, const double h, VectorType &y);
 
     /**
      * Provides user access to the internally used ARKODE memory.
@@ -1299,7 +1282,7 @@ namespace SUNDIALS
     /**
      * Set up the (non)linear solver and preconditioners in the ARKODE memory
      * object based on the user-specified functions.
-     * @param solution The solution vector whihc is used as a template to create
+     * @param solution The solution vector which is used as a template to create
      *   new vectors.
      */
     void
@@ -1335,12 +1318,12 @@ namespace SUNDIALS
     /**
      * ARKode solution vector.
      */
-    N_Vector yy;
+    internal::NVectorView<VectorType> yy;
 
     /**
      * ARKode absolute tolerances vector.
      */
-    N_Vector abs_tolls;
+    internal::NVectorView<VectorType> abs_tolls;
 
     /**
      * MPI communicator. SUNDIALS solver runs happily in
index 0a52dfeced16155544c1e7dc5a13bdaf432adb7f..7cf6d94a4d39c8cc09d7e52e3d168b14f5a7ca8a 100644 (file)
@@ -33,7 +33,7 @@
 #    include <deal.II/lac/petsc_vector.h>
 #  endif
 
-#  include <deal.II/sundials/copy.h>
+#  include <deal.II/sundials/n_vector.h>
 
 #  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
 #    include <arkode/arkode_impl.h>
@@ -72,21 +72,11 @@ namespace SUNDIALS
       Assert(user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
-      GrowingVectorMemory<VectorType> mem;
-
-      typename VectorMemory<VectorType>::Pointer src_yy(mem);
-      solver.reinit_vector(*src_yy);
-
-      typename VectorMemory<VectorType>::Pointer dst_yp(mem);
-      solver.reinit_vector(*dst_yp);
-
-      copy(*src_yy, yy);
 
-      int err = solver.explicit_function(tt, *src_yy, *dst_yp);
+      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+      auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
 
-      copy(yp, *dst_yp);
-
-      return err;
+      return solver.explicit_function(tt, *src_yy, *dst_yp);
     }
 
 
@@ -101,21 +91,11 @@ namespace SUNDIALS
       Assert(user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
-      GrowingVectorMemory<VectorType> mem;
 
-      typename VectorMemory<VectorType>::Pointer src_yy(mem);
-      solver.reinit_vector(*src_yy);
+      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+      auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
 
-      typename VectorMemory<VectorType>::Pointer dst_yp(mem);
-      solver.reinit_vector(*dst_yp);
-
-      copy(*src_yy, yy);
-
-      int err = solver.implicit_function(tt, *src_yy, *dst_yp);
-
-      copy(yp, *dst_yp);
-
-      return err;
+      return solver.implicit_function(tt, *src_yy, *dst_yp);
     }
 
 
@@ -135,17 +115,9 @@ namespace SUNDIALS
       Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-      GrowingVectorMemory<VectorType> mem;
-
-      typename VectorMemory<VectorType>::Pointer src_ypred(mem);
-      solver.reinit_vector(*src_ypred);
-
-      typename VectorMemory<VectorType>::Pointer src_fpred(mem);
-      solver.reinit_vector(*src_fpred);
-
-      copy(*src_ypred, ypred);
-      copy(*src_fpred, fpred);
 
+      auto *src_ypred = internal::unwrap_nvector_const<VectorType>(ypred);
+      auto *src_fpred = internal::unwrap_nvector_const<VectorType>(fpred);
       // avoid reinterpret_cast
       bool jcurPtr_tmp = false;
       int  err         = solver.setup_jacobian(convfail,
@@ -178,31 +150,21 @@ namespace SUNDIALS
       Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-      GrowingVectorMemory<VectorType> mem;
-
-      typename VectorMemory<VectorType>::Pointer src(mem);
-      solver.reinit_vector(*src);
 
-      typename VectorMemory<VectorType>::Pointer src_ycur(mem);
-      solver.reinit_vector(*src_ycur);
+      auto *dst      = internal::unwrap_nvector<VectorType>(b);
+      auto *src_ycur = internal::unwrap_nvector_const<VectorType>(ycur);
+      auto *src_fcur = internal::unwrap_nvector_const<VectorType>(fcur);
 
-      typename VectorMemory<VectorType>::Pointer src_fcur(mem);
-      solver.reinit_vector(*src_fcur);
-
-      typename VectorMemory<VectorType>::Pointer dst(mem);
-      solver.reinit_vector(*dst);
-
-      copy(*src, b);
-      copy(*src_ycur, ycur);
-      copy(*src_fcur, fcur);
+      // make a temporary copy to work on in the user call
+      VectorType src = *dst;
 
       int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
                                              arkode_mem->ark_gamma,
                                              *src_ycur,
                                              *src_fcur,
-                                             *src,
+                                             src,
                                              *dst);
-      copy(b, *dst);
+
 
       return err;
     }
@@ -236,18 +198,13 @@ namespace SUNDIALS
       Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-      GrowingVectorMemory<VectorType> mem;
-
-      typename VectorMemory<VectorType>::Pointer src(mem);
-      solver.reinit_vector(*src);
 
-      typename VectorMemory<VectorType>::Pointer dst(mem);
-      solver.reinit_vector(*dst);
+      auto *dst = internal::unwrap_nvector<VectorType>(b);
 
-      copy(*src, b);
+      // make a temporary copy to work on in the user call
+      VectorType src = *dst;
 
-      int err = solver.solve_mass_system(*src, *dst);
-      copy(b, *dst);
+      int err = solver.solve_mass_system(src, *dst);
 
       return err;
     }
@@ -268,22 +225,13 @@ namespace SUNDIALS
         *static_cast<ARKode<VectorType> *>(user_data);
       GrowingVectorMemory<VectorType> mem;
 
-      typename VectorMemory<VectorType>::Pointer src_v(mem);
-      solver.reinit_vector(*src_v);
-      typename VectorMemory<VectorType>::Pointer dst_Jv(mem);
-      solver.reinit_vector(*dst_Jv);
-      typename VectorMemory<VectorType>::Pointer src_y(mem);
-      solver.reinit_vector(*src_y);
-      typename VectorMemory<VectorType>::Pointer src_fy(mem);
-      solver.reinit_vector(*src_fy);
-      copy(*src_v, v);
-      copy(*src_y, y);
-      copy(*src_fy, fy);
-
-      int err =
-        solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy);
-      copy(Jv, *dst_Jv);
-      return err;
+      auto *src_v  = internal::unwrap_nvector_const<VectorType>(v);
+      auto *src_y  = internal::unwrap_nvector_const<VectorType>(y);
+      auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+
+      auto *dst_Jv = internal::unwrap_nvector<VectorType>(Jv);
+
+      return solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy);
     }
 
 
@@ -298,15 +246,11 @@ namespace SUNDIALS
       Assert(user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer src_y(mem);
-      solver.reinit_vector(*src_y);
-      typename VectorMemory<VectorType>::Pointer src_fy(mem);
-      solver.reinit_vector(*src_fy);
-      copy(*src_y, y);
-      copy(*src_fy, fy);
-      int err = solver.jacobian_times_setup(t, *src_y, *src_fy);
-      return err;
+
+      auto *src_y  = internal::unwrap_nvector_const<VectorType>(y);
+      auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+
+      return solver.jacobian_times_setup(t, *src_y, *src_fy);
     }
 
 
@@ -327,23 +271,14 @@ namespace SUNDIALS
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
 
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer src_y(mem);
-      solver.reinit_vector(*src_y);
-      typename VectorMemory<VectorType>::Pointer src_fy(mem);
-      solver.reinit_vector(*src_fy);
-      typename VectorMemory<VectorType>::Pointer src_r(mem);
-      solver.reinit_vector(*src_r);
-      typename VectorMemory<VectorType>::Pointer dst_z(mem);
-      solver.reinit_vector(*dst_z);
-      copy(*src_y, y);
-      copy(*src_fy, fy);
-      copy(*src_r, r);
-
-      int status = solver.jacobian_preconditioner_solve(
+      auto *src_y  = internal::unwrap_nvector_const<VectorType>(y);
+      auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+      auto *src_r  = internal::unwrap_nvector_const<VectorType>(r);
+
+      auto *dst_z = internal::unwrap_nvector<VectorType>(z);
+
+      return solver.jacobian_preconditioner_solve(
         t, *src_y, *src_fy, *src_r, *dst_z, gamma, delta, lr);
-      copy(z, *dst_z);
-      return status;
     }
 
 
@@ -361,13 +296,10 @@ namespace SUNDIALS
       Assert(user_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer src_y(mem);
-      solver.reinit_vector(*src_y);
-      typename VectorMemory<VectorType>::Pointer src_fy(mem);
-      solver.reinit_vector(*src_fy);
-      copy(*src_y, y);
-      copy(*src_fy, fy);
+
+      auto *src_y  = internal::unwrap_nvector_const<VectorType>(y);
+      auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+
       return solver.jacobian_preconditioner_setup(
         t, *src_y, *src_fy, jok, *jcurPtr, gamma);
     }
@@ -384,17 +316,11 @@ namespace SUNDIALS
       Assert(mtimes_data != nullptr, ExcInternalError());
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(mtimes_data);
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer src_v(mem);
-      solver.reinit_vector(*src_v);
-      typename VectorMemory<VectorType>::Pointer dst_Mv(mem);
-      solver.reinit_vector(*dst_Mv);
 
-      copy(*src_v, v);
-      int err = solver.mass_times_vector(t, *src_v, *dst_Mv);
-      copy(Mv, *dst_Mv);
+      auto *src_v  = internal::unwrap_nvector_const<VectorType>(v);
+      auto *dst_Mv = internal::unwrap_nvector<VectorType>(Mv);
 
-      return err;
+      return solver.mass_times_vector(t, *src_v, *dst_Mv);
     }
 
 
@@ -425,17 +351,10 @@ namespace SUNDIALS
       ARKode<VectorType> &solver =
         *static_cast<ARKode<VectorType> *>(user_data);
 
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer src_r(mem);
-      solver.reinit_vector(*src_r);
-      typename VectorMemory<VectorType>::Pointer dst_z(mem);
-      solver.reinit_vector(*dst_z);
-      copy(*src_r, r);
-
-      int status =
-        solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr);
-      copy(z, *dst_z);
-      return status;
+      auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
+      auto *dst_z = internal::unwrap_nvector<VectorType>(z);
+
+      return solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr);
     }
 
 
@@ -506,23 +425,17 @@ namespace SUNDIALS
       auto                content = access_content<VectorType>(LS);
       ARKode<VectorType> &solver  = *(content->solver);
 
-      GrowingVectorMemory<VectorType>            mem;
-      typename VectorMemory<VectorType>::Pointer rhs(mem);
-      solver.reinit_vector(*rhs);
-      typename VectorMemory<VectorType>::Pointer dst(mem);
-      solver.reinit_vector(*dst);
+      auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
+      auto *dst_x = internal::unwrap_nvector<VectorType>(x);
 
-      copy(*rhs, b);
-
-      SundialsOperator<VectorType>       op(solver,
+      SundialsOperator<VectorType> op(solver,
                                       content->A_data,
                                       content->a_times_fn);
+
       SundialsPreconditioner<VectorType> preconditioner(
         solver, content->P_data, content->preconditioner_solve, tol);
 
-      int err = content->lsolve(op, preconditioner, *dst, *rhs, tol);
-      copy(x, *dst);
-      return err;
+      return content->lsolve(op, preconditioner, *dst_x, *src_b, tol);
     }
 
 
@@ -630,8 +543,6 @@ namespace SUNDIALS
                              const MPI_Comm &      mpi_comm)
     : data(data)
     , arkode_mem(nullptr)
-    , yy(nullptr)
-    , abs_tolls(nullptr)
     , communicator(is_serial_vector<VectorType>::value ?
                      MPI_COMM_SELF :
                      Utilities::MPI::duplicate_communicator(mpi_comm))
@@ -665,30 +576,6 @@ namespace SUNDIALS
   unsigned int
   ARKode<VectorType>::solve_ode(VectorType &solution)
   {
-    const unsigned int system_size = solution.size();
-
-    // The solution is stored in solution. Here we take only a view of it.
-#  ifdef DEAL_II_WITH_MPI
-    if (is_serial_vector<VectorType>::value == false)
-      {
-        const IndexSet    is                = solution.locally_owned_elements();
-        const std::size_t local_system_size = is.n_elements();
-
-        yy = N_VNew_Parallel(communicator, local_system_size, system_size);
-
-        abs_tolls =
-          N_VNew_Parallel(communicator, local_system_size, system_size);
-      }
-    else
-#  endif
-      {
-        Assert(is_serial_vector<VectorType>::value,
-               ExcInternalError(
-                 "Trying to use a serial code with a parallel vector."));
-        yy        = N_VNew_Serial(system_size);
-        abs_tolls = N_VNew_Serial(system_size);
-      }
-
     DiscreteTime time(data.initial_time,
                       data.final_time,
                       data.initial_step_size);
@@ -711,8 +598,6 @@ namespace SUNDIALS
         (void)status;
         AssertARKode(status);
 
-        copy(solution, yy);
-
         time.set_next_step_size(actual_next_time - time.get_current_time());
         time.advance_time();
 
@@ -727,70 +612,31 @@ namespace SUNDIALS
                       time.get_step_number());
       }
 
-      // Free the vectors which are no longer used.
-#  ifdef DEAL_II_WITH_MPI
-    if (is_serial_vector<VectorType>::value == false)
-      {
-        N_VDestroy_Parallel(yy);
-        N_VDestroy_Parallel(abs_tolls);
-      }
-    else
-#  endif
-      {
-        N_VDestroy_Serial(yy);
-        N_VDestroy_Serial(abs_tolls);
-      }
-
     return time.get_step_number();
   }
 
 #  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
   template <typename VectorType>
   void
-  ARKode<VectorType>::reset(const double      current_time,
-                            const double      current_time_step,
-                            const VectorType &solution)
+  ARKode<VectorType>::reset(const double current_time,
+                            const double current_time_step,
+                            VectorType & solution)
   {
-    unsigned int system_size;
-
     if (arkode_mem)
       ARKodeFree(&arkode_mem);
 
     arkode_mem = ARKodeCreate();
 
-    // Free the vectors which are no longer used.
-    if (yy)
-      {
-        free_vector(yy);
-        free_vector(abs_tolls);
-      }
-
     int status;
     (void)status;
-    system_size = solution.size();
-#    ifdef DEAL_II_WITH_MPI
-    if (is_serial_vector<VectorType>::value == false)
-      {
-        const IndexSet    is                = solution.locally_owned_elements();
-        const std::size_t local_system_size = is.n_elements();
-
-        yy = N_VNew_Parallel(communicator, local_system_size, system_size);
-
-        abs_tolls =
-          N_VNew_Parallel(communicator, local_system_size, system_size);
-      }
-    else
-#    endif
-      {
-        yy        = N_VNew_Serial(system_size);
-        abs_tolls = N_VNew_Serial(system_size);
-      }
-
-    copy(yy, solution);
 
     Assert(explicit_function || implicit_function,
            ExcFunctionNotProvided("explicit_function || implicit_function"));
 
+    // just a view on the memory in solution, all write operations on yy by
+    // ARKODE will automatically be mirrored to solution
+    yy = internal::make_nvector_view(solution);
+
     status = ARKodeInit(
       arkode_mem,
       explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
@@ -801,7 +647,7 @@ namespace SUNDIALS
 
     if (get_local_tolerances)
       {
-        copy(abs_tolls, get_local_tolerances());
+        abs_tolls = make_nvector_view(get_local_tolerances());
         status =
           ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
         AssertARKode(status);
@@ -885,23 +731,17 @@ namespace SUNDIALS
   void
   ARKode<VectorType>::reset(const double current_time,
                             const double current_time_step,
-                            const VectorType &solution)
+                            VectorType &solution)
   {
     if (arkode_mem)
       ARKStepFree(&arkode_mem);
 
-    // Free the vectors which are no longer used.
-    if (yy)
-      {
-        free_vector(yy);
-        free_vector(abs_tolls);
-      }
-
     int status;
     (void)status;
-    yy = create_vector(solution);
-    abs_tolls = create_vector(solution);
-    copy(yy, solution);
+
+    // just a view on the memory in solution, all write operations on yy by
+    // ARKODE will automatically be mirrored to solution
+    yy = internal::make_nvector_view(solution);
 
     Assert(explicit_function || implicit_function,
            ExcFunctionNotProvided("explicit_function || implicit_function"));
@@ -912,9 +752,11 @@ namespace SUNDIALS
       current_time,
       yy);
 
+    Assert(arkode_mem != nullptr, ExcInternalError());
+
     if (get_local_tolerances)
       {
-        copy(abs_tolls, get_local_tolerances());
+        abs_tolls = internal::make_nvector_view(get_local_tolerances());
         status =
           ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
         AssertARKode(status);
@@ -1009,13 +851,12 @@ namespace SUNDIALS
       }
     else
       {
-        N_Vector y_template = create_vector(solution);
+        auto y_template = internal::make_nvector_view(solution);
 
         SUNNonlinearSolver fixed_point_solver =
           SUNNonlinSol_FixedPoint(y_template,
                                   data.anderson_acceleration_subspace);
 
-        free_vector(y_template);
         status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
         AssertARKode(status);
       }
@@ -1098,43 +939,6 @@ namespace SUNDIALS
 
 
 
-  template <typename VectorType>
-  N_Vector
-  ARKode<VectorType>::create_vector(const VectorType &template_vector) const
-  {
-    N_Vector new_vector;
-    auto     system_size = template_vector.size();
-
-#  ifdef DEAL_II_WITH_MPI
-    if (is_serial_vector<VectorType>::value == false)
-      {
-        new_vector =
-          N_VNew_Parallel(communicator,
-                          template_vector.locally_owned_elements().n_elements(),
-                          system_size);
-      }
-    else
-#  endif
-      {
-        new_vector = N_VNew_Serial(system_size);
-      }
-    return new_vector;
-  }
-
-  template <typename VectorType>
-  void
-  ARKode<VectorType>::free_vector(N_Vector vector) const
-  {
-#  ifdef DEAL_II_WITH_MPI
-    if (is_serial_vector<VectorType>::value == false)
-      N_VDestroy_Parallel(vector);
-    else
-#  endif
-      N_VDestroy_Serial(vector);
-  }
-
-
-
   template <typename VectorType>
   void *
   ARKode<VectorType>::get_arkode_memory() const
@@ -1163,17 +967,11 @@ namespace SUNDIALS
   SundialsOperator<VectorType>::vmult(VectorType &      dst,
                                       const VectorType &src) const
   {
-    N_Vector sun_dst = solver.create_vector(dst);
-    N_Vector sun_src = solver.create_vector(src);
-    copy(sun_src, src);
-    int status = a_times_fn(A_data, sun_src, sun_dst);
+    auto sun_dst = internal::make_nvector_view(dst);
+    auto sun_src = internal::make_nvector_view(src);
+    int  status  = a_times_fn(A_data, sun_src, sun_dst);
     (void)status;
     AssertARKode(status);
-
-    copy(dst, sun_dst);
-
-    solver.free_vector(sun_dst);
-    solver.free_vector(sun_src);
   }
 
 
@@ -1204,20 +1002,14 @@ namespace SUNDIALS
         return;
       }
 
-    N_Vector sun_dst = solver.create_vector(dst);
-    N_Vector sun_src = solver.create_vector(src);
-    copy(sun_src, src);
+    auto sun_dst = internal::make_nvector_view(dst);
+    auto sun_src = internal::make_nvector_view(src);
     // for custom preconditioners no distinction between left and right
     // preconditioning is made
     int status =
       p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
     (void)status;
     AssertARKode(status);
-
-    copy(dst, sun_dst);
-
-    solver.free_vector(sun_dst);
-    solver.free_vector(sun_src);
   }
 #  endif
 
index 9e2ed6e8a6bb0529664faaf7948351155a4a96f2..503f3a999d0607944d2d78b05296a52b7119c05e 100644 (file)
@@ -1,12 +1,12 @@
 
-DEAL::0 3.9 1.1 2.8
-DEAL::1 2.075274578 1.059601971 2.499947283
-DEAL::2 1.100460018 1.723610111 2.499972388
-DEAL::3 0.7896202372 2.326668043 2.499979776
-DEAL::4 0.8117726433 2.725583965 2.499987084
-DEAL::5 1.187002916 2.595541878 2.499973736
-DEAL::6 1.887311324 1.494098057 2.499952698
-DEAL::7 1.323271341 1.619016372 2.499965777
-DEAL::8 0.9217874641 2.130219386 2.499970637
-DEAL::9 0.8495425771 2.539475069 2.499991366
-DEAL::10 1.064969831 2.595955794 2.499973383
\ No newline at end of file
+DEAL::0.00000 3.900000000 1.100000000 2.800000000
+DEAL::1.000000000 2.075274578 1.059601971 2.499947291
+DEAL::2.000000000 1.100460018 1.723610111 2.499972248
+DEAL::3.000000000 0.7896202372 2.326668043 2.499980243
+DEAL::4.000000000 0.8117726433 2.725583965 2.499982784
+DEAL::5.000000000 1.187002916 2.595541878 2.499973234
+DEAL::6.000000000 1.887311324 1.494098057 2.499934582
+DEAL::7.000000000 1.323271340 1.619016372 2.499963670
+DEAL::8.000000000 0.9217874641 2.130219386 2.499972986
+DEAL::9.000000000 0.8495425770 2.539475069 2.499983024
+DEAL::10.00000000 1.064969831 2.595955794 2.499973386

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.