diff --git a/CMakeLists.txt b/CMakeLists.txt index fddf7b1..a6f0183 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -69,7 +69,7 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -march=native -mtune=native") # For importing find directives for Cmake -set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/common/cmake_tools) +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake_tools) # ########################## # Python ##################### @@ -91,20 +91,6 @@ MESSAGE("Python major version: ${TaSMET_PY_MAJOR_VERSION}") # Find the site_packages directory of python execute_process(COMMAND python${TaSMET_PY_MAJOR_VERSION} -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())" OUTPUT_VARIABLE PYTHON_SITE_PACKAGES OUTPUT_STRIP_TRAILING_WHITESPACE) -# ################################ -# Initialize swig -# ################################ - -find_package(SWIG REQUIRED) -include(${SWIG_USE_FILE}) - -SET(CMAKE_SWIG_FLAGS -Wall -DSWIG_PYTHON) -if(${TaSMET_PY_MAJOR_VERSION}=="3") - SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} -py${TaSMET_PY_MAJOR_VERSION}) -endif(${TaSMET_PY_MAJOR_VERSION}=="3") - -include_directories(common/src) - # Armadillo find_package(Armadillo REQUIRED) add_definitions(-DARMA_USE_SUPERLU -DARMA_USE_CXX11) @@ -114,26 +100,13 @@ add_definitions(-DARMA_USE_SUPERLU -DARMA_USE_CXX11) # This is the common code (gas and solid libs, etc) # link_directories(common) -aux_source_directory(common/src/gas gas) include_directories( - ${PYTHON_INCLUDE_DIRS} - common/src/swig - common/src/gas + ${PYTHON_INCLUDE_DIRS} + src + src/solver + src/material ) # Add the code subdirectory add_subdirectory(src) +add_subdirectory(testing) - -# set_source_files_properties( ${swig_generated_file_fullname} - # PROPERTIES COMPILE_FLAGS "${SWIG_COMMON_COMPILE_FLAGS} ") - - -# ================================== Installation - -# Install common files -install(FILES ${PROJECT_SOURCE_DIR}/common/__init__.py - DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}/common) -install(FILES ${PROJECT_SOURCE_DIR}/common/common.py - DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}/common) - -# Rest of the files is installed from src/CMakeLists.txt diff --git a/doc/usg.lyx b/doc/usg.lyx index ced7676..6928247 100644 --- a/doc/usg.lyx +++ b/doc/usg.lyx @@ -1881,6 +1881,7 @@ Create a PhaseConstraint object: \end_layout \begin_layout Verbatim + pc=PhaseConstraint(Varnr, freqnr, left) \end_layout @@ -1889,6 +1890,7 @@ Apply this contraint to a segment which can accept them, for example a Tube: \end_layout \begin_layout Verbatim + t1.setPhaseConstraint(pc) \end_layout @@ -3768,10 +3770,54 @@ Momentum equation (prescribes pressure \end_layout \begin_layout Standard -And +And the boundary condition for the temperature is computed assuming adiabatic + compression-expansion. + Currently, this is implemented for thermally perfect gases only: +\end_layout + +\begin_layout Standard \begin_inset Formula \begin{equation} -T_{p}=T_{0}\left(\frac{p_{0}+p_{p}}{p_{0}}\right)^{\frac{\gamma_{0}-1}{\gamma_{0}}} +c_{p}\left(T\right)\mathrm{d}T=\frac{\mathrm{d}p}{\rho} +\end{equation} + +\end_inset + +FIlling in the perfect gas law and a bit of bookkeeping results in +\end_layout + +\begin_layout Standard +\begin_inset Formula +\begin{equation} +\frac{1}{R_{s}}\int\limits _{T_{0}}^{T_{p}}\frac{c_{p}\left(T\right)}{T}\mathrm{d}T-\ln\left(\frac{p_{p}-p_{0}}{p_{0}}\right)=0 +\end{equation} + +\end_inset + +Note that in general, to solve this equation for the temperature requires + a numerical integration, however for the currently implemented gases, +\begin_inset Formula $c_{p}$ +\end_inset + + is a polynomial function of +\begin_inset Formula $T$ +\end_inset + +: +\begin_inset Formula +\begin{equation} +c_{p}(T)=\sum_{i=0}^{N_{c_{p}}}c_{p,i}T^{i} +\end{equation} + +\end_inset + +In that case +\end_layout + +\begin_layout Standard +\begin_inset Formula +\begin{equation} +\frac{1}{R_{s}}\left(c_{p,0}\ln\left(\frac{T_{p}}{T_{0}}\right)+\sum_{i=1}^{N_{c_{p}}}\frac{c_{p,i}\left(T_{p}-T_{0}\right)^{i}}{i}\right)-\ln\left(\frac{p_{p}-p_{0}}{p_{0}}\right)=0.\label{eq:T-p-adiabatic} \end{equation} \end_inset @@ -3779,6 +3825,47 @@ T_{p}=T_{0}\left(\frac{p_{0}+p_{p}}{p_{0}}\right)^{\frac{\gamma_{0}-1}{\gamma_{0 \end_layout +\begin_layout Standard +Equation +\begin_inset CommandInset ref +LatexCommand ref +reference "eq:T-p-adiabatic" + +\end_inset + + is solved using a one-dimensional root finding algorithm (see Section +\begin_inset CommandInset ref +LatexCommand ref +reference "sec:One-dimensional-function-solvers" + +\end_inset + +). + Note that for an ideal gas an explicit formula is available: +\end_layout + +\begin_layout Standard +\begin_inset Formula +\begin{equation} +T_{p,\mathrm{ideal}}=T_{0}\left(\frac{p_{0}+p_{p}}{p_{0}}\right)^{\frac{\gamma_{0}-1}{\gamma_{0}}}.\label{eq:ideal-gas-isentropic-p-T} +\end{equation} + +\end_inset + +Looking closely at Equation +\begin_inset CommandInset ref +LatexCommand ref +reference "eq:ideal-gas-isentropic-p-T" + +\end_inset + +, we find that +\begin_inset Formula $T_{p,\mathrm{ideal}}$ +\end_inset + + provides a good guess for the final solution. +\end_layout + \begin_layout Subsection AdiabaticWall \end_layout @@ -5044,8 +5131,45 @@ Minor loss Systems \end_layout -\begin_layout Section +\begin_layout Chapter +Solvers +\end_layout +\begin_layout Section +One-dimensional function solvers +\begin_inset CommandInset label +LatexCommand label +name "sec:One-dimensional-function-solvers" + +\end_inset + + +\end_layout + +\begin_layout Subsection +Gradient-based +\end_layout + +\begin_layout Subsection +Gradient free +\begin_inset CommandInset label +LatexCommand label +name "subsec:Gradient-free" + +\end_inset + + +\end_layout + +\begin_layout Standard +As a gradient free one-dimensional function solver, we use Brent's mehhod. + Brent's method combines root bracketing, bisection and inverse quadratic + interpolation to find the root of the function without using the gradient. + See Wikipedia for more information. +\end_layout + +\begin_layout Section +Minimizers \end_layout \end_body diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2966575..9cb5666 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,7 +18,8 @@ include_directories( # python # seg # sol - # sys + sys + solver # var # volume ) @@ -28,19 +29,25 @@ add_subdirectory(material) # add_subdirectory(duct) # add_subdirectory(mech) # add_subdirectory(seg) -# add_subdirectory(sol) +add_subdirectory(solver) add_subdirectory(sys) # add_subdirectory(var) -add_library(tasmet_src tasmet_tracer.cpp tasmet_exception.cpp tasmet_assert.cpp) -target_link_libraries(tasmet_src - # duct - # mech - # seg - # sol - # sys - # var +add_library(tasmet_src_src + tasmet_tracer.cpp + tasmet_exception.cpp + tasmet_assert.cpp ) +add_library(tasmet_src) +target_link_libraries(tasmet_src + funcs + material + solver + sys + # This one should be last as other parts link to these utilities + tasmet_src_src + ) + # set_source_files_properties(swig/nonlin.i # PROPERTIES CPLUSPLUS ON) diff --git a/src/material/air.cpp b/src/material/air.cpp index 5a8f321..c4c0f6d 100644 --- a/src/material/air.cpp +++ b/src/material/air.cpp @@ -8,47 +8,24 @@ #include "air.h" -namespace { - // Air-specific data - d kappac[]={-0.00227583562,1.15480022E-4, \ - -7.90252856E-8,4.11702505E-11, \ - -7.43864331E-15}; - - d cpc[]={1047.63657,-0.372589265, \ - 9.45304214E-4,-6.02409443E-7, \ - 1.2858961E-10}; - d muc[]={-8.38278E-7,8.35717342E-8, \ - -7.69429583E-11,4.6437266E-14, \ - -1.06585607E-17}; -} // namespace - -d Air::h(d T,d p) const { - return cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5); -} -vd Air::h(const vd& T,const vd& p) const { - return cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5); -} -vd Air::cp(const vd& T,const vd& p) const { - return cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4); -} -d Air::cp(d T,d p) const { - return cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4); -} - -vd Air::kappa(const vd& T,const vd& p) const { - return kappac[0]+kappac[1]*T+kappac[2]*pow(T,2)+kappac[3]*pow(T,3)+kappac[4]*pow(T,4); -} d Air::kappa(d T,d p) const { - return kappac[0]+kappac[1]*T+kappac[2]*pow(T,2)+kappac[3]*pow(T,3)+kappac[4]*pow(T,4); -} -vd Air::mu(const vd& T,const vd& p) const { - return muc[1]*T+muc[2]*pow(T,2)+muc[3]*pow(T,3)+muc[4]*pow(T,4)+muc[0]; + us nfac = _kappac.size(); + vd powfac(nfac); + + for(us i=0;i _cpc = {1047.63657,-0.372589265, \ + 9.45304214E-4,-6.02409443E-7, \ + 1.2858961E-10}; + + vdi<5> _kappac = {-0.00227583562,1.15480022E-4, \ + -7.90252856E-8,4.11702505E-11, \ + -7.43864331E-15}; + + vdi<5> _muc = {-8.38278E-7,8.35717342E-8, \ + -7.69429583E-11,4.6437266E-14, \ + -1.06585607E-17}; protected: d Rs() const {return 287;} public: Air(d T0,d p0):PerfectGas(air,T0,p0){} + const vd& cpc() const {return _cpc;} d cp(d T,d p) const; - vd cp(const vd& T,const vd& p) const; d h(d T,d p) const; - vd h(const vd& T,const vd& p) const; d mu(d T,d p) const; - vd mu(const vd& T,const vd& p) const; d kappa(d T,d p) const; - vd kappa(const vd& T,const vd& p) const; ~Air(){} }; diff --git a/src/material/gas.h b/src/material/gas.h index 3aba24b..ecb1d14 100644 --- a/src/material/gas.h +++ b/src/material/gas.h @@ -14,6 +14,11 @@ #include "tasmet_types.h" #include "tasmet_constants.h" +#define element_wise(varname) \ + vd varname_(T.size()); \ + for(us i=0;i _cpc = {5195}; protected: d Rs() const {return 2070;} public: Helium(d T0,d p0):PerfectGas(helium,T0,p0){} - d cp(d T,d p) const { return 5195;} - vd cp(const vd& T,const vd& p) const; - - d h(d T,d p) const {return cp(0.0,0.0)*T;} - vd h(const vd& T,const vd& p) const; + const vd& cpc() const { return _cpc; } + d cp(d T,d p) const { return _cpc(0);} d mu(d T,d p) const; - vd mu(const vd& T,const vd& p) const; d kappa(d T,d p) const; - vd kappa(const vd& T,const vd& p) const; virtual ~Helium(){} }; diff --git a/src/material/nitrogen.cpp b/src/material/nitrogen.cpp index 8ee8bb8..fe4d96f 100644 --- a/src/material/nitrogen.cpp +++ b/src/material/nitrogen.cpp @@ -8,42 +8,13 @@ #include "nitrogen.h" -namespace { - // Nitrogen-specific data - d cpc[]={3.29868,0.00140824, \ - -3.96322e-6,5.64152e-9, \ - -2.44486e-12}; - - d kappavals[]={6.995e-3,0.0631e-3}; -} // namespace - -d Nitrogen::h(d T,d p) const { - return Rs()*(cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5)); -} -vd Nitrogen::h(const vd& T,const vd& p) const { - return Rs()*(cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5)); -} -vd Nitrogen::cp(const vd& T,const vd& p) const { - return Rs()*(cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4)); -} -d Nitrogen::cp(d T,d p) const { - return Rs()*(cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4)); -} -// Document Mina -vd Nitrogen::kappa(const vd& T,const vd& p) const { - return kappavals[1]*T+kappavals[0]; -} d Nitrogen::kappa(d T,d p) const { - return kappavals[1]*T+kappavals[0]; -} -// http://www.lmnoeng.com/Flow/GasViscosity.php -vd Nitrogen::mu(const vd& T,const vd& p) const { - return (0.01781/1000)*(411.55/(T+111))%pow(T/300.55,1.5); + return _kappac(1)*T+_kappac(0); } + d Nitrogen::mu(d T,d p) const { return (0.01781/1000)*(411.55/(T+111))*pow(T/300.55,1.5); } - ////////////////////////////////////////////////////////////////////// diff --git a/src/material/nitrogen.h b/src/material/nitrogen.h index 9e5068f..fb7c9be 100644 --- a/src/material/nitrogen.h +++ b/src/material/nitrogen.h @@ -13,20 +13,22 @@ class Nitrogen : public PerfectGas { + vdi<5> _cpc={3.29868,0.00140824, + -3.96322e-6, + 5.64152e-9, + -2.44486e-12}; + + vdi<2> _kappac={6.995e-3,0.0631e-3}; + protected: d Rs() const {return 297;} public: Nitrogen(d T0,d p0):PerfectGas(nitrogen,T0,p0){} - - d cp(d T,d p) const; - vd cp(const vd& T,const vd& p) const; - d h(d T,d p) const; - vd h(const vd& T,const vd& p) const; + + const vd& cpc() const { return _cpc; } d mu(d T,d p) const; - vd mu(const vd& T,const vd& p) const; d kappa(d T,d p) const; - vd kappa(const vd& T,const vd& p) const; - ~Nitrogen(){} + }; diff --git a/src/material/perfectgas.h b/src/material/perfectgas.h index 226be17..38f1caf 100644 --- a/src/material/perfectgas.h +++ b/src/material/perfectgas.h @@ -2,8 +2,11 @@ // // Author: J.A. de Jong // -// Description: -// +// Description: A thermally perfect gas obeys the ideal-gas law +// equation of state (p=rho*R*T) but has TEMPERATURE-DEPENDENT and +// PRESSURE-DEPENDENT specific heat constants. Our definition of a +// perfect gas has only temperature-dependent specific heat, viscosity +// and thermal conductivity values. ////////////////////////////////////////////////////////////////////// #pragma once #ifndef PERFECTGAS_H @@ -28,33 +31,33 @@ protected: // Not implemented for a perfect gas: // mu, kappa, h, cp, Rs +public: virtual d Rs() const=0; -public: + // Vector of polynomial coefficients to compute the + // specific heat at constant pressure + virtual const vd& cpc() const=0; + ~PerfectGas(){} - vd rho(const vd&T,const vd&p) const { - checkzero(T); - return p/Rs()/T; + d h(d T,d p) const { + const vd& cpc_ = this->cpc(); + us nfac = cpc_.size(); + vd powfac(nfac); + + for(us i=0;icpc(); + us nfac = cpc_.size(); + vd powfac(nfac); + + for(us i=0;i + +inline int sign(d& x){ + return x < 0? -1: 1; +} +const us maxiter = 100; + +inline std::pair bracket_root(NoGradientNonlinearSystem& sys,const d guess) { + + TRACE(10,"bracket_root()"); + + d xa = guess; + d fa = sys.residual(xa); + + d fb = fa; + // We add here a small amount, such that, in case the first guess + // is 0, we + d xb = xa+std::numeric_limits::epsilon(); + + d factor = 1.001; + d delta = xa/10 + std::numeric_limits::epsilon(); + + us iter = 0; + + int step = 64; + + bool switched = false; + + while (iter < maxiter) { + + VARTRACE(5,xb); + + xb += delta*factor; + + fb = sys.residual(xb); + + if(sign(fa) != sign(fb)) { + // Found it! + VARTRACE(5,xa); + VARTRACE(5,xb); + VARTRACE(5,fa); + VARTRACE(5,fb); + return std::make_pair(xa,xb); + } + + if (abs(fb) < abs(fa)){ + // We found a point which is closer to the real root + xa = xb; + fa = fb; + } + else if(!switched) { + TRACE(5,"Switching search direction"); + switched = true; + // We searched in the wrong direction + factor*=-1; + } + + // Heuristic: normally it's best not to increase the step sizes as we'll just end up + // with a really wide range to search for the root. However, if the initial guess was *really* + // bad then we need to speed up the search otherwise we'll take forever if we're orders of + // magnitude out. This happens most often if the guess is a small value (say 1) and the result + // we're looking for is close to std::numeric_limits::min(). + // + + if((10000 - iter) % step == 0) { + factor *=2; + if(step > 1) step/=2; + } + + iter++; + } + + WARN("Failed to bracket root"); + + return std::make_pair(xa,xb); + +} +#endif // BRACKET_LIN_H +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/brent.cpp b/src/solver/brent.cpp new file mode 100644 index 0000000..550968a --- /dev/null +++ b/src/solver/brent.cpp @@ -0,0 +1,173 @@ +// brent.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Brent's root finding algorithm as implemented from Wikipedia: +// https://en.wikipedia.org/wiki/Brent's_method +////////////////////////////////////////////////////////////////////// + +#include "brent.h" +#include "tasmet_tracer.h" +#include "tasmet_exception.h" +#include +#include "tasmet_constants.h" +#include // std::swap +#include "bracket_root.h" + +const d eps = std::numeric_limits::epsilon(); + +Brent::Brent(const NoGradientNonlinearSystem& sys,us maxiter,d reltol): + Solver(sys), + _reltol(reltol), + _maxiter(maxiter) +{ + + TRACE(15,"Brent::Brent"); + #ifdef TASMET_DEBUG + bool ok=true; + ok&=(maxiter>0); + ok&=(reltol >= eps); + + if(!ok){ + throw TaSMETError("Brent: invalid arguments"); + } + #endif +} +namespace { + inline bool is_between(d x,d a,d b) { + const d& mi = std::min(a,b); + const d& ma = std::max(a,b); + return ((mi <= x) && (x <=ma)); + } +} + +void Brent::start_implementation(NoGradientNonlinearSystem& system, + progress_callback* callback) { + + TRACE(15,"Brent::start_implementation"); + + us niter = 0; + + d a = system.getSolution(); + + std::pair brackets = bracket_root(system,a); + a = brackets.first; + d b = brackets.second; + + d fa = system.residual(a); + d fb = system.residual(b); + + d c = a; + d fc = fa; + + d s; // Test point + d fs; // Function at test point + + d t,u,v; + d d_; + + if(fa*fb>=0){ + WARN("Brent solver failed: root is not bracketed"); + return; + } + + if(abs(fa) < abs(fb)) { + std::swap(a,b); + std::swap(fa,fb); + } + + bool mflag = true; + bool bisec_flag; + SolverProgress progress; + + while(_running && progress.iteration <=_maxiter) { + + if((fa!=fc) && (fb!=fc)){ + // Inverse quadratic interpolation + TRACE(15,"IQI"); + t = a*fb*fc/((fa-fb)*(fa-fc)); + u = b*fa*fc/((fb-fa)*(fb-fc)); + v = c*fa*fb/((fc-fa)*(fc-fb)); + + s = t+u+v; + } + else { + // Secant method + TRACE(15,"Secant"); + s = b-fb*(b-a)/(fb-fa); + } + + bisec_flag = false; + + d abssmb = abs(s-b); + d absbmc = abs(b-c); + d abscmd = abs(c-d_); + + VARTRACE(15,s); + + if((bisec_flag |= (!is_between(s,(3*a+b)/4,b)))) goto bflag; + TRACE(15,"Survived 1"); + if(bisec_flag |= (mflag && (abssmb >= absbmc/2))) goto bflag; + TRACE(15,"Survived 2"); + if(bisec_flag |= ((!mflag) && (abssmb >= abscmd/2))) goto bflag; + TRACE(15,"Survived 3"); + if(bisec_flag |= (mflag && (absbmc < abs(_reltol)))) goto bflag;; + TRACE(15,"Survived 4"); + bflag: + if(bisec_flag || ((!mflag) && (abscmd < abs(_reltol)))) { + TRACE(15,"Bisection"); + s = (a+b)/2; + mflag = true; + } + else { + mflag = false; + } + + fs = system.residual(s); + d_ = c; + + if(fa*fs < 0) { + c = b; + fc = fb; + b=s; + fb = fs; + } + else { + c = a; + fc = fa; + a = s; + fa = fs; + } + if(abs(fa),d> { + d _reltol; + us _maxiter; +public: + Brent(const NoGradientNonlinearSystem& sys,us maxiter=10000,d reltol=1e-6); +protected: + void start_implementation(NoGradientNonlinearSystem& sys,progress_callback*); +}; + + +#endif // BRENT_H +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/newton_raphson.cpp b/src/solver/newton_raphson.cpp new file mode 100644 index 0000000..975bc13 --- /dev/null +++ b/src/solver/newton_raphson.cpp @@ -0,0 +1,57 @@ +// newton_raphson.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "newton_raphson.h" +#include "tasmet_tracer.h" + +void NewtonRaphson::start_implementation(GradientNonlinearSystem& system, + progress_callback* callback) { + + assert(callback); + + us niter = 0; + + // Relative error, function error + d reler,funer; + + vd guess = system.getSolution(); + + vd residual=system.residual(); + + SolverProgress progress; + SolverAction action; + + while (_running && progress.iteration<=_maxiter) { + + sdmat jac=system.jacobian(); + + assert(jac.n_cols==residual.size()); + assert(jac.n_rows==residual.size()); + + vd dx = -1*_dampfac*arma::spsolve(jac,residual,"superlu"); + + guess += dx; + + system.updateSolution(guess); + + residual = system.residual(); + + progress.rel_err = norm(dx); + progress.fun_err = norm(residual); + progress.iteration++; + + action = (*callback)(progress); + + if(action==Stop){ + break; + } + + } + +} +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/newton_raphson.h b/src/solver/newton_raphson.h new file mode 100644 index 0000000..d55e595 --- /dev/null +++ b/src/solver/newton_raphson.h @@ -0,0 +1,30 @@ +// newton_raphson.h +// +// Author: J.A. de Jong +// +// Description: +// Implementation of the Newton-Raphson method +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef NEWTON_RAPHSON_H +#define NEWTON_RAPHSON_H +#include "solver.h" + +class NewtonRaphson: public Solver { + + d _dampfac = 1.0; + d _maxiter = 10000; +public: + NewtonRaphson(const GradientNonlinearSystem&sys,d dampfac,us maxiter): + Solver(sys), + _dampfac(dampfac), + _maxiter(maxiter){} + +protected: + void start_implementation(GradientNonlinearSystem& system,progress_callback* callback); + +}; + + +#endif // NEWTON_RAPHSON_H +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/solver.cpp b/src/solver/solver.cpp new file mode 100644 index 0000000..267c0ba --- /dev/null +++ b/src/solver/solver.cpp @@ -0,0 +1,101 @@ +// solver.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// +////////////////////////////////////////////////////////////////////// + +#include "solver.h" +#include "tasmet_exception.h" +#include + +template +static void SolverThread(Solver* solver, + system_T* system, + progress_callback* callback); + + +template +Solver::Solver(const system_T& sys){ + _sys = sys.copy(); + if(!_sys) + throw TasMETBadAlloc(); + _running = false; +} + +template +Solver::~Solver(){ + + stop(); + assert(!_running); + assert(!_solver_thread); + delete _sys; + +} + +template +void Solver::start(progress_callback* callback,bool wait){ + + if(_running){ + assert(_solver_thread); + throw TaSMETError("Solver already running"); + } + + assert(_solver_thread == nullptr); + + this->_solver_thread = new std::thread(SolverThread, + this, + _sys, + callback); + if(!_solver_thread) + throw TasMETBadAlloc(); + + if(wait){ + while (_running){ + std::this_thread::sleep_for(std::chrono::milliseconds(1)); + } + // Cleanup resources + stop(); + } +} + +template +void Solver::stop() { + + _running = false; + if(_solver_thread){ + + _solver_thread->join(); + delete _solver_thread; + _solver_thread = nullptr; + + } +} +template +result_T Solver::getSolution() const{ + if(_running){ + throw TaSMETError("Solver is running"); + } + // Cleanup thread resources + stop(); + return _sys->getSolution(); +} +template +void SolverThread(Solver* solver,system_T* system,progress_callback* callback) { + assert(system); + + solver->_running = true; + + solver->start_implementation(*system,callback); + + solver->_running = false; + +} + + +// Explicit instantiation for a +template class Solver,vd>; +template class Solver; +template class Solver,d>; +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/solver.h b/src/solver/solver.h new file mode 100644 index 0000000..25085af --- /dev/null +++ b/src/solver/solver.h @@ -0,0 +1,64 @@ +// solver.h +// +// Author: J.A. de Jong +// +// Description: +// When solve() is called, a copy of itself is created on the heap, +// for which a thread will be run. After the thread is done, the +// updated (solved) TaSystem will overwrite the old one. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef SOLVER_H +#define SOLVER_H +#include +#include +#include + +#include "tasmet_types.h" +#include "tasmet_assert.h" +#include "system.h" + +struct SolverProgress +{ + size_t iteration = 0; + d fun_err; + d rel_err; + bool done = false; +}; + +enum SolverAction{ + Continue=0, + Stop = 1 +}; + +typedef std::function progress_callback; + +template +class Solver { + +protected: + system_T* _sys = nullptr; + std::thread* _solver_thread = nullptr; + std::atomic _running; +public: + Solver(const system_T& sys); + Solver(const Solver&)=delete; + Solver& operator=(const Solver&)=delete; + + void start(progress_callback* callback=nullptr,bool wait=true); + void stop(); // Stops the solver + + // Returns the solution of the problem + result_T getSolution() const; + + virtual ~Solver(); + template + friend void SolverThread(Solver*,Y*,progress_callback*); +protected: + virtual void start_implementation(system_T& sys,progress_callback*)=0; +}; + + +#endif /* _SOLVER_H_ */ + +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/system.cpp b/src/solver/system.cpp new file mode 100644 index 0000000..d79a30b --- /dev/null +++ b/src/solver/system.cpp @@ -0,0 +1,61 @@ +// system.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: Instantiations for the SingleDofNoGradient class for +// the float and complex type +////////////////////////////////////////////////////////////////////// + +#include "system.h" + + +template +SingleDofNoGradient::SingleDofNoGradient(const NoGradientNonlinearSystem& wrapped_system) { + + this->wrapped_system = wrapped_system.copy(); +} +template +SingleDofNoGradient::~SingleDofNoGradient(){ + delete this->wrapped_system; +} +template +SingleDofNoGradient* SingleDofNoGradient::copy() const { + return new SingleDofNoGradient(*wrapped_system); +} + +template<> +vd SingleDofNoGradient::residual() const { + d val = wrapped_system->residual(); + return vd({val}); +} +template<> +vd SingleDofNoGradient::residual() const { + c val = wrapped_system->residual(); + return vd({val.real(),val.imag()}); +} +template<> +vd SingleDofNoGradient::getSolution() const { + d val = wrapped_system->getSolution(); + return vd({val}); +} +template<> +vd SingleDofNoGradient::getSolution() const { + c val = wrapped_system->getSolution(); + return vd({val.real(),val.imag()}); +} + +template<> +void SingleDofNoGradient::updateSolution(const vd& new_guess){ + assert(new_guess.size()==1); + wrapped_system->updateSolution(new_guess(0)); +} +template<> +void SingleDofNoGradient::updateSolution(const vd& new_guess){ + assert(new_guess.size()==2); + wrapped_system->updateSolution(c(new_guess(0),new_guess(1))); +} + +// Explicit intantiation for float-type and complex-type +template class SingleDofNoGradient; +template class SingleDofNoGradient; +////////////////////////////////////////////////////////////////////// diff --git a/src/solver/system.h b/src/solver/system.h new file mode 100644 index 0000000..3a0b555 --- /dev/null +++ b/src/solver/system.h @@ -0,0 +1,74 @@ +// system.h +// +// Author: J.A. de Jong +// +// Description: +// This header file describes an interface to a system of equations +// that can be solved using a (non)linear solver. +////////////////////////////////////////////////////////////////////// +#pragma once +#ifndef SYSTEM_H +#define SYSTEM_H +#include "tasmet_types.h" +#include "tasmet_assert.h" + +template +class NoGradientNonlinearSystem{ + // A nonlinear system for which the gradient of the residual to + // the solution vector cannot be obtained. +public: + virtual T residual() const=0; + + virtual T residual(const T& guess) { updateSolution(guess); return residual();} + + // Obtain an initial guess of the solution + virtual T getSolution() const=0; + + // Create a copy of the system + virtual NoGradientNonlinearSystem* copy() const=0; + + virtual void updateSolution(const T& new_guess)=0; + + virtual ~NoGradientNonlinearSystem(){} +}; + + +// Wrapper around a vector-valued system, such that a solver for a +// multi-dimensional system can be called for a single-DOF system +template +class SingleDofNoGradient: public NoGradientNonlinearSystem { + + static_assert(is_same::value || is_same::value,"Instaniation must be complex or double"); + + NoGradientNonlinearSystem* wrapped_system = nullptr; + +public: + + SingleDofNoGradient(const NoGradientNonlinearSystem& wrapped_system); + + virtual vd residual() const; + + virtual vd getSolution() const; + + // Create a copy of the system + virtual SingleDofNoGradient* copy() const; + + virtual void updateSolution(const vd& new_guess); + + ~SingleDofNoGradient(); + +}; + +class GradientNonlinearSystem: public NoGradientNonlinearSystem { + +public: + + virtual sdmat jacobian() const=0; + + GradientNonlinearSystem* copy() const=0; +}; + + + +#endif // SYSTEM_H +////////////////////////////////////////////////////////////////////// diff --git a/src/tasmet_assert.h b/src/tasmet_assert.h index 3b208c9..e9f4a0e 100644 --- a/src/tasmet_assert.h +++ b/src/tasmet_assert.h @@ -8,12 +8,21 @@ #pragma once #ifndef TASMET_ASSERT_H #define TASMET_ASSERT_H - -#ifndef TASMET_DEBUG #include +#ifndef TASMET_DEBUG + static_assert(false, "TASMET_DEBUG macro not defined. Please set it to 1 or 0"); + #endif +// Compile-time type checking for template instantiation. +template +struct is_same : std::false_type {}; + +template +struct is_same : std::true_type {}; + + #if TASMET_DEBUG == 1 #include #include "tasmet_io.h" diff --git a/src/tasmet_tracer.cpp b/src/tasmet_tracer.cpp index d7f1d5f..cf5164d 100644 --- a/src/tasmet_tracer.cpp +++ b/src/tasmet_tracer.cpp @@ -8,6 +8,6 @@ #include "tasmet_tracer.h" - +int tasmet_tracer_level = MAXTRACELEVEL; ////////////////////////////////////////////////////////////////////// diff --git a/src/tasmet_types.h b/src/tasmet_types.h index 94e1729..aa276dd 100644 --- a/src/tasmet_types.h +++ b/src/tasmet_types.h @@ -29,7 +29,7 @@ typedef size_t us; /* Size type I always use */ // To change the whole code to 32-bit floating points, change this to // float. #if TASMET_FLOAT == 32 -// typedef float d; /* Shortcut for double */ +typedef float d; /* Shortcut for double */ #elif TASMET_FLOAT == 64 typedef double d; /* Shortcut for double */ #else @@ -67,6 +67,10 @@ const c sqmI=sqrt(-1.0*I); const c minI(0,-1); //Minus complex unity const d number_pi=arma::datum::pi; // The number 3.14159265359.. +template +using vdi = arma::Col::fixed; +template +using vci = arma::Col::fixed; typedef arma::Col vd; /* Column vector of doubles */ typedef arma::Col vc; /* Column vector of complex numbers */ diff --git a/testing/CMakeLists.txt b/testing/CMakeLists.txt new file mode 100644 index 0000000..7b53451 --- /dev/null +++ b/testing/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(solver) diff --git a/testing/solver/CMakeLists.txt b/testing/solver/CMakeLists.txt new file mode 100644 index 0000000..e81a9f9 --- /dev/null +++ b/testing/solver/CMakeLists.txt @@ -0,0 +1,4 @@ +add_executable(newton_raphson_test newton_raphson_test.cpp) +add_executable(brent_test brent_test.cpp) +target_link_libraries(brent_test tasmet_src armadillo openblas pthread) +target_link_libraries(newton_raphson_test tasmet_src armadillo openblas pthread) diff --git a/testing/solver/brent_test.cpp b/testing/solver/brent_test.cpp new file mode 100644 index 0000000..2deaab7 --- /dev/null +++ b/testing/solver/brent_test.cpp @@ -0,0 +1,69 @@ +// newton_raphson_test.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Test program to test the Newton-Raphson solver implementation +////////////////////////////////////////////////////////////////////// +#include "system.h" +#include "brent.h" +#include "tasmet_io.h" +#include "tasmet_tracer.h" +#include "helium.h" +#include "nitrogen.h" +#include "air.h" + +SolverAction solver_callback(SolverProgress p){ + + cout << "Iteration: " << p.iteration << ". Function error: " << p.fun_err <<". Relative error: " << p.rel_err << endl; + if(p.fun_err < 1e-7){ + TRACE(9, "Returning stop"); + return Stop; + } + else { + TRACE(9, "Returning continue"); + return Continue; + } +} + +class TestFunction: public NoGradientNonlinearSystem { + d cur_sol; +public: + TestFunction(d guess){cur_sol = guess;} + d residual() const { + return pow(cur_sol,2)-3; + } + + void updateSolution(const d& res){ cur_sol = res;} + d getSolution() const { return cur_sol;} + TestFunction* copy() const {return new TestFunction(cur_sol);} + sdmat jacobian() const { + sdmat res(1,1); + res(0,0) = 2*cur_sol; + return res; + } +}; + +int main(){ + + INITTRACE(5); + + TestFunction t(100); + Brent solver(t); + + std::function p = solver_callback; + solver.start(&p); + + d res = solver.getSolution(); + + cout << "Final solution: x = " << res << endl; + + Nitrogen nit(293.15,101325); + Air air(293.15,101325); + + +} + + + +////////////////////////////////////////////////////////////////////// diff --git a/testing/solver/newton_raphson_test.cpp b/testing/solver/newton_raphson_test.cpp new file mode 100644 index 0000000..8e62551 --- /dev/null +++ b/testing/solver/newton_raphson_test.cpp @@ -0,0 +1,62 @@ +// newton_raphson_test.cpp +// +// last-edit-by: J.A. de Jong +// +// Description: +// Test program to test the Newton-Raphson solver implementation +////////////////////////////////////////////////////////////////////// +#include "system.h" +#include "newton_raphson.h" +#include "tasmet_io.h" +#include "tasmet_tracer.h" + +SolverAction solver_callback(SolverProgress p){ + + cout << "Iteration: " << p.iteration << ". Function error: " << p.fun_err << endl; + if(p.fun_err < 1e-7){ + TRACE(9, "Returning stop"); + return Stop; + } + else { + TRACE(9, "Returning continue"); + return Continue; + } +} + +class TestFunction: public GradientNonlinearSystem { + d cur_sol; +public: + TestFunction(d guess){cur_sol = guess;} + vd residual() const { + d val = pow(cur_sol,2)-3; + return vd({val}); + } + + void updateSolution(const vd& res){ cur_sol = res(0);} + vd getSolution() const { return vd({cur_sol});} + TestFunction* copy() const {return new TestFunction(cur_sol);} + sdmat jacobian() const { + sdmat res(1,1); + res(0,0) = 2*cur_sol; + return res; + } +}; + +int main(){ + + INITTRACE(10); + + TestFunction t(100); + NewtonRaphson solver(t,1.0,1000); + + std::function p = solver_callback; + solver.start(&p); + solver.stop(); + + cout << "Final solution: x = " << solver.getSolution(); + +} + + + +//////////////////////////////////////////////////////////////////////