@@ -231,6 +231,7 @@ namespace aspect
231231 Py_DECREF (pArgs);
232232
233233 const ArrayView<const double > data = PythonHelper::numpy_to_array_view (pValue);
234+ const double one_over_dt = 1.0 / ((this ->get_timestep () > 0.0 ) ? this ->get_timestep () : 1.0 );
234235 for (size_t i=0 ; i<data.size (); ++i)
235236 {
236237 const Tensor<1 ,dim> gravity = this ->get_gravity_model ().gravity_vector (this ->evaluation_points [i]);
@@ -239,7 +240,7 @@ namespace aspect
239240 if (gravity.norm () > 0.0 )
240241 topography_direction = -gravity / gravity.norm ();
241242
242- velocities[i] = topography_direction * data[i];
243+ velocities[i] = topography_direction * data[i] * one_over_dt ;
243244 }
244245
245246 Py_DECREF (pValue);
@@ -254,12 +255,13 @@ namespace aspect
254255 //const auto &mapping = this->get_mapping();
255256 std::vector<Point<dim>> real_evaluation_points(this->evaluation_points.size());
256257 std::vector<std::vector<double>> data(this->evaluation_points.size(), std::vector<double>(dim, 0.0));
258+ const double one_over_dt = 1.0 / ((this->get_timestep() > 0.0) ? this->get_timestep() : 1.0);
257259
258260 for (unsigned int i=0; i<this->evaluation_points.size(); ++i)
259261 {
260262 real_evaluation_points[i] = this->evaluation_points[i]; // TODO: use mapping to compute real position
261263 for (unsigned int c=0; c<dim; ++c)
262- data[i][c] = velocities[i][c];
264+ data[i][c] = velocities[i][c] * one_over_dt ;
263265 }
264266
265267 const std::vector<std::string> data_component_names(dim, "velocity");
0 commit comments