nvs_bdf_q.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
#include <kernel/base_header.hpp>
#include <control/time/base.hpp>
 
#include "vector"
 
namespace FEAT
{
  namespace Control
  {
    namespace Time
    {
       template<typename DT_>
       class NvsBdfQ
       {
         public:
           typedef DT_ DataType;
 
         private:
           Index _num_steps;
           Index _num_steps_startup;
           Index _p_extrapolation_steps;
           Index _p_extrapolation_steps_startup;
           const DataType _reynolds;
           DataType _delta_t;
           bool _use_deformation;
           bool _use_rotational_form;
           bool _startup;
 
         public:
           std::vector<DataType> coeff_lhs_v;
           std::vector<DataType> coeff_rhs_v_1;
           std::vector<DataType> coeff_rhs_v_phi;
           std::vector<DataType> coeff_extrapolation_p;
 
           DataType coeff_rhs_f;
           DataType coeff_rhs_proj_D_v;
           DataType coeff_rhs_phi;
           DataType coeff_rhs_v_2;
           DataType coeff_rhs_v_p;
 
           TermHandling ale_handling;
           TermHandling conv_handling;
           TermHandling visc_handling;
 
           NvsBdfQ(PropertyMap* config_section, const DataType reynolds, const bool use_deformation,
           const bool startup):
             _num_steps(~Index(0)),
             _num_steps_startup(1),
             _p_extrapolation_steps(~Index(0)),
             _p_extrapolation_steps_startup(~Index(0)),
             _reynolds(reynolds),
             _delta_t(0),
             _use_deformation(use_deformation),
             _use_rotational_form(false),
             _startup(startup),
             coeff_lhs_v(3),
             coeff_rhs_v_1(3),
             coeff_rhs_v_phi(0),
             coeff_extrapolation_p(0),
             coeff_rhs_f(0),
             coeff_rhs_proj_D_v(0),
             coeff_rhs_phi(0),
             coeff_rhs_v_2(0),
             coeff_rhs_v_p(0),
             ale_handling(TermHandling::off),
             conv_handling(TermHandling::off),
             visc_handling(TermHandling::off)
             {
               XASSERT(config_section != nullptr);
               // Get time discretization order
               auto num_steps_p = config_section->query("num_steps");
               XASSERTM(num_steps_p.second, "TimeDiscretization section is missing the mandatory num_steps entry!");
               _num_steps = Index(std::stoul(num_steps_p.first));
               XASSERTM(_num_steps == Index(1) || _num_steps == Index(2),"Only BDF1 and BDF2 are implemented!");
 
               // Get the order for extrapolation of the pressure in time
               auto p_extrapolation_p = config_section->query("p_extrapolation_steps");
               XASSERTM(p_extrapolation_p.second,
               "TimeDiscretization section is missing the mandatory p_extrapolation_steps entry!");
               _p_extrapolation_steps = Index(std::stoul(p_extrapolation_p.first));
               XASSERTM(_p_extrapolation_steps < Index(3),"p extrapolation is only implemented for order 0, 1, 2!");
               _p_extrapolation_steps_startup = Math::min(_p_extrapolation_steps, Index(1));
 
               // Use the rotational form?
               auto use_rotational_form_p = config_section->query("use_rotational_form");
               XASSERTM(use_rotational_form_p.second,
               "TimeDiscretization section is missing the mandatory use_rotational_form entry!");
               XASSERTM(std::stoi(use_rotational_form_p.first) == 0 || std::stoi(use_rotational_form_p.first) == 1,
               "use_rotational has to be set to 0 or 1");
               _use_rotational_form = std::stoi(use_rotational_form_p.first) == 1;
 
               // Get timestep size
               auto delta_t_p = config_section->query("delta_t");
               XASSERTM(delta_t_p.second, "TimeDiscretization section is missing the mandatory delta_t entry!");
               _delta_t = DataType(std::stod(delta_t_p.first));
               XASSERT(_delta_t > DataType(0));
 
               // Use ale (meaning solve Navier-Stokes instead of plain Stokes)?
               auto ale_p = config_section->query("ALE");
               XASSERTM(ale_p.second,
               "TimeDiscretization section is missing the mandatory ALE entry!");
               ale_p.first.parse(ale_handling);
 
               // Use convection (meaning solve Navier-Stokes instead of plain Stokes)?
               auto convection_p = config_section->query("convection");
               XASSERTM(convection_p.second,
               "TimeDiscretization section is missing the mandatory convection entry!");
               convection_p.first.parse(conv_handling);
 
               // Use viscous (meaning solve Navier-Stokes instead of plain Stokes)?
               auto viscous_p = config_section->query("viscous");
               XASSERTM(viscous_p.second,
               "TimeDiscretization section is missing the mandatory viscous entry!");
               viscous_p.first.parse(visc_handling);
 
               coeff_rhs_v_phi.clear();
               coeff_rhs_v_phi = std::vector<DataType>(_num_steps, DataType(0));
 
               coeff_extrapolation_p.clear();
               coeff_extrapolation_p = std::vector<DataType>(_p_extrapolation_steps, DataType(0));
 
               if(startup)
               {
                 set_coeffs(_num_steps_startup, _p_extrapolation_steps_startup);
               }
               else
               {
                 set_coeffs(_num_steps, _p_extrapolation_steps);
               }
             }
 
           virtual ~NvsBdfQ()
           {
           }
 
           DataType delta_t()
           {
             return _delta_t;
           }
 
           bool is_rotational()
           {
             return _use_rotational_form;
           }
 
           Index get_num_steps()
           {
             return _startup ? _num_steps_startup : _num_steps;
           }
 
           Index get_max_num_steps()
           {
             return Math::max(_num_steps, _num_steps_startup);
           }
 
           Index get_max_p_extrapolation_steps()
           {
             return Math::max(_p_extrapolation_steps, _p_extrapolation_steps_startup);
           }
 
           Index get_p_extrapolation_steps()
           {
             return _startup ? _p_extrapolation_steps_startup : _p_extrapolation_steps;
           }
 
           void finish_startup()
           {
             _startup = false;
 
             set_coeffs(_num_steps, _p_extrapolation_steps);
           }
 
           void set_coeffs(const Index num_steps, const Index p_extrapolation_steps)
           {
 
             if(num_steps == Index(1))
             {
               // Mass matrix
               coeff_lhs_v.at(0) = DataType(1);
               // Convection matrix
               coeff_lhs_v.at(1) = (conv_handling == TermHandling::impl || ale_handling == TermHandling::impl)
                 ? _delta_t : DataType(0);
               // Stiffness matrix
               coeff_lhs_v.at(2) = (visc_handling == TermHandling::impl) ? _delta_t/_reynolds : DataType(0);
 
               // Mass matrix
               coeff_rhs_v_1.at(0) = DataType(1);
               // Convection matrix
               coeff_rhs_v_1.at(1) = (conv_handling == TermHandling::expl || ale_handling == TermHandling::expl)
                 ? -_delta_t : DataType(0);
               // Stiffness matrix
               coeff_rhs_v_1.at(2) = (visc_handling == TermHandling::expl) ? -_delta_t/_reynolds : DataType(0);
 
               // Mass matrix coefficient for previous time step vector
               coeff_rhs_v_2 = DataType(0);
 
               coeff_rhs_f = _delta_t;
 
               // coefficient for -(p*[k], div phi) on the rhs for the velocity eq.
               coeff_rhs_v_p = -_delta_t;
               // coefficient for -(phi[k], div phi) on the rhs for the velocity eq.
               coeff_rhs_v_phi.at(0) = -_delta_t;
 
               coeff_rhs_phi = DataType(1)/_delta_t;
 
             }
             else if(num_steps == Index(2))
             {
               // Mass matrix
               coeff_lhs_v.at(0) = DataType(1);
               // Convection matrix
               coeff_lhs_v.at(1) = (conv_handling == TermHandling::impl) ? DataType(2)/DataType(3)*_delta_t : DataType(0);
               // Stiffness matrix
               coeff_lhs_v.at(2) = (visc_handling == TermHandling::impl || ale_handling == TermHandling::impl)
                 ? DataType(2)/DataType(3)*_delta_t/_reynolds : DataType(0);
 
               // Mass matrix
               coeff_rhs_v_1.at(0) = DataType(4)/DataType(3);
               // Convection matrix
               coeff_rhs_v_1.at(1) = (conv_handling == TermHandling::expl || ale_handling == TermHandling::expl)
                 ? -DataType(2)/DataType(3)*_delta_t : DataType(0);
               // Stiffness matrix
               coeff_rhs_v_1.at(2) =
                 (visc_handling == TermHandling::expl) ? -DataType(2)/DataType(3)*_delta_t/_reynolds : DataType(0);
 
               coeff_rhs_f = DataType(2)/DataType(3)*_delta_t;
 
               // Mass matrix coefficient for previous time step vector
               coeff_rhs_v_2 = -DataType(1)/DataType(3);
 
               // coefficient for -(p*[k], div phi) on the rhs for the velocity eq.
               coeff_rhs_v_p = -DataType(2)/DataType(3)*_delta_t;
 
               // coefficient for -(phi[k], div phi) on the rhs for the velocity eq.
               coeff_rhs_v_phi.at(0) =  -DataType(8)/DataType(9)*_delta_t;
               // coefficient for -(phi[k-1], div phi) on the rhs for the velocity eq.
               coeff_rhs_v_phi.at(1) = DataType(2)/DataType(9)*_delta_t;
 
               coeff_rhs_phi = DataType(3)/(DataType(2)*_delta_t);
 
             }
             else
             {
               XABORTM("Invalid number of BDF steps: "+stringify(num_steps));
             }
 
             // Coefficient for the div v term for the projection step
             coeff_rhs_proj_D_v = _use_rotational_form ? DataType(1)/_reynolds : DataType(0);
 
             if(_use_deformation)
             {
               coeff_rhs_proj_D_v *= DataType(2);
             }
 
             // Coefficients for the extrapolation of the pressure
             if(p_extrapolation_steps == Index(1))
             {
               coeff_extrapolation_p.at(0) = DataType(1);
             }
             else if(p_extrapolation_steps == Index(2))
             {
               coeff_extrapolation_p.at(0) = DataType(2);
               coeff_extrapolation_p.at(1) = -DataType(1);
             }
             else
             {
               XABORTM("Invalid number of pressure extrapolation steps: "+stringify(p_extrapolation_steps));
             }
 
           }
       };
    } // namespace Time
  } // namespace Control
} // name FEAT