### C++ Runge-Kutta-Method l1-norm returns nan

122
views
0
6 months ago by

Hi,
I tried to implement a runge-kutta-method to get the optimal stepsize for a time-diskrete problem.
I still got an implementation working for a fixed stepsize. My implementation of runge-Kutta is:

double rungeKuttaFifthOrder(std::shared_ptr<dolfin::Function> y,
double max_stepsize,
double tol,
double start_time,
std::shared_ptr<convectionDiffusion3D::Form_L> L,
std::shared_ptr<convectionDiffusion3D::Form_a> a)
{
const unsigned int N = 6;

std::array<std::shared_ptr<dolfin::FunctionAXPY>, N> k;
//initialize factor arrays
std::array<std::array<double, 5>, N> factors = {
std::array<double, 5>{0.0, 0.0, 0.0, 0.0, 0.0},
std::array<double, 5>{0.25, 0.0, 0.0, 0.0, 0.0},
std::array<double, 5>{3.0 / 32.0, 9.0 / 32.0, 0.0, 0.0, 0.0},
std::array<double, 5>{
1932.0 / 2197.0, -7200.0 / 2197.0, 7296.0 / 2197.0, 0.0, 0.0},
std::array<double, 5>{
439.0 / 216.0, 8.0, 3680.0 / 513.0, 845.0 / 4104.0, 0.0},
std::array<double, 5>{
-8.0 / 27.0, 2.0, -3544.0 / 2565.0, 1859.0 / 4104, -11.0 / 40.0}};
std::array<double, N> eval_factors = {
0, 1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 0.5};
std::array<double, N> y_factors = {
25.0 / 216.0, 0.0, 1408.0 / 2565.0, 2197.0 / 4104.0, -1.0 / 5.0, 0.0};
std::array<double, N> z_factors = {16.0 / 135.0,
0.0,
6656.0 / 12825.0,
28561.0 / 56430.0,
-9.0 / 50.0,
2.0 / 55.0};
//create FunctionAXPY objects with same content as y
auto z_next = dolfin::FunctionAXPY(y, 1.0);
auto y_next = dolfin::FunctionAXPY(y, 1.0);

//compute functions to be evaluated and evaluate them at given time
for (unsigned int i = 0; i < N; i++) {
auto ev_func_axpy = dolfin::FunctionAXPY(y, 1.0);
if (i > 0) {
//compute function
for (unsigned int j = 0; j < i - 1; j++) {
if (j > 0)
ev_func_axpy =
ev_func_axpy + (*k.at(j) * factors.at(i).at(j));
}
}
//create dolfin::Function from dolfin::FunctionAXPY
std::shared_ptr<dolfin::Function> ev_func =
std::make_shared<dolfin::Function>(y->function_space());
*ev_func = ev_func_axpy;

//compute evaluation time
auto eval_offset = max_stepsize * eval_factors.at(i);
auto eval_time = start_time + eval_offset;

//evaluate Function and store value in k[i]
k.at(i) = std::make_shared<dolfin::FunctionAXPY>(
evaluate_func(ev_func, start_time, eval_time, L, a), 1.0);

//compute y_next and z_next
y_next = y_next + (*k.at(i) * y_factors.at(i));
z_next = z_next + (*k.at(i) * z_factors.at(i));
}

auto diff = z_next - y_next;
auto diff_func = dolfin::Function(y->function_space());
diff_func = diff;
//compute l1-norm (returns NaN)
double norm_diff = dolfin::norm(*diff_func.vector(), std::string("l1"));
double s = tol / (2.0 * norm_diff);
s = pow(s, 0.25);

return s;

Unfortunately the line

double norm_diff = dolfin::norm(*diff_func.vector(), std::string("l1));

returns "NaN" and therefore the whole function returns "NaN". The
evaluate_func function solves the equation-system at a given time and returns the function as shared_ptr.
I was not able to verify that the steps which were done before the computation of the l1-norm were correct because the debugger does not show me the content of the function objects.

Has anyone an idea where I messed up my implementation?

I use "fenics:all/xenial 1:2017.1.0.1~ppa1~xenial1 uptodate" installed via apt-get on Ubuntu 16.04

Thanks in forward,
Justus

PS: my implementation is based on this pdf

Community: FEniCS Project