ODEint.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#include "ODEint.h"
 
namespace limpet {
 
void ODEint_RK(GlobalData_t *y, GlobalData_t *ydot, void f(GlobalData_t*,GlobalData_t*,void*), void *f_data, int len, GlobalData_t dt)
{
  int     i;
  GlobalData_t  k1[ODEint_MaxVecLen];
  GlobalData_t *ynpk1 = k1;
  GlobalData_t *k2    = k1;
 
  for (i=0;i<len;i++)  {
    k1[i]    = dt*ydot[i];
    ynpk1[i] = y[i]+0.5*k1[i];
  }
 
  f(ynpk1,ydot,f_data);
 
  for (i=0;i<len;i++) {
    k2[i] = dt*ydot[i];
    y[i] += k2[i];
  }
}
 
void ODEint_FE(GlobalData_t *y, GlobalData_t *ydot, void *f, void *f_data, int len, GlobalData_t dt)
{
  for (int i=0;i<len;i++)
    y[i] += ydot[i]*dt;
}
 
 
void d_ODEint_RK(d_OdeVec *y, d_OdeVec *ydot, void f(GlobalData_t*,d_OdeVec*,void*),
                 void *f_data, GlobalData_t dt)
{
  int     i;
  GlobalData_t  k1[ODEint_MaxVecLen];
  GlobalData_t *ynpk1 = k1;
  GlobalData_t *k2    = k1;
 
  for (i=0;i<y->len;i++)  {
    k1[i]    = dt*ydot->vec[i];
    ynpk1[i] = y->vec[i]+0.5*k1[i];
  }
 
  f(ynpk1,ydot,f_data);
 
  for (i=0;i<y->len;i++) {
    k2[i]      = dt*ydot->vec[i];
    y->vec[i] += k2[i];
  }
}
 
void d_ODEint_FE(d_OdeVec*y, d_OdeVec *ydot, void *f, void *f_data, GlobalData_t dt)
{
  for (int i=0;i<y->len;i++)
    y->vec[i] += ydot->vec[i]*dt;
}
 
}  // namespace limpet