anne/tasmet
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Little further

Browse Source
This commit is contained in:
Anne de Jong 2017-01-08 20:55:03 +01:00
parent 7edd65fa0d
commit 696ba04422
29 changed files with 1290 additions and 189 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -21,3 +21,4 @@ doc/usg.pdf
brent_test
tasmet
tasmet_automoc.dir
tasmet_solvemodel

7
CMakeLists.txt
View File

@ -126,11 +126,14 @@ include_directories(
src/material
src/sys
)
link_directories(/usr/lib/chaiscript)
# Add the code subdirectory
add_subdirectory(src)
add_subdirectory(testing)
add_executable(tasmet src/main.cpp src/gui/tasmet_resources.qrc)
add_executable(tasmet_solvemodel src/tasmet_solvemodel.cpp)
target_link_libraries(tasmet tasmet_gui tasmet_src messages PythonQt Qt5::Widgets openblas tcmalloc)
target_link_libraries(tasmet_solvemodel tasmet_src messages PythonQt openblas)
target_link_libraries(tasmet tasmet_gui tasmet_src messages PythonQt Qt5::Widgets chaiscript_stdlib-5.7.0 openblas)
target_link_libraries(tasmet_solvemodel tasmet_src messages PythonQt openblas chaiscript_stdlib-5.7.0)

2
Doxyfile
View File

@ -187,7 +187,7 @@ SHORT_NAMES = NO
# description.)
# The default value is: NO.
JAVADOC_AUTOBRIEF = NO
JAVADOC_AUTOBRIEF = YES
# If the QT_AUTOBRIEF tag is set to YES then doxygen will interpret the first
# line (until the first dot) of a Qt-style comment as the brief description. If

13
doc/usg.lyx
View File

@ -439,19 +439,6 @@ Solver
object.
Often, the user creates a Python scripts which has approximately the following
shape:
\begin_inset ERT
status open
\begin_layout Plain Layout
\backslash
lstinputlisting{tutorial/program_overview.py}
\end_layout
\end_inset
\end_layout
\begin_layout Standard

6
src/CMakeLists.txt
View File

@ -27,10 +27,12 @@ include_directories(
)
add_library(tasmet_src
chaiscript.cpp
tasmet_tracer.cpp
tasmet_exception.cpp
tasmet_assert.cpp
tasmet_pyeval.cpp
tasmet_evalscript.cpp
protobuf/message_tools.cpp
@ -72,4 +74,4 @@ add_library(tasmet_src
add_subdirectory(gui)
add_subdirectory(protobuf)
target_link_libraries(tasmet_src Qt5::Widgets)
target_link_libraries(tasmet_src Qt5::Widgets superlu dl pthread)

28
src/chaiscript.cpp Normal file
View File

@ -0,0 +1,28 @@
// chaiscript.cpp
//
// last-edit-by: J.A. de Jong
//
// Description:
//
//////////////////////////////////////////////////////////////////////
#include "chaiscript.h"
#include <chaiscript/chaiscript_stdlib.hpp>
#include <chaiscript/chaiscript.hpp>
#include "chaiscript/math.hpp"
#include "tasmet_types.h"
using chaiscript::ChaiScript;
std::unique_ptr<chaiscript::ChaiScript> getChaiScriptInstance() {
auto mathlib = chaiscript::extras::math::bootstrap();
std::unique_ptr<ChaiScript> chai(new ChaiScript());
// std::unique_ptr<ChaiScript> chai(new ChaiScript(chaiscript::Std_Lib::library()));
chai->add(mathlib);
chai->add_global(chaiscript::var(number_pi),"pi");
return chai;
}
//////////////////////////////////////////////////////////////////////

20
src/chaiscript.h Normal file
View File

@ -0,0 +1,20 @@
// chaiscript.h
//
// Author: J.A. de Jong
//
// Description:
// Slow compiling chaiscript therefore we provide this wrapper
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef CHAISCRIPT_H
#define CHAISCRIPT_H
namespace chaiscript {
class ChaiScript;
}
#include <memory>
std::unique_ptr<chaiscript::ChaiScript> getChaiScriptInstance();
#endif // CHAISCRIPT_H
//////////////////////////////////////////////////////////////////////

809
src/chaiscript/math.hpp Normal file
View File

@ -0,0 +1,809 @@
#include <cmath>
#include <memory>
#include <chaiscript/chaiscript.hpp>
namespace chaiscript {
namespace extras {
namespace math {
// TRIG FUNCTIONS
template<typename Ret, typename Param>
ModulePtr cos(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::cos)), "cos");
return m;
}
template<typename Ret, typename Param>
ModulePtr sin(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::sin)), "sin");
return m;
}
template<typename Ret, typename Param>
ModulePtr tan(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::tan)), "tan");
return m;
}
template<typename Ret, typename Param>
ModulePtr acos(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::acos)), "acos");
return m;
}
template<typename Ret, typename Param>
ModulePtr asin(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::asin)), "asin");
return m;
}
template<typename Ret, typename Param>
ModulePtr atan(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::atan)), "atan");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr atan2(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::atan2)), "atan2");
return m;
}
// HYPERBOLIC FUNCTIONS
template<typename Ret, typename Param>
ModulePtr cosh(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::cosh)), "cosh");
return m;
}
template<typename Ret, typename Param>
ModulePtr sinh(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::sinh)), "sinh");
return m;
}
template<typename Ret, typename Param>
ModulePtr tanh(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::tanh)), "tanh");
return m;
}
template<typename Ret, typename Param>
ModulePtr acosh(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::acosh)), "acosh");
return m;
}
template<typename Ret, typename Param>
ModulePtr asinh(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::asinh)), "asinh");
return m;
}
template<typename Ret, typename Param>
ModulePtr atanh(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::atanh)), "atanh");
return m;
}
// EXPONENTIAL AND LOGARITHMIC FUNCTIONS
template<typename Ret, typename Param>
ModulePtr exp(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::exp)), "exp");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr frexp(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::frexp)), "frexp");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr ldexp(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::ldexp)), "ldexp");
return m;
}
template<typename Ret, typename Param>
ModulePtr log(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log)), "log");
return m;
}
template<typename Ret, typename Param>
ModulePtr log10(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log10)), "log10");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr modf(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::modf)), "modf");
return m;
}
template<typename Ret, typename Param>
ModulePtr exp2(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::exp2)), "exp2");
return m;
}
template<typename Ret, typename Param>
ModulePtr expm1(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::expm1)), "expm1");
return m;
}
template<typename Ret, typename Param>
ModulePtr ilogb(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::ilogb)), "ilogb");
return m;
}
template<typename Ret, typename Param>
ModulePtr log1p(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log1p)), "log1p");
return m;
}
template<typename Ret, typename Param>
ModulePtr log2(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::log2)), "log2");
return m;
}
template<typename Ret, typename Param>
ModulePtr logb(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::logb)), "logb");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr scalbn(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::scalbn)), "scalbn");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr scalbln(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::scalbln)), "scalbln");
return m;
}
// POWER FUNCTIONS
template<typename Ret, typename Param1, typename Param2>
ModulePtr pow(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::pow)), "pow");
return m;
}
template<typename Ret, typename Param>
ModulePtr sqrt(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::sqrt)), "sqrt");
return m;
}
template<typename Ret, typename Param>
ModulePtr cbrt(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::cbrt)), "cbrt");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr hypot(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::hypot)), "hypot");
return m;
}
// ERROR AND GAMMA FUNCTIONS
template<typename Ret, typename Param>
ModulePtr erf(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::erf)), "erf");
return m;
}
template<typename Ret, typename Param>
ModulePtr erfc(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::erfc)), "erfc");
return m;
}
template<typename Ret, typename Param>
ModulePtr tgamma(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::tgamma)), "tgamma");
return m;
}
template<typename Ret, typename Param>
ModulePtr lgamma(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::lgamma)), "lgamma");
return m;
}
// ROUNDING AND REMAINDER FUNCTIONS
template<typename Ret, typename Param>
ModulePtr ceil(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::ceil)), "ceil");
return m;
}
template<typename Ret, typename Param>
ModulePtr floor(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::floor)), "floor");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr fmod(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fmod)), "fmod");
return m;
}
template<typename Ret, typename Param>
ModulePtr trunc(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::trunc)), "trunc");
return m;
}
template<typename Ret, typename Param>
ModulePtr round(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::round)), "round");
return m;
}
template<typename Ret, typename Param>
ModulePtr lround(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::lround)), "lround");
return m;
}
// long long ints do not work
template<typename Ret, typename Param>
ModulePtr llround(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::llround)), "llround");
return m;
}
template<typename Ret, typename Param>
ModulePtr rint(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::rint)), "rint");
return m;
}
template<typename Ret, typename Param>
ModulePtr lrint(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::lrint)), "lrint");
return m;
}
// long long ints do not work
template<typename Ret, typename Param>
ModulePtr llrint(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::llrint)), "llrint");
return m;
}
template<typename Ret, typename Param>
ModulePtr nearbyint(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::nearbyint)), "nearbyint");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr remainder(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::remainder)), "remainder");
return m;
}
template<typename Ret, typename Param1, typename Param2, typename Param3>
ModulePtr remquo(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2, Param3)>(&std::remquo)), "remquo");
return m;
}
// FLOATING-POINT MANIPULATION FUNCTIONS
template<typename Ret, typename Param1, typename Param2>
ModulePtr copysign(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::copysign)), "copysign");
return m;
}
template<typename Ret, typename Param>
ModulePtr nan(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::nan)), "nan");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr nextafter(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::nextafter)), "nextafter");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr nexttoward(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::nexttoward)), "nexttoward");
return m;
}
// MINIMUM, MAXIMUM, DIFFERENCE FUNCTIONS
template<typename Ret, typename Param1, typename Param2>
ModulePtr fdim(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fdim)), "fdim");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr fmax(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fmax)), "fmax");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr fmin(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::fmin)), "fmin");
return m;
}
// OTHER FUNCTIONS
template<typename Ret, typename Param>
ModulePtr fabs(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::fabs)), "fabs");
return m;
}
template<typename Ret, typename Param>
ModulePtr abs(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::abs)), "abs");
return m;
}
template<typename Ret, typename Param1, typename Param2, typename Param3>
ModulePtr fma(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2, Param3)>(&std::fma)), "fma");
return m;
}
// CLASSIFICATION FUNCTIONS
template<typename Ret, typename Param>
ModulePtr fpclassify(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::fpclassify)), "fpclassify");
return m;
}
template<typename Ret, typename Param>
ModulePtr isfinite(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isfinite)), "isfinite");
return m;
}
template<typename Ret, typename Param>
ModulePtr isinf(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isinf)), "isinf");
return m;
}
template<typename Ret, typename Param>
ModulePtr isnan(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isnan)), "isnan");
return m;
}
template<typename Ret, typename Param>
ModulePtr isnormal(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::isnormal)), "isnormal");
return m;
}
template<typename Ret, typename Param>
ModulePtr signbit(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param)>(&std::signbit)), "signbit");
return m;
}
// COMPARISON FUNCTIONS
template<typename Ret, typename Param1, typename Param2>
ModulePtr isgreater(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isgreater)), "isgreater");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr isgreaterequal(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isgreaterequal)), "isgreaterequal");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr isless(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isless)), "isless");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr islessequal(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::islessequal)), "islessequal");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr islessgreater(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::islessgreater)), "islessgreater");
return m;
}
template<typename Ret, typename Param1, typename Param2>
ModulePtr isunordered(ModulePtr m = std::make_shared<Module>())
{
m->add(chaiscript::fun(static_cast<Ret (*)(Param1, Param2)>(&std::isunordered)), "isunordered");
return m;
}
ModulePtr bootstrap(ModulePtr m = std::make_shared<Module>())
{
// TRIG FUNCTIONS
cos<double, double>(m);
cos<float, float>(m);
cos<long double, long double>(m);
sin<double, double>(m);
sin<float, float>(m);
sin<long double, long double>(m);
tan<double, double>(m);
tan<float, float>(m);
tan<long double, long double>(m);
acos<double, double>(m);
acos<float, float>(m);
acos<long double, long double>(m);
asin<double, double>(m);
asin<float, float>(m);
asin<long double, long double>(m);
atan<double, double>(m);
atan<float, float>(m);
atan<long double, long double>(m);
atan2<double, double, double>(m);
atan2<float, float, float>(m);
atan2<long double, long double, long double>(m);
// HYPERBOLIC FUNCTIONS
cosh<double, double>(m);
cosh<float, float>(m);
cosh<long double, long double>(m);
sinh<double, double>(m);
sinh<float, float>(m);
sinh<long double, long double>(m);
tanh<double, double>(m);
tanh<float, float>(m);
tanh<long double, long double>(m);
acosh<double, double>(m);
acosh<float, float>(m);
acosh<long double, long double>(m);
asinh<double, double>(m);
asinh<float, float>(m);
asinh<long double, long double>(m);
atanh<double, double>(m);
atanh<float, float>(m);
atanh<long double, long double>(m);
// EXPONENTIAL AND LOGARITHMIC FUNCTIONS
exp<double, double>(m);
exp<float, float>(m);
exp<long double, long double>(m);
frexp<double, double, int *>(m);
frexp<float, float, int *>(m);
frexp<long double, long double, int *>(m);
ldexp<double, double, int>(m);
ldexp<float, float, int>(m);
ldexp<long double, long double, int>(m);
log<double, double>(m);
log<float, float>(m);
log<long double, long double>(m);
log10<double, double>(m);
log10<float, float>(m);
log10<long double, long double>(m);
modf<double, double, double *>(m);
modf<float, float, float *>(m);
modf<long double, long double, long double *>(m);
exp2<double, double>(m);
exp2<float, float>(m);
exp2<long double, long double>(m);
expm1<double, double>(m);
expm1<float, float>(m);
expm1<long double, long double>(m);
ilogb<int, double>(m);
ilogb<int, float>(m);
ilogb<int, long double>(m);
log1p<double, double>(m);
log1p<float, float>(m);
log1p<long double, long double>(m);
log2<double, double>(m);
log2<float, float>(m);
log2<long double, long double>(m);
logb<double, double>(m);
logb<float, float>(m);
logb<long double, long double>(m);
scalbn<double, double, int>(m);
scalbn<float, float, int>(m);
scalbn<long double, long double, int>(m);
scalbln<double, double, long int>(m);
scalbln<float, float, long int>(m);
scalbln<long double, long double, long int>(m);
// POWER FUNCTIONS
pow<double, double, double>(m);
pow<float, float, float>(m);
pow<long double, long double, long double>(m);
sqrt<double, double>(m);
sqrt<float, float>(m);
sqrt<long double, long double>(m);
cbrt<double, double>(m);
cbrt<float, float>(m);
cbrt<long double, long double>(m);
hypot<double, double, double>(m);
hypot<float, float, float>(m);
hypot<long double, long double, long double>(m);
// ERROR AND GAMMA FUNCTIONS
erf<double, double>(m);
erf<float, float>(m);
erf<long double, long double>(m);
erfc<double, double>(m);
erfc<float, float>(m);
erfc<long double, long double>(m);
tgamma<double, double>(m);
tgamma<float, float>(m);
tgamma<long double, long double>(m);
lgamma<double, double>(m);
lgamma<float, float>(m);
lgamma<long double, long double>(m);
// ROUNDING AND REMAINDER FUNCTIONS
ceil<double, double>(m);
ceil<float, float>(m);
ceil<long double, long double>(m);
floor<double, double>(m);
floor<float, float>(m);
floor<long double, long double>(m);
fmod<double, double, double>(m);
fmod<float, float, float>(m);
fmod<long double, long double, long double>(m);
trunc<double, double>(m);
trunc<float, float>(m);
trunc<long double, long double>(m);
round<double, double>(m);
round<float, float>(m);
round<long double, long double>(m);
lround<long int, double>(m);
lround<long int, float>(m);
lround<long int, long double>(m);
// long long ints do not work
llround<long long int, double>(m);
llround<long long int, float>(m);
llround<long long int, long double>(m);
rint<double, double>(m);
rint<float, float>(m);
rint<long double, long double>(m);
lrint<long int, double>(m);
lrint<long int, float>(m);
lrint<long int, long double>(m);
// long long ints do not work
llrint<long long int, double>(m);
llrint<long long int, float>(m);
llrint<long long int, long double>(m);
nearbyint<double, double>(m);
nearbyint<float, float>(m);
nearbyint<long double, long double>(m);
remainder<double, double, double>(m);
remainder<float, float, float>(m);
remainder<long double, long double, long double>(m);
remquo<double, double, double, int *>(m);
remquo<float, float, float, int *>(m);
remquo<long double, long double, long double, int *>(m);
// FLOATING-POINT MANIPULATION FUNCTIONS
copysign<double, double, double>(m);
copysign<float, float, float>(m);
copysign<long double, long double, long double>(m);
nan<double, const char*>(m);
nextafter<double, double, double>(m);
nextafter<float, float, float>(m);
nextafter<long double, long double, long double>(m);
nexttoward<double, double, long double>(m);
nexttoward<float, float, long double>(m);
nexttoward<long double, long double, long double>(m);
// MINIMUM, MAXIMUM, DIFFERENCE FUNCTIONS
fdim<double, double, double>(m);
fdim<float, float, float>(m);
fdim<long double, long double, long double>(m);
fmax<double, double, double>(m);
fmax<float, float, float>(m);
fmax<long double, long double, long double>(m);
fmin<double, double, double>(m);
fmin<float, float, float>(m);
fmin<long double, long double, long double>(m);
// OTHER FUNCTIONS
fabs<double, double>(m);
fabs<float, float>(m);
fabs<long double, long double>(m);
abs<double, double>(m);
abs<float, float>(m);
abs<long double, long double>(m);
fma<double, double, double, double>(m);
fma<float, float, float, float>(m);
fma<long double, long double, long double, long double>(m);
// CLASSIFICATION FUNCTIONS
fpclassify<int, float>(m);
fpclassify<int, double>(m);
fpclassify<int, long double>(m);
isfinite<bool, float>(m);
isfinite<bool, double>(m);
isfinite<bool, long double>(m);
isinf<bool, float>(m);
isinf<bool, double>(m);
isinf<bool, long double>(m);
isnan<bool, float>(m);
isnan<bool, double>(m);
isnan<bool, long double>(m);
isnormal<bool, float>(m);
isnormal<bool, double>(m);
isnormal<bool, long double>(m);
signbit<bool, float>(m);
signbit<bool, double>(m);
signbit<bool, long double>(m);
// COMPARISON FUNCTIONS
isgreater<bool, double, double>(m);
isgreater<bool, float, float>(m);
isgreater<bool, long double, long double>(m);
isgreaterequal<bool, double, double>(m);
isgreaterequal<bool, float, float>(m);
isgreaterequal<bool, long double, long double>(m);
isless<bool, double, double>(m);
isless<bool, float, float>(m);
isless<bool, long double, long double>(m);
islessequal<bool, double, double>(m);
islessequal<bool, float, float>(m);
islessequal<bool, long double, long double>(m);
islessgreater<bool, double, double>(m);
islessgreater<bool, float, float>(m);
islessgreater<bool, long double, long double>(m);
isunordered<bool, double, double>(m);
isunordered<bool, float, float>(m);
isunordered<bool, long double, long double>(m);
return m;
}
}
}
}

61
src/duct/duct.cpp
View File

@ -8,24 +8,27 @@
#include "duct.h"
#include "tasystem.h"
#include "tasmet_assert.h"
#include "tasmet_pyeval.h"
#include "tasmet_evalscript.h"
Duct::Duct(const us id,const pb::Duct& duct):
Segment(id,duct.name()),
Geom(duct)
Duct::Duct(const us id,const pb::Duct& ductpb):
Segment(id,ductpb.name()),
Geom(ductpb),
_ductpb(ductpb)
{
TRACE(15,"Duct::Duct()");
const char* invTsfun = "Invalid solid-temperature prescribing function";
EvaluateFun Tsfun(duct.stempfunc(),invTsfun);
Tsfun.addGlobalDef("L",duct.length());
EvaluateFun Tsfun(ductpb.stempfunc(),invTsfun);
Tsfun.addGlobalDef("L",ductpb.length());
_Tsprescribed = Tsfun(x);
if(min(_Tsprescribed) < constants::min_T0 ||
max(_Tsprescribed) > constants::max_T0) {
throw TaSMETError(invTsfun);
}
switch (duct.htmodel()) {
switch (ductpb.htmodel()) {
case pb::Isentropic: {
break;
}
@ -34,25 +37,31 @@ Duct::Duct(const us id,const pb::Duct& duct):
break;
}
}
Duct::Duct(const Duct& other):
Segment(other),
Geom(other),
_ductpb(other._ductpb),
_Tsprescribed(other._Tsprescribed)
{
// Do something with the equations here
TRACE(15,"Duct::~Duct");
}
Duct* Duct::copy() const {
return new Duct(*this);
}
Duct::~Duct() {
TRACE(15,"Duct::~Duct");
// for(Equation* eq: _eqs){
// delete eq;
// }
}
void Duct::residual(const TaSystem& sys,arma::subview_col<d> && residual) const {
TRACE(15,"Duct::residual()");
const arma::subview_col<d> sol = sys.getSolution(_id);
VARTRACE(15,sol(0));
@ -61,16 +70,22 @@ void Duct::residual(const TaSystem& sys,arma::subview_col<d> && residual) const
}
vd Duct::initialSolution(const TaSystem& sys) const {
vd initsol(getNDofs(sys));
TRACE(15,"Duct::initialSolution()");
vd initsol(getNDofs(sys));
VARTRACE(15,initsol.size());
us vars_per_gp = 4;
vars_per_gp+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
us vars_per_gp = 5;
us Ns = sys.Ns();
const Gas& gas = sys.gas();
for(us i=0;i<ngp();i++){
VARTRACE(15,i);
// Initial density
initsol.subvec((i*vars_per_gp+0)*Ns,(i*vars_per_gp+1)*Ns-1) =
gas.rho0();
@ -79,17 +94,18 @@ vd Duct::initialSolution(const TaSystem& sys) const {
0;
// Initial Temperature
initsol.subvec((i*vars_per_gp+2)*Ns,(i*vars_per_gp+2)*Ns-1) =
initsol.subvec((i*vars_per_gp+2)*Ns,(i*vars_per_gp+3)*Ns-1) =
_Tsprescribed(i);
// Initial pressure
initsol.subvec((i*vars_per_gp+3)*Ns,(i*vars_per_gp+4)*Ns-1) =
gas.p0();
// Initial solid temperature
// Initial solid temperature, if not prescribed
if(_ductpb.stempmodel() != pb::Prescribed) {
initsol.subvec((i*vars_per_gp+4)*Ns,(i*vars_per_gp+5)*Ns-1) =
_Tsprescribed(i);
}
}
return initsol;
@ -97,18 +113,23 @@ vd Duct::initialSolution(const TaSystem& sys) const {
}
us Duct::getNEqs(const TaSystem& sys) const {
TRACE(15,"Duct::getNEqs()");
us Ns = sys.Ns();
us number_eqs = _eqs.size();
us number_eqs = 4;
// When we have to solve a solid heat balance
number_eqs+= (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
return Ns*number_eqs*(ngp()-1);
return Ns*number_eqs*ngp();
}
us Duct::getNDofs(const TaSystem& sys) const {
TRACE(15,"Duct::getNDofs()");
us Ns = sys.Ns();
// rho,u,T,p,Ts
us nvars_per_gp = 5;
us nvars_per_gp = 4;
nvars_per_gp += (_ductpb.stempmodel() == pb::HeatBalance ? 1 : 0);
return Ns*nvars_per_gp*ngp();
}

5
src/duct/duct.h
View File

@ -22,8 +22,11 @@ class Duct : public Segment, public Geom {
// Drag* _drag = nullptr;
// Heat* _heat = nullptr;
std::vector<Equation*> _eqs;
pb::Duct _ductpb;
vd _Tsprescribed;
protected:
Duct(const Duct&);
public:
@ -38,6 +41,8 @@ public:
virtual Duct* copy() const;
const Geom& geom() const;
const pb::Duct& getDuctPb() const { return _ductpb;}
// Solving
virtual void residual(const TaSystem&,arma::subview_col<d>&& residual) const;

20
src/duct/geom.cpp
View File

@ -8,12 +8,14 @@
#include "geom.h"
#include "grid.h"
#include "duct.pb.h"
#include "tasmet_pyeval.h"
#include "tasmet_evalscript.h"
#include <memory>
#include "tasmet_exception.h"
#include "tasmet_tracer.h"
Geom::Geom(const pb::Duct& duct) {
TRACE(15,"Geom::Geom");
// Step one: create the grid
std::unique_ptr<Grid> grid;
switch (duct.gridtype()) {
@ -35,28 +37,32 @@ Geom::Geom(const pb::Duct& duct) {
x = grid->getx();
const char* invS = "Invalid cross-sectional area function";
{
EvaluateFun Sfun(duct.area(),invS);
Sfun.addGlobalDef("L",duct.length());
S = Sfun(x);
}
if(min(S) <= 0) {
throw TaSMETError(invS);
}
const char* invpor = "Invalid porosity function";
{
EvaluateFun phifun(duct.phi(),invpor);
Sfun.addGlobalDef("L",duct.length());
phi = Sfun(x);
phifun.addGlobalDef("L",duct.length());
phi = phifun(x);
}
if(min(phi) <= 0 || max(phi) > 1) {
throw TaSMETError(invpor);
}
const char* invrh = "Invalid hydraulic radius function";
{
EvaluateFun rhfun(duct.rh(),invrh);
Sfun.addGlobalDef("L",duct.length());
rh = Sfun(x);
rhfun.addGlobalDef("L",duct.length());
rh = rhfun(x);
}
if(min(rh) <= 0 ) {
throw TaSMETError(invrh);
}

24
src/ductbc/adiabaticwall.cpp
View File

@ -9,14 +9,14 @@
#include "ductbc.pb.h"
#include "tasmet_tracer.h"
#include "tasmet_assert.h"
#include "tasystem.h"
#include "duct.h"
AdiabaticWall::AdiabaticWall(const us id,
const TaSystem& sys,
const pb::DuctBc& dbc):
Segment(id,dbc.name()),
_side(dbc.side()),
_duct_id(dbc.duct_id())
DuctBc(id,dbc)
{
TRACE(15,"AdiabaticWall(id,sys,dbc)");
tasmet_assert(dbc.type() == pb::AdiabaticWall,"Wrong type given to constructor");
@ -24,8 +24,8 @@ AdiabaticWall::AdiabaticWall(const us id,
}
AdiabaticWall::AdiabaticWall(const AdiabaticWall& o):
Segment(o) {
DuctBc(o)
{
TRACE(15,"AdiabaticWall(o)");
}
@ -39,19 +39,21 @@ vd AdiabaticWall::initialSolution(const TaSystem& sys) const {
return vd();
}
void AdiabaticWall::residual(const TaSystem&,
void AdiabaticWall::residual(const TaSystem& sys,
arma::subview_col<d>&& residual
) const {
TRACE(15,"AdiabaticWall::residual()");
}
us AdiabaticWall::getNEqs(const TaSystem&) const {
us AdiabaticWall::getNEqs(const TaSystem& sys) const {
TRACE(15,"AdiabaticWall::getNEqs()");
// u = 0
// dT/dx = 0
// dTs/dx = 0 => 3 eqs
bool has_solideq = getDuct(sys).getDuctPb().stempmodel() != pb::Prescribed;
return sys.Ns()*(has_solideq ? 3: 2);
}
void AdiabaticWall::show(const TaSystem&,us verbosity_level) const {
TRACE(15,"AdiabaticWall::show()");

4
src/ductbc/adiabaticwall.h
View File

@ -9,7 +9,7 @@
#ifndef ADIABATICWALL_H
#define ADIABATICWALL_H
#include "segment.h"
#include "ductbc.pb.h"
#include "ductbc.h"
class TaSystem;
class Variable;
@ -19,7 +19,7 @@ class Variable;
* blocks the flow and does not allow any axial heat conduction.
*/
class AdiabaticWall: public Segment {
class AdiabaticWall: public DuctBc {
pb::DuctSide _side; /**< Duct side at which this b.c. works */
us _duct_id; /**< ID of Duct for this b.c. */

10
src/ductbc/ductbc.cpp
View File

@ -11,8 +11,16 @@
#include "tasmet_assert.h"
#include "pressurebc.h"
#include "adiabaticwall.h"
#include "tasystem.h"
#include "duct.h"
Segment* DuctBc::newDuctBc(const us id,
const Duct& DuctBc::getDuct(const TaSystem& sys) const {
return dynamic_cast<const Duct&>(sys.getSegment(_dbc.duct_id()));
}
DuctBc* DuctBc::newDuctBc(const us id,
const TaSystem& sys,
const pb::DuctBc& dbc) {

22
src/ductbc/ductbc.h
View File

@ -9,17 +9,27 @@
#ifndef DUCTBC_H
#define DUCTBC_H
#include "segment.h"
#include "ductbc.pb.h"
namespace pb{
class DuctBc;
}
class DuctBc {
class TaSystem;
class Duct;
class DuctBc :public Segment {
pb::DuctBc _dbc;
public:
static Segment* newDuctBc(const us id,
DuctBc(const us id,
const pb::DuctBc& dbc):
Segment(id,dbc.name()),
_dbc(dbc) {}
DuctBc(const DuctBc& o): Segment(o),_dbc(o._dbc) {}
static DuctBc* newDuctBc(const us id,
const TaSystem& sys,
const pb::DuctBc&);
const Duct& getDuct(const TaSystem&) const;
};

19
src/ductbc/pressurebc.cpp
View File

@ -9,20 +9,19 @@
#include "pressurebc.h"
#include "ductbc.pb.h"
#include "tasmet_tracer.h"
#include "tasystem.h"
#include "duct.h"
PressureBc::PressureBc(const us id,
const TaSystem& sys,
const pb::DuctBc& dbc):
Segment(id,dbc.name())
DuctBc(id,dbc)
{
TRACE(15,"PressureBc(id,sys,dbc)");
}
PressureBc::PressureBc(const PressureBc& o):
Segment(o) {
DuctBc(o) {
TRACE(15,"PressureBc(o)");
@ -46,10 +45,14 @@ void PressureBc::residual(const TaSystem&,
}
us PressureBc::getNEqs(const TaSystem&) const {
us PressureBc::getNEqs(const TaSystem& sys) const {
TRACE(15,"PressureBc::getNEqs()");
// p = x
// T = x
// This one only if the duct solves for solid
// Ts = x => 3 equations
bool has_solideq = getDuct(sys).getDuctPb().stempmodel() != pb::Prescribed;
return sys.Ns()*(has_solideq ? 3: 2);
}
void PressureBc::show(const TaSystem&,us verbosity_level) const {
TRACE(15,"PressureBc::show()");

7
src/ductbc/pressurebc.h
View File

@ -8,14 +8,15 @@
#pragma once
#ifndef PRESSUREBC_H
#define PRESSUREBC_H
#include "segment.h"
#include "ductbc.h"
#include "ductbc.pb.h"
class TaSystem;
class Variable;
class PressureBc: public Segment {
class PressureBc: public DuctBc {
Variable *_p,*_T,*_Ts;
us _duct_id; /**< ID of Duct for this b.c. */
pb::DuctSide _side; /**< Duct side at which this b.c. works */
protected:
PressureBc(const PressureBc&);

67
src/gui/mainwindow.ui
View File

@ -1,5 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<author>J.A. de Jong</author>
<class>MainWindow</class>
<widget class="QMainWindow" name="MainWindow">
<property name="geometry">
@ -7,7 +8,7 @@
<x>0</x>
<y>0</y>
<width>683</width>
<height>488</height>
<height>554</height>
</rect>
</property>
<property name="windowTitle">
@ -63,21 +64,30 @@
<item row="5" column="0">
<widget class="QLabel" name="label_7">
<property name="text">
<string>System type:</string>
<string>&amp;System type:</string>
</property>
<property name="buddy">
<cstring>systemtype</cstring>
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QLabel" name="label_5">
<property name="text">
<string>Gas type:</string>
<string>&amp;Gas type:</string>
</property>
<property name="buddy">
<cstring>gastype</cstring>
</property>
</widget>
</item>
<item row="3" column="0">
<widget class="QLabel" name="label_15">
<property name="text">
<string>Reference temperature</string>
<string>&amp;Reference temperature</string>
</property>
<property name="buddy">
<cstring>T0</cstring>
</property>
</widget>
</item>
@ -114,14 +124,20 @@
<item row="4" column="0">
<widget class="QLabel" name="label_17">
<property name="text">
<string>Reference pressure</string>
<string>&amp;Reference pressure</string>
</property>
<property name="buddy">
<cstring>p0</cstring>
</property>
</widget>
</item>
<item row="0" column="0">
<widget class="QLabel" name="label">
<property name="text">
<string>Number of harmonics: </string>
<string>&amp;Number of harmonics: </string>
</property>
<property name="buddy">
<cstring>nf</cstring>
</property>
</widget>
</item>
@ -135,7 +151,10 @@
<item row="1" column="0">
<widget class="QLabel" name="label_4">
<property name="text">
<string>Fundamental frequency: </string>
<string>&amp;Fundamental frequency: </string>
</property>
<property name="buddy">
<cstring>freq</cstring>
</property>
</widget>
</item>
@ -192,7 +211,10 @@
<item row="0" column="0">
<widget class="QLabel" name="label_9">
<property name="text">
<string>Segment type</string>
<string>S&amp;egment type</string>
</property>
<property name="buddy">
<cstring>segmenttype</cstring>
</property>
</widget>
</item>
@ -202,7 +224,10 @@
<item row="1" column="0">
<widget class="QLabel" name="label_10">
<property name="text">
<string>Segment ID</string>
<string>Segment &amp;ID</string>
</property>
<property name="buddy">
<cstring>segmentid</cstring>
</property>
</widget>
</item>
@ -216,7 +241,10 @@
<item row="2" column="0">
<widget class="QLabel" name="label_8">
<property name="text">
<string>Segment name</string>
<string>Segment &amp;name</string>
</property>
<property name="buddy">
<cstring>segmentname</cstring>
</property>
</widget>
</item>
@ -237,7 +265,7 @@
<item row="3" column="0" colspan="2">
<widget class="QPushButton" name="addsegment">
<property name="text">
<string>Add/edit segment...</string>
<string>&amp;Add/edit segment...</string>
</property>
</widget>
</item>
@ -320,7 +348,7 @@
<x>0</x>
<y>0</y>
<width>683</width>
<height>18</height>
<height>21</height>
</rect>
</property>
<widget class="QMenu" name="menuFile">
@ -432,6 +460,21 @@
</property>
</action>
</widget>
<tabstops>
<tabstop>nf</tabstop>
<tabstop>freq</tabstop>
<tabstop>gastype</tabstop>
<tabstop>T0</tabstop>
<tabstop>p0</tabstop>
<tabstop>systemtype</tabstop>
<tabstop>segmenttype</tabstop>
<tabstop>segmentid</tabstop>
<tabstop>segmentname</tabstop>
<tabstop>addsegment</tabstop>
<tabstop>removesegment</tabstop>
<tabstop>backlog</tabstop>
<tabstop>segoverview</tabstop>
</tabstops>
<resources>
<include location="tasmet_resources.qrc"/>
</resources>

18
src/main.cpp
View File

@ -9,10 +9,6 @@
#include <QCommandLineParser>
#include <QCommandLineOption>
// For using python from within Qt
#include <PythonQt.h>
#include "tasmet_config.h"
#include "tasmet_tracer.h"
#include "gui/mainwindow.h"
@ -45,18 +41,6 @@ int main(int argc, char *argv[]) {
// Q_INIT_RESOURCE(application);
INITTRACE(15);
// Initialize PythonQt
// PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut);
PythonQt::init();
PythonQt* pyqt = PythonQt::self();
PythonQtObjectPtr context = pyqt->getMainModule();
QVariant rv = context.evalScript("from math import *\n");
if(pyqt->hadError()) {
return -1;
}
QApplication app(argc, argv);
catchUnixSignals({SIGQUIT, SIGINT, SIGTERM, SIGHUP});
@ -74,8 +58,6 @@ int main(int argc, char *argv[]) {
// if (!parser.positionalArguments().isEmpty())
// mainWin.loadFile(parser.positionalArguments().first());
TaSMETMainWindow win;
win.show();

29
src/protobuf/message_tools.cpp
View File

@ -12,9 +12,9 @@
#include "tasmet_exception.h"
#include "model.pb.h"
#include <google/protobuf/text_format.h>
using google::protobuf::TextFormat;
#include <google/protobuf/util/message_differencer.h>
// #include <google/protobuf/text_format.h>
// using google::protobuf::TextFormat;
template<typename T>
T loadMessage(const string& filepath) {
@ -32,11 +32,11 @@ T loadMessage(const string& filepath) {
T sys;
std::stringstream strStream;
strStream << myfile.rdbuf(); //read the file
string data = strStream.str();//str holds the content of the file
// std::stringstream strStream;
// strStream << myfile.rdbuf(); //read the file
// string data = strStream.str();//str holds the content of the file
if(!TextFormat::ParseFromString(data,&sys)) {
if(!sys.ParseFromIstream(&myfile)) {
string error = "Invalid TaSMET Model file: ";
error += filepath;
throw TaSMETError(error);
@ -53,27 +53,20 @@ void saveMessage(const string& filepath,const T& sys) {
if(!sfile.good()){
throw TaSMETError("Could not write to file");
}
string data;
if(TextFormat::PrintToString(sys,&data)) {
// Can maybe assign to itself. Which is no problem
// according to C++ docs
TRACE(15,"Saving file succeeded");
sfile << data << std::flush;
// Close file here, such that in can be opened to compare
// whether the file is still dirty
sfile.close();
// if(TextFormat::PrintToString(sys,&data)) {
if(sys.SerializeToOstream(&sfile)) {
TRACE(15,"Saving file succeeded");
}
else {
throw TaSMETError("Could not save model to file");
}
}
// // Returns true when the two messages are equal
template<typename T>
bool compareMessage(const T& s1,const T& s2) {
return (s1.SerializeAsString()==s2.SerializeAsString());
return google::protobuf::util::MessageDifferencer::Equals(s1,s2);
}
// Explicit instantiation for pb::Model

4
src/protobuf/model.proto
View File

@ -8,11 +8,7 @@ message Model {
optional System system = 1;
optional SolverParams sparams = 2;
optional string backlog = 3 [default="Type your info here..."];
repeated double solution = 4;
}

14
src/solver/newton_raphson.h
View File

@ -10,16 +10,18 @@
#define NEWTON_RAPHSON_H
#include "solver.h"
/**
* Newton Raphson Solver implementation.
*
*/
class NewtonRaphson: public Solver<GradientNonlinearSystem,vd> {
d _dampfac = 1.0;
d _maxiter = 10000;
d _dampfac = 1.0; /**< This is the applied damping factor */
public:
NewtonRaphson(const GradientNonlinearSystem&sys,d dampfac,us maxiter):
NewtonRaphson(const GradientNonlinearSystem&sys,d dampfac=1.0):
Solver(sys),
_dampfac(dampfac),
_maxiter(maxiter){}
_dampfac(dampfac)
{}
protected:
void start_implementation(GradientNonlinearSystem& system,progress_callback* callback);

42
src/solver/solver.cpp
View File

@ -8,13 +8,32 @@
#include "tasmet_tracer.h"
#include "solver.h"
#include "tasmet_exception.h"
#include <chrono>
#include "tasmet_io.h"
SolverAction ostream_progress_callback(SolverProgress pg,std::ostream* out, d funtol,d reltol) {
// Create reference to ostream
std::ostream& myout = (out) ? *out : std::cout;
myout << "Iteration: " << pg.iteration;
myout << " , Residual: " << pg.fun_err << " . Solution error: ";
myout << pg.rel_err << endl;
if(pg.fun_err <= funtol && pg.rel_err <= reltol) {
myout << "Solution found within required tolerance" << endl;
return Stop;
}
return Continue;
}
template<typename system_T,typename result_T>
static void SolverThread(Solver<system_T,result_T>* solver,
system_T* system,
progress_callback* callback);
progress_callback simple_progress_callback(d funtol,d reltol,std::ostream* out) {
using namespace std::placeholders; // for _1, _2, _3...
using std::bind;
return bind(ostream_progress_callback,_1,out,funtol,reltol);
}
template<typename system_T,typename result_T>
@ -33,28 +52,15 @@ Solver<system_T,result_T>::~Solver(){
template<typename system_T,typename result_T>
void Solver<system_T,result_T>::start(progress_callback* callback){
TRACE(15,"Waiting for solver...");
start_implementation(*_sys,callback);
}
template<typename system_T,typename result_T>
result_T Solver<system_T,result_T>::getSolution() {
return _sys->getSolution();
}
template<typename system_T,typename result_T>
void SolverThread(Solver<system_T,result_T>* solver,system_T* system,progress_callback* callback) {
assert(system);
solver->_running = true;
solver->start_implementation(*system,callback);
solver->_running = false;
}
// Explicit instantiation for some types of systems and results

58
src/solver/solver.h
View File

@ -13,44 +13,88 @@
#include <functional>
#include <thread>
#include <atomic>
#include <ostream>
#include "tasmet_types.h"
#include "tasmet_assert.h"
#include "system.h"
/**
* Passed to the progress_callback function to keep track of the
* Solver progress.
*
*/
struct SolverProgress
{
size_t iteration = 0;
size_t iteration = 0; /**< The current iteration no. */
d fun_err = 1e0;
d rel_err = 1e0;
bool done = false;
bool error = false; /**< Tells the callback that the
Solver stops as an action of
itself. Probably due to an internal
error. */
};
/**
* This enum should be return by the Solver callback function.
*
*/
enum SolverAction{
Continue=0,
Stop=1
Continue=0, /**< Tells the Solver to continue */
Stop=1 /**< Tells the Solver to stop. */
};
namespace pb{
class SolverParams;
}
/// Typedef for the used callback which is called by a Solver instance.
typedef std::function<SolverAction(SolverProgress)> progress_callback;
/*!
* Returns a progress calback function that prints the progress
* to the given ostream. This is a simple CLI callback function which
* prints the solver progress to the given ostream as input. If no
* ostream is given, the callback will print to stdout.
*
* @param funtol The required residual tolerance to stop the solver
* @param reltol The required solution tolerance to stop the solver
* @param out The ostream to output the progress to.
*
* @return A std::function object that can be passed to Solver::Start().
*/
progress_callback simple_progress_callback(d funtol,d reltol,
std::ostream* out=nullptr);
template<typename system_T,typename result_T>
class Solver {
protected:
system_T* _sys = nullptr;
system_T* _sys = nullptr; /**< Stores and owns a system (creates
a copy from the given System in the
constructor */
public:
Solver(const system_T& sys);
Solver(const Solver&)=delete;
Solver& operator=(const Solver&)=delete;
/// Start the solver, using the given progress callback
void start(progress_callback* callback);
// Returns the solution of the problem
/// Returns the solution of the problem after the Solver is done
/// solving.
result_T getSolution();
virtual ~Solver();
/// Create a new Solver instance based on protobuf
/// parameters. Throws a TaSMETError on error.
static Solver* newSolver(const pb::SolverParams& params);
protected:
/// This member fcn should be implemented by the Solver instance.
virtual void start_implementation(system_T& sys,progress_callback*)=0;
};

9
src/sys/tasystem.cpp
View File

@ -129,8 +129,10 @@ vd TaSystem::residual() const {
us total_neqs = arma::sum(neqs);
if(total_neqs>constants::maxndofs) {
throw TaSMETError("Too many DOFS required."
" Problem too large.");
stringstream error;
error << "Too many DOFS required. Problem too large. Number of equations computed: ";
error << total_neqs;
throw TaSMETError(error);
}
// This vector of indices stores the last equation number + 1 for
@ -237,11 +239,12 @@ vus TaSystem::getNDofs() const {
return Ndofs;
}
vus TaSystem::getNEqs() const {
TRACE(0,"TaSystem::getNDofs()");
TRACE(15,"TaSystem::getNEqs()");
vus Neqs(_segs.size());
us i=0;
for (auto seg :_segs) {
Neqs(i)=seg.second->getNEqs(*this);
VARTRACE(15,Neqs(i));
i++;
}
return Neqs;

89
src/tasmet_evalscript.cpp Normal file
View File

@ -0,0 +1,89 @@
// tasmet_evalscript.cpp
//
// last-edit-by: J.A. de Jong
//
// Description:
//
//////////////////////////////////////////////////////////////////////
#include "tasmet_evalscript.h"
#include <chaiscript/chaiscript.hpp>
#include "chaiscript.h"
#include "tasmet_exception.h"
#include "tasmet_tracer.h"
#include <sstream>
using chaiscript::ChaiScript;
using std::stringstream;
template<typename T>
inline T wrap_eval(ChaiScript* chai,const string& script) {
try {
return chai->eval<T>(script);
}
catch(chaiscript::exception::eval_error &e) {
TRACE(15,script);
throw TaSMETError(e.what());
}
catch(std::bad_cast &e) {
TRACE(15,script);
throw TaSMETError("Cannot get return value from script");
}
}
template<>
inline void wrap_eval<void>(ChaiScript* chai,const string& script) {
try {
chai->eval(script);
}
catch(std::runtime_error &e) {
TRACE(15,script);
throw TaSMETError(e.what());
}
catch(std::bad_cast &e) {
TRACE(15,script);
throw TaSMETError("Cannot get return value from script");
}
}
EvaluateFun::EvaluateFun(const string& fun_return,
const string& err_msg):
_err_msg(err_msg),
_fun_return(fun_return)
{
TRACE(15,"EvaluateFun::EvaluateFun()");
_chai = getChaiScriptInstance();
if(!_chai) throw TaSMETBadAlloc();
string script = "def myfun(x) {\n";
script += "return " + fun_return + "; }\n";
wrap_eval<void>(_chai.get(),script);
}
void EvaluateFun::addGlobalDef(const string& name,const d value) {
TRACE(15,"EvaluateFun::addGlobalDef()");
_chai->add_global(chaiscript::var(value),name);
}
vd EvaluateFun::operator()(const vd& x) {
TRACE(15,"EvaluateFun::operator(vd)");
vd y(x.size());
for(us i=0;i<x.size();i++) {
stringstream value;
value << "myfun(" << std::scientific << x(i) << ");";
y(i) = wrap_eval<d>(_chai.get(),value.str());
}
return y;
}
EvaluateFun::~EvaluateFun() {
}
//////////////////////////////////////////////////////////////////////

41
src/tasmet_evalscript.h Normal file
View File

@ -0,0 +1,41 @@
// tasmet_evalscript.h
//
// Author: J.A. de Jong
//
// Description:
//
// Description: This file implements a helper class to evaluate simple
// 1D single-parameter python math functions
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef TASMET_EVALSCRIPT_H
#define TASMET_EVALSCRIPT_H
#include "tasmet_types.h"
#include <memory>
namespace chaiscript {
class ChaiScript;
}
// Helper class to evaluate simple 1D math function evaluations in Python
class EvaluateFun {
std::unique_ptr<chaiscript::ChaiScript> _chai;
string _err_msg;
string _fun_return;
public:
EvaluateFun(const string& fun_return,const string& err_msg = "Script error");
// Add a global definition to the namespace
void addGlobalDef(const string& name,
const d value);
// Evaluate a single function at multiple points and return
// the result for each point
vd operator()(const vd& points);
// Used to cleanup the python namespace
~EvaluateFun();
};
#endif // TASMET_EVALSCRIPT_H
//////////////////////////////////////////////////////////////////////

2
src/tasmet_io.h
View File

@ -10,6 +10,7 @@
#define TASMET_IO_H
#include <iostream>
#include <iomanip>
#include <sstream>
//Spoiling global namespace with often used functions and variables
using std::cout; /* Output to stdout */
@ -20,6 +21,7 @@ using std::cin; /* Input from stdin */
using std::ostream;
using std::cerr;
using std::ios;
using std::stringstream;
#endif // TASMET_IO_H
//////////////////////////////////////////////////////////////////////

26
src/tasmet_solvemodel.cpp
View File

@ -9,9 +9,10 @@
#include "tasmet_io.h"
#include "tasmet_exception.h"
#include "sys/tasystem.h"
#include "solver/newton_raphson.h"
#include "protobuf/model.pb.h"
// For using python from within Qt
#include <PythonQt.h>
#include <memory>
void usage(const char* fn) {
cout << "Usage: " << endl;
@ -22,18 +23,6 @@ int main(int argc, char *argv[]) {
INITTRACE(10);
// Initialize PythonQt
// PythonQt::init(PythonQt::IgnoreSiteModule | PythonQt::RedirectStdOut);
PythonQt::init();
PythonQt* pyqt = PythonQt::self();
PythonQtObjectPtr context = pyqt->getMainModule();
QVariant rv = context.evalScript("from math import *\n");
if(pyqt->hadError()) {
return -1;
}
if(argc != 2) {
usage(argv[0]);
return -1;
@ -53,10 +42,10 @@ int main(int argc, char *argv[]) {
cout << "Loaded model: " << endl;
cout << model.system().DebugString();
TaSystem* system;
std::unique_ptr<TaSystem> system;
try {
system = new TaSystem(model.system());
system = std::unique_ptr<TaSystem>(new TaSystem(model.system()));
}
catch(TaSMETError &e) {
cerr << "Model initialization error: " << endl;
@ -64,6 +53,11 @@ int main(int argc, char *argv[]) {
return -1;
}
NewtonRaphson solver(*system.get());
progress_callback cb = simple_progress_callback(1e-6,1e-6);
solver.start(&cb);
}