[clib] Create common interface for Func1 base functions

This commit is contained in:
Ingmar Schoegl 2023-06-23 11:00:48 -06:00 committed by Ray Speth
parent 36faa940cd
commit 851c5b7dd7
4 changed files with 197 additions and 37 deletions

View File

@ -45,6 +45,10 @@ class Func1
public:
Func1() = default;
//! Universal constructor; interpretation of n and definition of parameters depends
//! on specialization.
Func1(size_t n, const vector<double>& params) {}
virtual ~Func1() = default;
Func1(const Func1& right);
@ -137,6 +141,9 @@ public:
m_c = omega;
}
//! Constructor uses single parameter (frequency)
Sin1(size_t n, const vector<double>& params);
Sin1(const Sin1& b) :
Func1(b) {
}
@ -179,6 +186,9 @@ public:
m_c = omega;
}
//! Constructor uses single parameter (frequency)
Cos1(size_t n, const vector<double>& params);
Cos1(const Cos1& b) :
Func1(b) {
}
@ -214,6 +224,9 @@ public:
m_c = A;
}
//! Constructor uses single parameter (exponent factor)
Exp1(size_t n, const vector<double>& params);
Exp1(const Exp1& b) :
Func1(b) {
}
@ -247,6 +260,9 @@ public:
m_c = n;
}
//! Constructor uses single parameter (exponent)
Pow1(size_t n, const vector<double>& params);
Pow1(const Pow1& b) :
Func1(b) {
}
@ -320,6 +336,9 @@ public:
m_c = A;
}
//! Constructor uses single parameter (constant)
Const1(size_t n, const vector<double>& params);
Const1(const Const1& b) :
Func1(b) {
}
@ -808,22 +827,27 @@ public:
* @param A peak value
* @param t0 offset
* @param fwhm full width at half max
* @since New in Cantera 3.0.
*/
class Gaussian : public Func1
class Gaussian1 : public Func1
{
public:
Gaussian(double A, double t0, double fwhm) {
Gaussian1(double A, double t0, double fwhm) {
m_A = A;
m_t0 = t0;
m_tau = fwhm/(2.0*std::sqrt(std::log(2.0)));
}
Gaussian(const Gaussian& b) :
//! Constructor uses 3 parameters in the following order:
//! [A, t0, fwhm]
Gaussian1(size_t n, const vector<double>& params);
Gaussian1(const Gaussian1& b) :
Func1(b) {
*this = Gaussian::operator=(b);
*this = Gaussian1::operator=(b);
}
Gaussian& operator=(const Gaussian& right) {
Gaussian1& operator=(const Gaussian1& right) {
if (&right == this) {
return *this;
}
@ -836,7 +860,7 @@ public:
}
virtual Func1& duplicate() const {
Gaussian* np = new Gaussian(*this);
Gaussian1* np = new Gaussian1(*this);
return *((Func1*)np);
}
@ -850,6 +874,25 @@ protected:
};
/**
* A Gaussian.
* \f[
* f(t) = A e^{-[(t - t_0)/\tau]^2}
* \f]
* where \f[ \tau = \frac{fwhm}{2\sqrt{\ln 2}} \f]
* @param A peak value
* @param t0 offset
* @param fwhm full width at half max
* @deprecated To be removed after Cantera 3.0; replaced by Gaussian1.
*/
class Gaussian : public Gaussian1
{
Gaussian(double A, double t0, double fwhm);
Gaussian(const Gaussian& b);
};
/**
* Polynomial of degree n.
*/
@ -861,6 +904,10 @@ public:
std::copy(c, c+m_cpoly.size(), m_cpoly.begin());
}
//! Constructor uses n + 1 parameters in the following order:
//! [an, ..., a1, a0]
Poly1(size_t n, const vector<double>& params);
Poly1(const Poly1& b) :
Func1(b) {
*this = Poly1::operator=(b);
@ -915,6 +962,10 @@ public:
std::copy(b, b+n, m_csin.begin());
}
//! Constructor uses 2 * n + 2 parameters in the following order:
//! [a0, a1, ... an, omega, b1, ... bn]
Fourier1(size_t n, const vector<double>& params);
Fourier1(const Fourier1& b) :
Func1(b) {
*this = Fourier1::operator=(b);
@ -976,6 +1027,10 @@ public:
}
}
//! Constructor uses 3 * n parameters in the following order:
//! [A1, b1, E1, A2, b2, E2, ... An, bn, En]
Arrhenius1(size_t n, const vector<double>& params);
Arrhenius1(const Arrhenius1& b) :
Func1() {
*this = Arrhenius1::operator=(b);

View File

@ -90,7 +90,7 @@ void runexample()
double A = 0.1;
double FWHM = 0.2;
double t0 = 0.5;
Gaussian igniter_mdot(A, t0, FWHM);
Gaussian1 igniter_mdot(A, t0, FWHM);
MassFlowController m3;
m3.install(igniter, combustor);

View File

@ -31,44 +31,36 @@ extern "C" {
func_t* r=0;
size_t m = lenp;
if (type == SinFuncType) {
r = new Sin1(params[0]);
vector<double> par = {params[0]};
r = new Sin1(0, par);
} else if (type == CosFuncType) {
r = new Cos1(params[0]);
vector<double> par = {params[0]};
r = new Cos1(0, par);
} else if (type == ExpFuncType) {
r = new Exp1(params[0]);
vector<double> par = {params[0]};
r = new Exp1(0, par);
} else if (type == PowFuncType) {
if (lenp < 1) {
throw CanteraError("func_new",
"exponent for pow must be supplied");
}
r = new Pow1(params[0]);
vector<double> par = {params[0]};
r = new Pow1(0, par);
} else if (type == ConstFuncType) {
r = new Const1(params[0]);
vector<double> par = {params[0]};
r = new Const1(0, par);
} else if (type == FourierFuncType) {
if (lenp < 2*n + 2) {
throw CanteraError("func_new",
"not enough Fourier coefficients");
}
r = new Fourier1(n, params[n+1], params[0], params + 1,
params + n + 2);
vector<double> par(lenp);
std::copy(params, params + lenp, par.data());
r = new Fourier1(n, par);
} else if (type == GaussianFuncType) {
if (lenp < 3) {
throw CanteraError("func_new",
"not enough Gaussian coefficients");
}
r = new Gaussian(params[0], params[1], params[2]);
vector<double> par(lenp);
std::copy(params, params + lenp, par.data());
r = new Gaussian1(n, par);
} else if (type == PolyFuncType) {
if (lenp < n + 1) {
throw CanteraError("func_new",
"not enough polynomial coefficients");
}
r = new Poly1(n, params);
vector<double> par(lenp);
std::copy(params, params + lenp, par.data());
r = new Poly1(n, par);
} else if (type == ArrheniusFuncType) {
if (lenp < 3*n) {
throw CanteraError("func_new",
"not enough Arrhenius coefficients");
}
r = new Arrhenius1(n, params);
vector<double> par(lenp);
std::copy(params, params + lenp, par.data());
r = new Arrhenius1(n, par);
} else if (type == PeriodicFuncType) {
r = new Periodic1(FuncCabinet::item(n), params[0]);
} else if (type == SumFuncType) {

View File

@ -4,6 +4,8 @@
// at https://cantera.org/license.txt for license and copyright information.
#include "cantera/numerics/Func1.h"
#include "cantera/base/global.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/base/stringUtils.h"
using namespace std;
@ -135,6 +137,15 @@ void Func1::setParent(Func1* p)
/*****************************************************************************/
Sin1::Sin1(size_t n, const vector<double>& params)
{
if (params.size() != 1) {
throw CanteraError("Sin1::Sin1",
"Constructor needs exactly one parameter (frequency).");
}
m_c = params[0];
}
string Sin1::write(const string& arg) const
{
if (m_c == 1.0) {
@ -152,6 +163,15 @@ Func1& Sin1::derivative() const
}
/*****************************************************************************/
Cos1::Cos1(size_t n, const vector<double>& params)
{
if (params.size() != 1) {
throw CanteraError("Cos1::Cos1",
"Constructor needs exactly one parameter (frequency).");
}
m_c = params[0];
}
Func1& Cos1::derivative() const
{
Func1* s = new Sin1(m_c);
@ -170,6 +190,15 @@ std::string Cos1::write(const std::string& arg) const
/**************************************************************************/
Exp1::Exp1(size_t n, const vector<double>& params)
{
if (params.size() != 1) {
throw CanteraError("Exp1::Exp1",
"Constructor needs exactly one parameter (exponent factor).");
}
m_c = params[0];
}
Func1& Exp1::derivative() const
{
Func1* f = new Exp1(m_c);
@ -191,6 +220,15 @@ std::string Exp1::write(const std::string& arg) const
/******************************************************************************/
Pow1::Pow1(size_t n, const vector<double>& params)
{
if (params.size() != 1) {
throw CanteraError("Pow1::Pow1",
"Constructor needs exactly one parameter (exponent).");
}
m_c = params[0];
}
Func1& Pow1::derivative() const
{
Func1* r;
@ -207,6 +245,69 @@ Func1& Pow1::derivative() const
/******************************************************************************/
Const1::Const1(size_t n, const vector<double>& params)
{
if (params.size() != 1) {
throw CanteraError("Const1::Const1",
"Constructor needs exactly one parameter (constant).");
}
m_c = params[0];
}
Poly1::Poly1(size_t n, const vector<double>& params)
{
if (params.size() != n + 1) {
throw CanteraError("Poly1::Poly1",
"Constructor needs exactly n + 1 = {} parameters (with n={}).", n + 1, n);
}
m_cpoly.resize(n + 1);
copy(params.data(), params.data() + m_cpoly.size(), m_cpoly.begin());
}
Fourier1::Fourier1(size_t n, const vector<double>& params)
{
if (params.size() != 2 * n + 2) {
throw CanteraError("Fourier1::Fourier1",
"Constructor needs exactly 2 * n + 2 = {} parameters (with n={}).",
2 * n + 2, n);
}
m_omega = params[n + 1];
m_a0_2 = 0.5 * params[0];
m_ccos.resize(n);
m_csin.resize(n);
copy(params.data() + 1, params.data() + n + 1, m_ccos.begin());
copy(params.data() + n + 2, params.data() + 2 * n + 2, m_csin.begin());
}
Gaussian1::Gaussian1(size_t n, const vector<double>& params)
{
if (params.size() != 3) {
throw CanteraError("Gaussian1::Gaussian1",
"Constructor needs exactly 3 parameters (amplitude, center, width).");
}
m_A = params[0];
m_t0 = params[1];
m_tau = params[2] / (2. * sqrt(log(2.)));
}
Arrhenius1::Arrhenius1(size_t n, const vector<double>& params)
{
if (params.size() != 3 * n) {
throw CanteraError("Arrhenius1::Arrhenius1",
"Constructor needs exactly 3 * n parameters grouped as (Ai, bi, Ei) for "
"i=0..n-1.");
}
m_A.resize(n);
m_b.resize(n);
m_E.resize(n);
for (size_t i = 0; i < n; i++) {
size_t loc = 3 * i;
m_A[i] = params[loc];
m_b[i] = params[loc + 1];
m_E[i] = params[loc + 2];
}
}
Tabulated1::Tabulated1(size_t n, const double* tvals, const double* fvals,
const std::string& method)
{
@ -282,6 +383,18 @@ Func1& Tabulated1::derivative() const {
/******************************************************************************/
Gaussian::Gaussian(double A, double t0, double fwhm) : Gaussian1(A, t0, fwhm)
{
warn_deprecated("Gaussian::Gaussian", "To be removed after Cantera 3.0. "
"Replaced by 'Gaussian1'.");
}
Gaussian::Gaussian(const Gaussian& b) : Gaussian1(b)
{
warn_deprecated("Gaussian::Gaussian", "To be removed after Cantera 3.0. "
"Replaced by 'Gaussian1'.");
}
string Func1::write(const std::string& arg) const
{
return fmt::format("<unknown {}>({})", ID(), arg);