cs_runge_kutta_integrator.h

Go to the documentation of this file.

#ifndef RK_INTEGRATOR_H
#define RK_INTEGRATOR_H
 
/*============================================================================
 * Explicit Runge-Kutta integrator utilities.
 *============================================================================*/
 
/*
  This file is part of code_saturne, a general-purpose CFD tool.
 
  Copyright (C) 1998-2026 EDF S.A.
 
  This program is free software; you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation; either version 2 of the License, or (at your option) any later
  version.
 
  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  details.
 
  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
  Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
 
/*----------------------------------------------------------------------------*/
 
#include "base/cs_defs.h"
 
/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/
 
#include <assert.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
/*----------------------------------------------------------------------------
 * Local headers
 *----------------------------------------------------------------------------*/
 
#include "alge/cs_balance.h"
#include "alge/cs_blas.h"
 
#include "base/cs_array.h"
#include "base/cs_base_accel.h"
#include "base/cs_boundary_conditions_set_coeffs.h"
#include "base/cs_dispatch.h"
#include "mesh/cs_mesh.h"
 
/*----------------------------------------------------------------------------
 * Header for the current file
 *----------------------------------------------------------------------------*/
 
#include "base/cs_runge_kutta_integrator_param.h"
 
/*============================================================================
 * Type definitions
 *============================================================================*/
 
/*----------------------------------------------------------------------------*/
 
namespace cs {
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
class runge_kutta_integrator {
 
  /*===========================================================================
   * Public methods
   *==========================================================================*/
 
public:
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  runge_kutta_integrator()
  {}
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  runge_kutta_integrator
  (
    cs_runge_kutta_scheme_t  scheme, 
    const char              *name,   
    const cs_real_t         *dt,     
    int                      dim,    
    cs_lnum_t                n_elts  
  )
  {
    _scheme = scheme;
    CS_MALLOC(_name, strlen(name) + 1, char);
    strcpy(_name, name);
 
    _dt.update_data(const_cast<cs_real_t *>(dt), n_elts);
 
    _stride = dim;
    _n_elts = n_elts;
    _n_stages = 1;
    _i_stage = 0;
 
    _a.reshape(RK_HIGHEST_ORDER, RK_HIGHEST_ORDER);
    _c.reshape(RK_HIGHEST_ORDER);
 
    init_scheme_();
 
    _rhs_stages.reshape(_n_stages, _n_elts, _stride);
    _u_old.reshape(_n_elts, _stride);
    _u_new.reshape(_n_elts, _stride);
    _mass.reshape(_n_elts);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  ~runge_kutta_integrator()
  {
    CS_FREE(_name);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  init_state
  (
    cs_dispatch_context &ctx,  
    const cs_real_t     *rho,  
    const cs_real_t     *vol,  
    cs_real_t           *pvara 
  )
  {
    _i_stage = 0;
 
    auto mass = _mass.view();
    auto u_old = _u_old.view();
    const int stride = _stride;
 
    ctx.parallel_for(_n_elts, CS_LAMBDA (cs_lnum_t e_id) {
      mass[e_id] = rho[e_id] * vol[e_id];
      for (cs_lnum_t i = 0; i < stride; i++)
        u_old(e_id, i) = pvara[stride*e_id + i];
    });
    ctx.wait();
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  solve_stage
  (
    cs_dispatch_context &ctx,       
    cs_real_t           *pvar_stage 
  )
  {
    auto dt = _dt.view_1d();
    auto mass = _mass.view();
 
    const int i_stg = _i_stage;
    const int stride = _stride;
 
    auto u0 = _u_old.view();
    auto u_new = _u_new.view();
    auto rhs_stage = _rhs_stages.view();
 
    const auto a = _a.sub_view(i_stg);
 
    ctx.parallel_for(_n_elts, CS_LAMBDA (cs_lnum_t e_id) {
      for (cs_lnum_t i = 0; i < stride; i++)
        u_new(e_id, i) = u0(e_id, i);
 
      for (int j_stg = 0; j_stg <= i_stg; j_stg++) {
        for (cs_lnum_t i = 0; i < stride; i++)
          u_new(e_id, i) += rhs_stage(j_stg, e_id, i)* a[j_stg] * dt[e_id]
                          / mass[e_id];
      }
 
      for (cs_lnum_t i = 0; i < stride; i++)
        pvar_stage[stride*e_id + i] = u_new(e_id, i);
    });
 
    ctx.wait();
 
    _i_stage += 1;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  stage_projection_rhs
  (
    cs_dispatch_context &ctx,    
    cs_real_t           *grad_dp 
  )
  {
    const int i_stg = _i_stage - 1;
    const auto a = _a.sub_view(i_stg);
 
    auto rhs = _rhs_stages.sub_view(i_stg);
    auto mass = _mass.view();
    const int stride = _stride;
 
    ctx.parallel_for(_n_elts, CS_LAMBDA (cs_lnum_t e_id) {
      for (cs_lnum_t i = 0; i < stride; i++) {
        rhs(e_id, i) -= grad_dp[stride*e_id + i] * mass[e_id];
        grad_dp[stride*e_id + i] *= a[i_stg];
      }
    });
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  stage_set_initial_rhs
  (
    cs_dispatch_context &ctx,     
    cs_real_t           *rhs_pvar 
  )
  {
    const int i_stg = _i_stage;
    auto rhs = _rhs_stages.sub_view(i_stg);
    const int stride = _stride;
 
    ctx.parallel_for(_n_elts, CS_LAMBDA (cs_lnum_t e_id) {
      for (cs_lnum_t i = 0; i < stride; i++)
        rhs(e_id, i) = rhs_pvar[stride*e_id + i];
    });
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline bool
  is_active() const
  {
    return (_scheme > CS_RK_NONE);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline bool
  is_staging() const
  {
    return (_i_stage < _n_stages);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline cs_runge_kutta_scheme_t
  scheme() const
  {
    return _scheme;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline cs_lnum_t
  n_elts() const
  {
    return _n_elts;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline int
  i_stage() const
  {
    return _i_stage;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline int
  stride() const
  {
    return _stride;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  cs_real_t *
  get_rhs_stage_sub_array
  (
    int stage_id 
  )
  {
    return _rhs_stages.sub_array(stage_id);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
  CS_F_HOST
  cs_span<cs_real_t>
  get_stage_coeff_a
  (
    int i 
  )
  {
    return _a.sub_view(i);
  }
 
private:
 
  /*===========================================================================
   * Private methods
   *==========================================================================*/
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  init_scheme_()
  {
    _a.zero();
    _c.zero();
 
    switch (_scheme) {
    case CS_RK_NONE:
      return;
 
    case CS_RK1:
      _n_stages = 1;
      _a(0, 0) = 1.;
      _c[0]    = 1.;
 
      break;
 
    case CS_RK2:
      [[fallthrough]];
    case CS_RK2_HEUN:
      _n_stages = 2;
      _a(0,0) = 1.;
      _a(1,0) = 0.5, _a(1,1) = 0.5;
 
      _c[0] = 1.;
      _c[1] = 1.;
 
      break;
 
    case CS_RK2_MP:
      _n_stages = 2;
      _a(0,0) = 0.5;
      _a(1,0) = 0. , _a(1,1) = 1.;
 
      _c[0] = 0.5;
      _c[1] = 1.;
 
      break;
 
    case CS_RK3:
      [[fallthrough]];
    case CS_RK3_WRAY:
      _n_stages = 3;
 
      _a(0,0) = 8. / 15.;
      _a(1,0) = 1./4., _a(1,1) = 5./12.;
      _a(2,0) = 1./4., _a(2,1) = 0., _a(2,2) = 3./4.;
 
      _c[0] = 8./15.;
      _c[1] = 2./3.;
      _c[2] = 1.;
 
      break;
 
    case CS_RK3_SSP:
      _n_stages = 3;
      _a(0,0) = 1.;
      _a(1,0) = 1./4., _a(1,1) = 1./4.;
      _a(2,0) = 1./6., _a(2,1) = 1./6., _a(2,2) = 2./3.;
 
      _c[0] = 1.;
      _c[1] = 1./2.;
      _c[2] = 1.;
 
      break;
 
    case CS_RK4:
      _n_stages = 4;
      _a(0,0) = 1.;
      _a(1,0) = 3./8., _a(1,1) = 1./8.;
      _a(2,0) = -1./8., _a(2,1) = -3./8., _a(2,2) = 3./2.;
      _a(3,0) = 1./6., _a(3,1) = -1./18., _a(3,2) = 2./3., _a(3,3) = 2./9.;
 
      _c[0] = 1.;
      _c[1] = 1./2.;
      _c[2] = 1.;
      _c[3] = 1.;
 
      break;
 
    default:
      bft_error(__FILE__, __LINE__, 0,
                "%s: Type of Runge-Kutta not available.\n"
                "%s: Stop building Runge-Kutta integrator.\n",
                __func__, __func__);
      break;
    }
  }
 
  /*===========================================================================
   * Private members
   *==========================================================================*/
 
  cs_runge_kutta_scheme_t  _scheme{CS_RK_NONE}; 
  char                    *_name{nullptr};      
  int                      _n_stages{0};        
  int                      _i_stage{0};         
  cs_span<cs_real_t>       _dt;                 
  int                      _stride{0};          
  cs_lnum_t                _n_elts{0};          
  /* variable storage */
  cs_array_2d<cs_real_t>   _u_old;              
  cs_array_2d<cs_real_t>   _u_new;              
  cs_array_3d<cs_real_t>   _rhs_stages;         
  cs_array<cs_real_t>      _mass;               
  /* Coeffients adapted from the Butcher tableau */
  cs_array_2d<cs_real_t>   _a;            
  cs_array<cs_real_t>      _c;            
};
}
 
using cs_runge_kutta_integrator_t = cs::runge_kutta_integrator;
 
/*============================================================================
 * Global variables
 *============================================================================*/
 
/*=============================================================================
 * Public function prototypes
 *============================================================================*/
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
int
cs_runge_kutta_integrator_create(cs_runge_kutta_scheme_t   scheme,
                                 const char               *name,
                                 const cs_real_t          *dt,
                                 int                       dim,
                                 cs_lnum_t                 n_elts);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_runge_kutta_integrators_initialize();
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_runge_kutta_staging_potential(cs_dispatch_context          &ctx,
                                 cs_runge_kutta_integrator_t  *rk,
                                 cs_real_t                    *phi_stage);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_runge_kutta_stage_complete_scalar_rhs
  (cs_dispatch_context         &ctx,
   cs_runge_kutta_integrator_t *rk,
   int                          idtvar,
   int                          f_id,
   int                          imucpp,
   cs_equation_param_t         *eqp,
   cs_field_bc_coeffs_t        *bc_coeffs,
   const cs_real_t              i_massflux[],
   const cs_real_t              b_massflux[],
   const cs_real_t              i_visc[],
   const cs_real_t              b_visc[],
   cs_real_t                    viscel[][6],
   const cs_real_t              weighf[][2],
   const cs_real_t              weighb[],
   int                          icvflb,
   const int                    icvfli[],
   cs_real_t                    pvar[],
   const cs_real_t              xcpp[]);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template<int stride>
void
cs_runge_kutta_stage_complete_rhs(cs_dispatch_context         &ctx,
                                  cs_runge_kutta_integrator_t *rk,
                                  int                          idtvar,
                                  int                          f_id,
                                  int                          ivisep,
                                  cs_equation_param_t         *eqp,
                                  cs_field_bc_coeffs_t        *bc_coeffs,
                                  const cs_real_t              i_massflux[],
                                  const cs_real_t              b_massflux[],
                                  const cs_real_t              i_visc[],
                                  const cs_real_t              b_visc[],
                                  const cs_real_t              i_secvis[],
                                  const cs_real_t              b_secvis[],
                                  cs_real_t                    viscel[][6],
                                  const cs_real_t              weighf[][2],
                                  const cs_real_t              weighb[],
                                  int                          icvflb,
                                  const int                    icvfli[],
                                  cs_real_t                    pvar[][stride])
{
  // Sanity check
  assert(rk != nullptr);
  assert(stride == rk->stride());
 
  // get the current stage index
  const int i_stg = rk->i_stage();
 
  cs_real_t *rhs = rk->get_rhs_stage_sub_array(i_stg);
 
  /* Allocate non reconstructed face value only if presence of limiter */
 
  /* We compute the total explicit balance. */
 
  int  inc = 1;
 
    /* The added convective scalar mass flux is:
     *      (thetex*Y_\face-imasac*Y_\celli)*mf.
     * When building the explicit part of the rhs, one
     * has to impose 0 on mass accumulation. */
  int  imasac = 0;
 
  eqp->theta = 1;
 
  if (stride == 3)
    cs_balance_vector(idtvar,
                      f_id,
                      imasac,
                      inc,       /* inc */
                      ivisep,
                      eqp,
                      nullptr, /* pvar == pvara */
                      (const cs_real_3_t *)pvar,
                      bc_coeffs,
                      i_massflux,
                      b_massflux,
                      i_visc,
                      b_visc,
                      i_secvis,
                      b_secvis,
                      viscel,
                      weighf,
                      weighb,
                      icvflb,
                      icvfli,
                      nullptr,
                      nullptr,
                      (cs_real_3_t *)rhs);
 
  else if (stride == 6)
    cs_balance_tensor(idtvar,
                      f_id,
                      imasac,
                      inc,       /* inc */
                      eqp,
                      nullptr, /* pvar == pvara */
                      (const cs_real_6_t *)pvar,
                      bc_coeffs,
                      i_massflux,
                      b_massflux,
                      i_visc,
                      b_visc,
                      viscel,
                      weighf,
                      weighb,
                      icvflb,
                      icvfli,
                      (cs_real_6_t *)rhs);
 
  eqp->theta = 0;
 
  ctx.wait();
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
inline bool
cs_runge_kutta_is_active(cs_runge_kutta_integrator_t  *rk)
{
  if (rk == nullptr)
    return false;
  else
    return rk->is_active();
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
inline bool
cs_runge_kutta_is_staging(cs_runge_kutta_integrator_t  *rk)
{
  if (rk == nullptr)
    return false;
 
  return rk->is_staging();
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
cs_runge_kutta_integrator_t *
cs_runge_kutta_integrator_by_id(int  rk_id);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
void
cs_runge_kutta_integrators_destroy();
 
/*----------------------------------------------------------------------------*/
 
#endif /* __RK_INTEGRATOR_H */