Added solvers for one-dimensional root finding

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}  
+  ${PYTHON_INCLUDE_DIRS}
-  common/src/swig
+  src
-  common/src/gas
+  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

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

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)
+add_library(tasmet_src_src
-target_link_libraries(tasmet_src
+  tasmet_tracer.cpp
-  # duct
+  tasmet_exception.cpp
-  # mech
+  tasmet_assert.cpp
-  # seg
-  # sol
-  # sys
-  # var
   )
+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)

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);
+    us nfac = _kappac.size();
-}
+    vd powfac(nfac);
-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];
+    for(us i=0;i<nfac;i++)
+        powfac(i) = pow(T,i);
+        
+    return dot(_kappac,powfac);
 }
 d Air::mu(d T,d 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 = _muc.size();
+    vd powfac(nfac);
+        
+    for(us i=0;i<nfac;i++)
+        powfac(i) = pow(T,i);
+        
+    return dot(_muc,powfac);
 }
 
-  
-
 
 //////////////////////////////////////////////////////////////////////

src/material/air.h

@@ -12,18 +12,26 @@
 
 
 class Air : public PerfectGas {
+    vdi<5> _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(){}
 };
 

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<T.size();i++) \
+        varname_ = varname(T(i),p(i));\
+    return varname_
 
 class Gas{
 public:
@@ -50,7 +55,9 @@ public:
 
     // Specific heat ratio
     d gamma(d T,d p) const {return cp(T,p)/cv(T,p);}
-    vd gamma(const vd& T,const vd& p) const {return cp(T,p)/cv(T,p);}
+    vd gamma(const vd& T,const vd& p) const {
+        element_wise(gamma);
+    }
 
     // Prandtl number
     vd pr(const vd& T,const vd& p) const {return mu(T,p)%cp(T,p)/kappa(T,p);}
@@ -59,49 +66,55 @@ public:
     // Virtuals that are part of the interface
     // Density [kg/m^3]
     virtual d rho(d T,d p) const=0;
-    virtual vd rho(const vd& T,const vd& p) const=0;
+    vd rho(const vd& T,const vd& p) const {
+        element_wise(rho);
+    }
 
     // Adiabatic speed of sound
     virtual d cm(d T,d p) const=0;
-    virtual vd cm(const vd& T,const vd& p) const=0;    
+    vd cm(const vd& T,const vd& p) const {
+        element_wise(cm);
+    }
 
     // Internal energy [J/kg]
     virtual d e(d T,d p) const=0;
-    virtual vd e(const vd& T,const vd& p) const=0;
+    vd e(const vd& T,const vd& p) const {
+        element_wise(e);
+    }
 
     // Static enthalpy [J/kg]
     virtual d h(d T,d p) const=0;
-    virtual vd h(const vd& T,const vd& p) const=0;
+    vd h(const vd& T,const vd& p) const {
+        element_wise(h);
+    }
 
     // Specific heat at constant pressure
     virtual d cp(d T,d p) const=0;
-    virtual vd cp(const vd& T,const vd& p) const=0;
+    vd cp(const vd& T,const vd& p) const {
+        element_wise(cp);
+    }
 
     // Specific heat at constant density
     virtual d cv(d T,d p) const=0;
-    virtual vd cv(const vd& T,const vd& p) const=0;
+    vd cv(const vd& T,const vd& p) const {
-
+        element_wise(cv);
-    virtual d beta(d T,d p) const=0;
+    }
-    virtual vd beta(const vd& T,const vd& p) const=0;
 
     // Dynamic viscosity [Pa*s]
     virtual d mu(d T,d p) const=0;
-    virtual vd mu(const vd& T,const vd& p) const=0;
+    vd mu(const vd& T,const vd& p) const {
+        element_wise(mu);
+    }
 
     // Thermal conductivity [W/mK]
     virtual d kappa(d T,d p) const=0;
-    virtual vd kappa(const vd& T,const vd& p) const=0;
+    vd kappa(const vd& T,const vd& p) const {
-
+        element_wise(kappa);
-
+    }
-
-
-
-
-
-
 
 };
 
+#undef element_wise
 
 
 #endif // GAS_H

src/material/helium.cpp

@@ -9,25 +9,12 @@
 #include "helium.h"
 
   
-vd Helium::cp(const vd& T,const vd& p) const {
-    return cp(0.0,0.0)*arma::ones(T.size());
-}
-vd Helium::h(const vd& T,const vd& p) const {
-    return cp(T,p)%T;
-}
 d Helium::mu(d T,d p) const {
     return 0.412e-6*pow(T,0.68014);
 }
 d Helium::kappa(d T,d p) const {
     return 0.0025672*pow(T,0.716);
 }
-vd Helium::mu(const vd& T,const vd& p) const {
-    return 0.412e-6*pow(T,0.68014);
-}
-vd Helium::kappa(const vd& T,const vd& p) const {
-    return 0.0025672*pow(T,0.716);
-}
-
   
 
 //////////////////////////////////////////////////////////////////////

src/material/helium.h

@@ -11,22 +11,18 @@
 #include "perfectgas.h"
 
 class Helium :public PerfectGas {
+    vdi<1> _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;}
+    const vd& cpc() const { return _cpc; }
-    vd cp(const vd& T,const vd& p) const;
+    d cp(d T,d p) const { return _cpc(0);}
-
-    d h(d T,d p) const {return cp(0.0,0.0)*T;}
-    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;
 
     virtual ~Helium(){}
 };

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];
+    return _kappac(1)*T+_kappac(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);
 }
+
 d Nitrogen::mu(d T,d p) const {
     return (0.01781/1000)*(411.55/(T+111))*pow(T/300.55,1.5);
 }
-
   
 
 //////////////////////////////////////////////////////////////////////

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;
+    const vd& cpc() const { return _cpc; }
-    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;
+
-    ~Nitrogen(){}
 };
 
   

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 {
+    d h(d T,d p) const {
-        checkzero(T);
+        const vd& cpc_ = this->cpc();
-        return p/Rs()/T;
+        us nfac = cpc_.size();
+        vd powfac(nfac);
+        
+        for(us i=0;i<nfac;i++)
+            powfac(i) = pow(T,i)/(1+i);
+        
+        return dot(cpc_,powfac);
     }
-    vd p(const vd& T,const vd& rho) const {
+    d cp(d T,d p) const {
-        return rho%(Rs()*T);
+        const vd& cpc_ = this->cpc();
-    }
+        us nfac = cpc_.size();
-    vd cv(const vd& T,const vd& p) const {
+        vd powfac(nfac);
-        return cp(T,p)-Rs();
+        
-    }
+        for(us i=0;i<nfac;i++)
-    vd e(const vd& T,const vd& p) const {
+            powfac(i) = pow(T,i);
-        return h(T,p)-Rs()*T;
+        return dot(cpc_,powfac);
-    }
-    d beta(d T,d p) const {
-        checkzero(T);
-        return 1/T;
-    }
-    vd beta(const vd& T,const vd& p) const {
-        return 1/T;
-    }
-    vd cm(const vd& T,const vd& p) const {
-        return sqrt(gamma(T,p)*Rs()%T);
     }
     d rho(d T,d p) const {
         checkzero(T);
@@ -73,6 +76,7 @@ public:
     
 };
 
+
 #endif // PERFECTGAS_H
 //////////////////////////////////////////////////////////////////////
 

src/solver/CMakeLists.txt

@@ -0,0 +1,7 @@
+include_directories(${PROJECT_SOURCE_DIR}/src/solver)
+add_library(solver
+  solver.cpp
+  system.cpp
+  newton_raphson.cpp
+  brent.cpp
+  )

src/solver/bracket_root.h

@@ -0,0 +1,91 @@
+// bracket_lin.h
+//
+// Author: J.A. de Jong 
+//
+// Description: Bracket a function's root given a guess value, returns
+// another guess where the solution has the opposite sign
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef BRACKET_LIN_H
+#define BRACKET_LIN_H
+#include "system.h"
+#include "tasmet_tracer.h"
+#include <limits>
+
+inline int sign(d& x){
+    return x < 0? -1: 1;
+}
+const us maxiter = 100;
+
+inline std::pair<d,d> bracket_root(NoGradientNonlinearSystem<d>& 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<d>::epsilon();
+
+    d factor = 1.001;
+    d delta = xa/10 + std::numeric_limits<d>::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<T>::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
+//////////////////////////////////////////////////////////////////////

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 <limits>
+#include "tasmet_constants.h"
+#include <algorithm>            // std::swap
+#include "bracket_root.h"
+
+const d eps = std::numeric_limits<d>::epsilon();
+
+Brent::Brent(const NoGradientNonlinearSystem<d>& 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<d>& system,
+                                 progress_callback* callback) {
+
+    TRACE(15,"Brent::start_implementation");
+
+    us niter = 0;
+
+    d a = system.getSolution();
+
+    std::pair<d,d> 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)<abs(fb)) {
+            std::swap(a,b);
+            std::swap(fa,fb);            
+        }
+        
+        progress.fun_err = abs(fs);
+        progress.rel_err = abs(b-a);
+        progress.iteration++;
+
+        VARTRACE(15,s);
+        VARTRACE(15,a);
+        VARTRACE(15,b); 
+        VARTRACE(15,c);
+        VARTRACE(15,fa);
+        VARTRACE(15,fb); 
+        VARTRACE(15,fc);
+        VARTRACE(15,fs);
+
+        SolverAction action = (*callback)(progress);
+
+        if(action==Stop){
+            break;
+        }
+
+
+    }
+    
+    system.updateSolution(b);
+}
+
+
+//////////////////////////////////////////////////////////////////////

src/solver/brent.h

@@ -0,0 +1,24 @@
+// brent.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Brent's root finding algorithm
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef BRENT_H
+#define BRENT_H
+#include "solver.h"
+
+class Brent: public Solver<NoGradientNonlinearSystem<d>,d> {
+    d _reltol;
+    us _maxiter;
+public:    
+    Brent(const NoGradientNonlinearSystem<d>& sys,us maxiter=10000,d reltol=1e-6);
+protected:
+    void start_implementation(NoGradientNonlinearSystem<d>& sys,progress_callback*);
+};
+
+
+#endif // BRENT_H
+//////////////////////////////////////////////////////////////////////

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;
+        }
+            
+    }
+
+}
+//////////////////////////////////////////////////////////////////////

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<GradientNonlinearSystem,vd> {
+
+    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
+//////////////////////////////////////////////////////////////////////

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 <chrono>
+
+template<typename system_T,typename result_T>
+static void SolverThread(Solver<system_T,result_T>* solver,
+                         system_T* system,
+                         progress_callback* callback);
+
+
+template<typename system_T,typename result_T>
+Solver<system_T,result_T>::Solver(const system_T& sys){
+    _sys = sys.copy();
+    if(!_sys)
+        throw TasMETBadAlloc();
+    _running = false;
+}
+
+template<typename system_T,typename result_T>
+Solver<system_T,result_T>::~Solver(){
+
+    stop();
+    assert(!_running);
+    assert(!_solver_thread);
+    delete _sys;
+    
+}
+
+template<typename system_T,typename result_T>
+void Solver<system_T,result_T>::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<system_T,result_T>,
+                                           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<typename system_T,typename result_T>
+void Solver<system_T,result_T>::stop() {
+
+    _running = false;
+    if(_solver_thread){
+
+        _solver_thread->join();
+        delete _solver_thread;
+        _solver_thread = nullptr;
+
+    }
+}
+template<typename system_T,typename result_T>
+result_T Solver<system_T,result_T>::getSolution() const{
+    if(_running){
+        throw TaSMETError("Solver is running");
+    }
+    // Cleanup thread resources
+    stop();
+    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 a
+template class Solver<NoGradientNonlinearSystem<vd>,vd>;
+template class Solver<GradientNonlinearSystem,vd>;
+template class Solver<NoGradientNonlinearSystem<d>,d>;
+//////////////////////////////////////////////////////////////////////

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 <functional>
+#include <thread>
+#include <atomic>
+
+#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<SolverAction(SolverProgress)> progress_callback;
+
+template<typename system_T,typename result_T>
+class Solver {
+
+protected:
+    system_T* _sys = nullptr;
+    std::thread* _solver_thread = nullptr;
+    std::atomic<bool> _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<typename Y,typename rT>
+    friend void SolverThread(Solver<Y,rT>*,Y*,progress_callback*);
+protected:
+    virtual void start_implementation(system_T& sys,progress_callback*)=0;
+};
+
+
+#endif /* _SOLVER_H_ */
+
+//////////////////////////////////////////////////////////////////////

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<typename T>
+SingleDofNoGradient<T>::SingleDofNoGradient(const NoGradientNonlinearSystem<T>& wrapped_system) {
+
+    this->wrapped_system = wrapped_system.copy();
+}
+template<typename T>
+SingleDofNoGradient<T>::~SingleDofNoGradient(){
+    delete this->wrapped_system;
+}
+template<typename T>
+SingleDofNoGradient<T>* SingleDofNoGradient<T>::copy() const {
+    return new SingleDofNoGradient<T>(*wrapped_system);
+}
+
+template<>
+vd SingleDofNoGradient<d>::residual() const {
+    d val = wrapped_system->residual();
+    return vd({val});
+}
+template<>
+vd SingleDofNoGradient<c>::residual() const {
+    c val = wrapped_system->residual();
+    return vd({val.real(),val.imag()});
+}
+template<>
+vd SingleDofNoGradient<d>::getSolution() const {
+    d val = wrapped_system->getSolution();
+    return vd({val});
+}
+template<>
+vd SingleDofNoGradient<c>::getSolution() const {
+    c val = wrapped_system->getSolution();
+    return vd({val.real(),val.imag()});
+}
+
+template<>
+void SingleDofNoGradient<d>::updateSolution(const vd& new_guess){
+    assert(new_guess.size()==1);
+    wrapped_system->updateSolution(new_guess(0));
+}
+template<>
+void SingleDofNoGradient<c>::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<d>;
+template class SingleDofNoGradient<c>;
+//////////////////////////////////////////////////////////////////////

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<typename T>
+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 <typename T>
+class SingleDofNoGradient: public NoGradientNonlinearSystem<vd> {
+
+    static_assert(is_same<T,d>::value || is_same<T,c>::value,"Instaniation must be complex or double");
+
+    NoGradientNonlinearSystem<T>* wrapped_system = nullptr;
+    
+public:
+    
+    SingleDofNoGradient(const NoGradientNonlinearSystem<T>& wrapped_system);
+    
+    virtual vd residual() const;
+
+    virtual vd getSolution() const;
+
+    // Create a copy of the system
+    virtual SingleDofNoGradient<T>* copy() const;
+
+    virtual void updateSolution(const vd& new_guess);
+
+    ~SingleDofNoGradient();
+
+};
+
+class GradientNonlinearSystem: public NoGradientNonlinearSystem<vd> {
+
+public:
+    
+    virtual sdmat jacobian() const=0;
+
+    GradientNonlinearSystem* copy() const=0;
+};
+
+
+
+#endif // SYSTEM_H
+//////////////////////////////////////////////////////////////////////

src/tasmet_assert.h

@@ -8,12 +8,21 @@
 #pragma once
 #ifndef TASMET_ASSERT_H
 #define TASMET_ASSERT_H
-
-#ifndef TASMET_DEBUG
 #include <type_traits>
+#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<class T, class U>
+struct is_same : std::false_type {};
+ 
+template<class T>
+struct is_same<T, T> : std::true_type {};
+
+
 #if TASMET_DEBUG == 1
 #include <cassert>
 #include "tasmet_io.h"

src/tasmet_tracer.cpp

@@ -8,6 +8,6 @@
 
 #include "tasmet_tracer.h"
 
-
+int tasmet_tracer_level = MAXTRACELEVEL;
 
 //////////////////////////////////////////////////////////////////////

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<size_t i>
+using vdi = arma::Col<d>::fixed<i>;
+template<size_t i>
+using vci = arma::Col<c>::fixed<i>;
 
 typedef arma::Col<d> vd;		  /* Column vector of doubles */		
 typedef arma::Col<c> vc;	  /* Column vector of complex numbers */	

testing/CMakeLists.txt

@@ -0,0 +1 @@
+add_subdirectory(solver)

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)

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> {
+    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<SolverAction(SolverProgress)> 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);
+
+
+}
+
+
+
+//////////////////////////////////////////////////////////////////////

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<SolverAction(SolverProgress)> p = solver_callback;
+    solver.start(&p);
+    solver.stop();
+    
+    cout << "Final solution: x = " << solver.getSolution();
+
+}
+
+
+
+//////////////////////////////////////////////////////////////////////