Forked from
Nektar / Nektar
184 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
DIRKTimeIntegrationSchemes.h 17.64 KiB
///////////////////////////////////////////////////////////////////////////////
//
// File: DIRKTimeIntegrationSchemes.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2018 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Combined header file for all basic DIRK based time integration
// schemes.
//
///////////////////////////////////////////////////////////////////////////////
// Note : If adding a new integrator be sure to register the
// integrator with the Time Integration Scheme Facatory in
// SchemeInitializor.cpp.
#ifndef NEKTAR_LIB_UTILITIES_TIME_INTEGRATION_DIRK_TIME_INTEGRATION_SCHEME
#define NEKTAR_LIB_UTILITIES_TIME_INTEGRATION_DIRK_TIME_INTEGRATION_SCHEME
#define LUE LIB_UTILITIES_EXPORT
#include <LibUtilities/TimeIntegration/TimeIntegrationSchemeGLM.h>
namespace Nektar::LibUtilities
{
///////////////////////////////////////////////////////////////////////////////
// DIRK Order N
class DIRKTimeIntegrationScheme : public TimeIntegrationSchemeGLM
{
public:
DIRKTimeIntegrationScheme(std::string variant, size_t order,
std::vector<NekDouble> freeParams)
: TimeIntegrationSchemeGLM(variant, order, freeParams)
{
// Currently up to 4th order is implemented.
ASSERTL0(1 <= order && order <= 4,
"Diagonally Implicit Runge Kutta integration scheme bad order "
"(1-4): " +
std::to_string(order));
m_integration_phases = TimeIntegrationAlgorithmGLMVector(1);
m_integration_phases[0] = TimeIntegrationAlgorithmGLMSharedPtr(
new TimeIntegrationAlgorithmGLM(this));
ASSERTL1((variant == "" || variant == "ES5" || variant == "ES6"), "DIRK Time integration scheme bad variant: " + variant +
". "
"Must blank, 'ES5', or 'ES6'");
if ("" == variant)
{
DIRKTimeIntegrationScheme::SetupSchemeData(m_integration_phases[0],
order);
}
else
{
DIRKTimeIntegrationScheme::SetupSchemeDataESDIRK(
m_integration_phases[0], variant, order, freeParams);
}
}
~DIRKTimeIntegrationScheme() override
{
}
static TimeIntegrationSchemeSharedPtr create(
std::string variant, size_t order, std::vector<NekDouble> freeParams)
{
TimeIntegrationSchemeSharedPtr p =
MemoryManager<DIRKTimeIntegrationScheme>::AllocateSharedPtr(
variant, order, freeParams);
return p;
}
static std::string className;
LUE static void SetupSchemeData(TimeIntegrationAlgorithmGLMSharedPtr &phase,
size_t order)
{
phase->m_schemeType = eDiagonallyImplicit;
phase->m_order = order;
phase->m_name =
std::string("DIRKOrder" + std::to_string(phase->m_order));
phase->m_numsteps = 1;
phase->m_numstages = phase->m_order;
phase->m_A = Array<OneD, Array<TwoD, NekDouble>>(1);
phase->m_B = Array<OneD, Array<TwoD, NekDouble>>(1);
phase->m_A[0] =
Array<TwoD, NekDouble>(phase->m_numstages, phase->m_numstages, 0.0);
phase->m_B[0] =
Array<TwoD, NekDouble>(phase->m_numsteps, phase->m_numstages, 0.0);
phase->m_U =
Array<TwoD, NekDouble>(phase->m_numstages, phase->m_numsteps, 1.0);
phase->m_V =
Array<TwoD, NekDouble>(phase->m_numsteps, phase->m_numsteps, 1.0);
switch (phase->m_order)
{
case 1:
{
// One-stage, 1st order, L-stable Diagonally Implicit
// Runge Kutta (backward Euler) method:
phase->m_A[0][0][0] = 1.0;
phase->m_B[0][0][0] = 1.0;
}
break;
case 2:
{
// Two-stage, 2nd order Diagonally Implicit Runge
// Kutta method. It is A-stable if and only if x ≥ 1/4 // and is L-stable if and only if x equals one of the
// roots of the polynomial x^2 - 2x + 1/2, i.e. if
// x = 1 ± sqrt(2)/2.
const NekDouble lambda = (2.0 - sqrt(2.0)) / 2.0;
phase->m_A[0][0][0] = lambda;
phase->m_A[0][1][0] = 1.0 - lambda;
phase->m_A[0][1][1] = lambda;
phase->m_B[0][0][0] = 1.0 - lambda;
phase->m_B[0][0][1] = lambda;
}
break;
case 3:
{
// Three-stage, 3rd order, L-stable Diagonally Implicit
// Runge Kutta method:
constexpr NekDouble lambda = 0.4358665215;
phase->m_A[0][0][0] = lambda;
phase->m_A[0][1][0] = 0.5 * (1.0 - lambda);
phase->m_A[0][2][0] =
0.25 * (-6.0 * lambda * lambda + 16.0 * lambda - 1.0);
phase->m_A[0][1][1] = lambda;
phase->m_A[0][2][1] =
0.25 * (6.0 * lambda * lambda - 20.0 * lambda + 5.0);
phase->m_A[0][2][2] = lambda;
phase->m_B[0][0][0] =
0.25 * (-6.0 * lambda * lambda + 16.0 * lambda - 1.0);
phase->m_B[0][0][1] =
0.25 * (6.0 * lambda * lambda - 20.0 * lambda + 5.0);
phase->m_B[0][0][2] = lambda;
}
break;
}
phase->m_numMultiStepValues = 1;
phase->m_numMultiStepImplicitDerivs = 0;
phase->m_numMultiStepExplicitDerivs = 0;
phase->m_timeLevelOffset = Array<OneD, size_t>(phase->m_numsteps);
phase->m_timeLevelOffset[0] = 0;
phase->CheckAndVerify();
}
LUE static void SetupSchemeDataESDIRK(
TimeIntegrationAlgorithmGLMSharedPtr &phase, std::string variant,
size_t order, std::vector<NekDouble> freeParams)
{
phase->m_schemeType = eDiagonallyImplicit;
phase->m_order = order;
phase->m_name =
std::string("DIRKOrder" + std::to_string(phase->m_order) + variant);
phase->m_numsteps = 1;
char const &stage = variant.back();
phase->m_numstages = std::atoi(&stage);
phase->m_A = Array<OneD, Array<TwoD, NekDouble>>(1);
phase->m_B = Array<OneD, Array<TwoD, NekDouble>>(1);
phase->m_A[0] =
Array<TwoD, NekDouble>(phase->m_numstages, phase->m_numstages, 0.0);
phase->m_B[0] =
Array<TwoD, NekDouble>(phase->m_numsteps, phase->m_numstages, 0.0);
phase->m_U = Array<TwoD, NekDouble>(phase->m_numstages, phase->m_numsteps, 1.0);
phase->m_V =
Array<TwoD, NekDouble>(phase->m_numsteps, phase->m_numsteps, 1.0);
const NekDouble ConstSqrt2 = sqrt(2.0);
switch (phase->m_order)
{
case 3:
{
// See: Kennedy, Christopher A., and Mark H. Carpenter.
// Diagonally implicit Runge-Kutta methods for ordinary
// differential equations. A review. No. NF1676L-19716. 2016.
ASSERTL0(5 == phase->m_numstages,
std::string("DIRKOrder3_ES" +
std::to_string(phase->m_numstages) +
" not defined"));
NekDouble lambda =
(freeParams.size()) ? freeParams[0] : 9.0 / 40.0;
phase->m_A[0][0][0] = 0.0;
phase->m_A[0][1][0] = lambda;
phase->m_A[0][2][0] = 9.0 * (1.0 + ConstSqrt2) / 80.0;
phase->m_A[0][3][0] =
(22.0 + 15.0 * ConstSqrt2) / (80.0 * (1.0 + ConstSqrt2));
phase->m_A[0][4][0] = (2398.0 + 1205.0 * ConstSqrt2) /
(2835.0 * (4.0 + 3.0 * ConstSqrt2));
phase->m_A[0][1][1] = phase->m_A[0][1][0];
phase->m_A[0][2][1] = phase->m_A[0][2][0];
phase->m_A[0][3][1] = phase->m_A[0][3][0];
phase->m_A[0][4][1] = phase->m_A[0][4][0];
phase->m_A[0][2][2] = lambda;
phase->m_A[0][3][2] = -7.0 / (40.0 * (1.0 + ConstSqrt2));
phase->m_A[0][4][2] = -2374.0 * (1.0 + 2.0 * ConstSqrt2) /
(2835.0 * (5.0 + 3.0 * ConstSqrt2));
phase->m_A[0][3][3] = lambda;
phase->m_A[0][4][3] = 5827.0 / 7560.0;
phase->m_A[0][4][4] = lambda;
phase->m_B[0][0][0] = phase->m_A[0][4][0];
phase->m_B[0][0][1] = phase->m_A[0][4][1];
phase->m_B[0][0][2] = phase->m_A[0][4][2];
phase->m_B[0][0][3] = phase->m_A[0][4][3];
phase->m_B[0][0][4] = phase->m_A[0][4][4];
}
break;
case 4:
{
// See: Kennedy, Christopher A., and Mark H. Carpenter.
// Diagonally implicit Runge-Kutta methods for ordinary
// differential equations. A review. No. NF1676L-19716. 2016.
ASSERTL0(6 == phase->m_numstages,
std::string("DIRKOrder4_ES" +
std::to_string(phase->m_numstages) +
" not defined"));
NekDouble lambda =
(freeParams.size()) ? freeParams[0] : 1.0 / 4.0;
phase->m_A[0][0][0] = 0.0;
phase->m_A[0][1][0] = lambda;
phase->m_A[0][2][0] = (1.0 - ConstSqrt2) / 8.0;
phase->m_A[0][3][0] = (5.0 - 7.0 * ConstSqrt2) / 64.0;
phase->m_A[0][4][0] =
(-13796.0 - 54539.0 * ConstSqrt2) / 125000.0; phase->m_A[0][5][0] = (1181.0 - 987.0 * ConstSqrt2) / 13782.0;
phase->m_A[0][1][1] = phase->m_A[0][1][0];
phase->m_A[0][2][1] = phase->m_A[0][2][0];
phase->m_A[0][3][1] = phase->m_A[0][3][0];
phase->m_A[0][4][1] = phase->m_A[0][4][0];
phase->m_A[0][5][1] = phase->m_A[0][5][0];
phase->m_A[0][2][2] = lambda;
phase->m_A[0][3][2] = 7.0 * (1.0 + ConstSqrt2) / 32.0;
phase->m_A[0][4][2] =
(506605.0 + 132109.0 * ConstSqrt2) / 437500.0;
phase->m_A[0][5][2] =
47.0 * (-267.0 + 1783.0 * ConstSqrt2) / 273343.0;
phase->m_A[0][3][3] = lambda;
phase->m_A[0][4][3] =
166.0 * (-97.0 + 376.0 * ConstSqrt2) / 109375.0;
phase->m_A[0][5][3] =
-16.0 * (-22922.0 + 3525.0 * ConstSqrt2) / 571953.0;
phase->m_A[0][4][4] = lambda;
phase->m_A[0][5][4] =
-15625.0 * (97.0 + 376.0 * ConstSqrt2) / 90749876.0;
phase->m_A[0][5][5] = lambda;
phase->m_B[0][0][0] = phase->m_A[0][5][0];
phase->m_B[0][0][1] = phase->m_A[0][5][1];
phase->m_B[0][0][2] = phase->m_A[0][5][2];
phase->m_B[0][0][3] = phase->m_A[0][5][3];
phase->m_B[0][0][4] = phase->m_A[0][5][4];
phase->m_B[0][0][5] = phase->m_A[0][5][5];
}
break;
default:
{
ASSERTL0(false, std::string("ESDIRK of order" +
std::to_string(phase->m_order) +
" not defined"));
break;
}
}
phase->m_numMultiStepValues = 1;
phase->m_numMultiStepImplicitDerivs = 0;
phase->m_numMultiStepExplicitDerivs = 0;
phase->m_timeLevelOffset = Array<OneD, size_t>(phase->m_numsteps);
phase->m_timeLevelOffset[0] = 0;
phase->CheckAndVerify();
}
protected:
LUE std::string v_GetName() const override
{
return std::string("DIRK");
}
LUE NekDouble v_GetTimeStability() const override
{
return 1.0;
}
}; // end class DIRKTimeIntegrationScheme
////////////////////////////////////////////////////////////////////////////////
// Backwards compatibility
class DIRKOrder1TimeIntegrationScheme : public DIRKTimeIntegrationScheme
{public:
DIRKOrder1TimeIntegrationScheme(std::string variant, size_t order,
std::vector<NekDouble> freeParams)
: DIRKTimeIntegrationScheme("", 1, freeParams)
{
boost::ignore_unused(variant, order);
}
static TimeIntegrationSchemeSharedPtr create(
[[maybe_unused]] std::string variant, [[maybe_unused]] size_t order,
std::vector<NekDouble> freeParams)
{
TimeIntegrationSchemeSharedPtr p =
MemoryManager<DIRKTimeIntegrationScheme>::AllocateSharedPtr(
"", 1, freeParams);
return p;
}
static std::string className;
protected:
static std::string TimeIntegrationMethodLookupId;
}; // end class DIRKOrder1TimeIntegrationScheme
class DIRKOrder2TimeIntegrationScheme : public DIRKTimeIntegrationScheme
{
public:
DIRKOrder2TimeIntegrationScheme(std::string variant, size_t order,
std::vector<NekDouble> freeParams)
: DIRKTimeIntegrationScheme("", 2, freeParams)
{
boost::ignore_unused(variant, order);
}
static TimeIntegrationSchemeSharedPtr create(
[[maybe_unused]] std::string variant, [[maybe_unused]] size_t order,
std::vector<NekDouble> freeParams)
{
TimeIntegrationSchemeSharedPtr p =
MemoryManager<DIRKTimeIntegrationScheme>::AllocateSharedPtr(
"", 2, freeParams);
return p;
}
static std::string className;
protected:
static std::string TimeIntegrationMethodLookupId;
}; // end class DIRKOrder2TimeIntegrationScheme
class DIRKOrder3TimeIntegrationScheme : public DIRKTimeIntegrationScheme
{
public:
DIRKOrder3TimeIntegrationScheme(std::string variant, size_t order,
std::vector<NekDouble> freeParams)
: DIRKTimeIntegrationScheme("", 3, freeParams)
{
boost::ignore_unused(variant, order);
}
static TimeIntegrationSchemeSharedPtr create(
[[maybe_unused]] std::string variant, [[maybe_unused]] size_t order,
std::vector<NekDouble> freeParams)
{
TimeIntegrationSchemeSharedPtr p =
MemoryManager<DIRKTimeIntegrationScheme>::AllocateSharedPtr( "", 3, freeParams);
return p;
}
static std::string className;
protected:
static std::string TimeIntegrationMethodLookupId;
}; // end class DIRKOrder3TimeIntegrationScheme
class DIRKOrder3_ES5TimeIntegrationScheme : public DIRKTimeIntegrationScheme
{
public:
DIRKOrder3_ES5TimeIntegrationScheme(std::string variant, size_t order,
std::vector<NekDouble> freeParams)
: DIRKTimeIntegrationScheme("ES5", 3, freeParams)
{
boost::ignore_unused(variant, order);
}
static TimeIntegrationSchemeSharedPtr create(
[[maybe_unused]] std::string variant, [[maybe_unused]] size_t order,
std::vector<NekDouble> freeParams)
{
TimeIntegrationSchemeSharedPtr p =
MemoryManager<DIRKTimeIntegrationScheme>::AllocateSharedPtr(
"ES5", 3, freeParams);
return p;
}
static std::string className;
protected:
static std::string TimeIntegrationMethodLookupId;
}; // end class DIRKOrder3_ES5TimeIntegrationScheme
class DIRKOrder4_ES6TimeIntegrationScheme : public DIRKTimeIntegrationScheme
{
public:
DIRKOrder4_ES6TimeIntegrationScheme(std::string variant, size_t order,
std::vector<NekDouble> freeParams)
: DIRKTimeIntegrationScheme("ES6", 4, freeParams)
{
boost::ignore_unused(variant, order);
}
static TimeIntegrationSchemeSharedPtr create(
[[maybe_unused]] std::string variant, [[maybe_unused]] size_t order,
std::vector<NekDouble> freeParams)
{
TimeIntegrationSchemeSharedPtr p =
MemoryManager<DIRKTimeIntegrationScheme>::AllocateSharedPtr(
"ES6", 4, freeParams);
return p;
}
static std::string className;
protected:
static std::string TimeIntegrationMethodLookupId;
}; // end class DIRKOrder4_ES6TimeIntegrationScheme
} // namespace Nektar::LibUtilities
#endif