Commit 86b33b9a authored by Malte Heithoff's avatar Malte Heithoff
Browse files

daecpp to native

parent acd214b5
FILE(GLOB DAE_CPP_SOURCES *.cpp)
FILE(GLOB DAE_CPP_INCLUDES *.h)
FILE(GLOB DAE_CPP_INCLUDES_PLOTTING external/matplotlib-cpp/*.h)
if(UNIX AND NOT APPLE)
add_library(daecpp SHARED ${DAE_CPP_SOURCES})
install(TARGETS daecpp DESTINATION lib)
endif(UNIX AND NOT APPLE)
add_library(daecpp_static STATIC ${DAE_CPP_SOURCES})
install(TARGETS daecpp_static DESTINATION lib)
install(FILES ${DAE_CPP_INCLUDES} DESTINATION include)
install(FILES ${DAE_CPP_INCLUDES_PLOTTING} DESTINATION include/external/matplotlib-cpp)
#cmakedefine DAE_FORTRAN_STYLE
#cmakedefine DAE_SINGLE
/*
* A set of helper functions to print on screen or write to files
* Jacobian matrix, Mass matrix, the RHS for debugging purposes.
*/
#include <iostream> // std::cout
#include <iomanip> // std::setw etc.
#include <fstream> // File output
#include <string> // std::string, std::to_string
#include <cmath> // std::abs
#include "RHS.h"
#include "mass_matrix.h"
#include "jacobian.h"
namespace daecpp_namespace_name
{
const char delimiter = '\t'; // Delimiter of columns in output text files
/*
* Helper function to write the RHS vector to a file
*/
void RHS::dump(const state_type &x, const double t)
{
std::cout << "RHS::dump() -- INFO: Writing the RHS at time t = "
<< t << "...\n";
const MKL_INT size = x.size();
state_type f(size); // the vector to be saved
this->operator()(x, f, t); // calls the RHS
std::ofstream outFile;
outFile.open("dump_RHS_" + std::to_string(m_dump_file_counter++) + ".txt");
outFile << "t = " << t << ":\n";
outFile << "i" << delimiter << "x[i]" << delimiter << "RHS[i]" << '\n';
for(MKL_INT i = 0; i < size; i++)
outFile << i << delimiter << x[i] << delimiter << f[i] << '\n';
outFile.close();
}
/*
* Helper function to write the Mass matrix to a file
*/
void MassMatrix::dump()
{
std::cout << "MassMatrix::dump() -- INFO: Writing the Mass matrix...\n";
sparse_matrix_holder M;
this->operator()(M); // calls the Mass matrix operator
m_matrix_converter(M); // converts the matrix if it is in simple form
const MKL_INT size =
M.ia.size() - 1; // derive the matrix size from ia index
if(size > 10000)
{
std::cout << "MassMatrix::dump() -- WARNING: the size of the Mass "
"matrix for writting is bigger than 10000x10000.\n";
}
std::ofstream outFile;
MKL_INT ja = 0;
outFile.open("dump_Mass_matrix.txt"); // Mass matrix is static - one file
outFile << "i,j:";
for(MKL_INT i = 0; i < size; i++)
{
outFile << delimiter << "i=" << i;
}
outFile << '\n';
for(MKL_INT j = 0; j < size; j++)
{
MKL_INT ent = M.ia[j + 1] - M.ia[j]; // Number of entries in line j
outFile << "j=" << j << delimiter;
for(MKL_INT i = 0; i < size; i++)
{
if(M.ja[ja] == i)
{
outFile << M.A[ja++];
if(!(--ent))
break;
}
outFile << delimiter;
}
outFile << '\n';
}
outFile.close();
}
/*
* Helper function to write the Jacbian matrix to a file (in dense format)
*/
void Jacobian::dump(const state_type &x, const double t)
{
std::cout << "Jacobian::dump() -- INFO: ";
sparse_matrix_holder M;
this->operator()(M, x, t); // calls the Jacobian matrix operator
m_matrix_converter(M); // converts the matrix if it is in simple form
if(m_jac_type)
std::cout << "Writing numerically estimated ";
else
std::cout << "Writing user-defined ";
std::cout << "Jacobian matrix at time t = " << t << "...\n";
const MKL_INT size =
M.ia.size() - 1; // derive the matrix size from ia index
if(size > 10000)
{
std::cout << "Jacobian::dump() -- WARNING: the size of the Jacobian "
"matrix for writting is bigger than 10000x10000.\n";
}
std::ofstream outFile;
MKL_INT ja = 0;
if(m_jac_type)
outFile.open("dump_Jacobian_" + std::to_string(m_dump_file_counter++) +
"_numerical.txt");
else
outFile.open("dump_Jacobian_" + std::to_string(m_dump_file_counter++) +
".txt");
outFile << "t=" << t;
for(MKL_INT i = 0; i < size; i++)
{
outFile << delimiter << "i=" << i;
}
outFile << '\n';
for(MKL_INT j = 0; j < size; j++)
{
MKL_INT ent = M.ia[j + 1] - M.ia[j]; // Number of entries in line j
outFile << "j=" << j << delimiter;
for(MKL_INT i = 0; i < size; i++)
{
if(M.ja[ja] == i)
{
outFile << M.A[ja++];
if(!(--ent))
break;
}
outFile << delimiter;
}
outFile << '\n';
}
outFile.close();
}
/*
* Helper function to show Jacobian structure (in sparse format)
*/
void Jacobian::print(const state_type &x, const double t)
{
if(x.size() > 1000)
{
std::cout << "\nJacobian::print() -- too much output. Skipped.\n";
return;
}
sparse_matrix_holder J;
this->operator()(J, x, t);
m_matrix_converter(J); // converts the matrix if it is in simple form
std::cout << std::right;
std::cout << "\nJacobian matrix at time t = " << t << ':';
std::cout << "\n-----------------------------------------\n";
std::cout << std::setw(7) << "i" << std::setw(16) << "J.A |"
<< std::setw(10) << "J.ja |" << std::setw(8) << "J.ia";
std::cout << "\n-----------------------------------------\n";
std::size_t size = (J.A.size() > J.ia.size()) ? J.A.size() : J.ia.size();
for(std::size_t i = 0; i < size; i++)
{
std::cout << std::setw(7) << i << ": ";
std::cout << std::setw(12);
if(i < J.A.size())
std::cout << J.A[i];
else
std::cout << ' ';
std::cout << " | " << std::setw(7);
if(i < J.ja.size())
std::cout << J.ja[i];
else
std::cout << ' ';
std::cout << " | ";
if(i < J.ia.size())
std::cout << std::setw(7) << J.ia[i];
std::cout << std::endl;
}
}
/*
* Helper function to compare two Jacobians and write the difference
*/
void Jacobian::compare(Jacobian jac, const state_type &x, const double t,
const double tol)
{
std::cout << "Jacobian::compare() -- INFO: Trying to compare two "
"Jacobians at time t = "
<< t << " and the tolerance tol = " << tol << "...\n";
sparse_matrix_holder M, J;
this->operator()(M, x, t); // calls the Jacobian matrix operator
jac(J, x, t); // external Jacobian to compare with
m_matrix_converter(M); // converts the matrix M if it is in simple form
m_matrix_converter(J); // converts the matrix J if it is in simple form
const MKL_INT size =
M.ia.size() - 1; // derive the matrix size from ia index
if((std::size_t)(size) != (J.ia.size() - 1))
{
std::cout << "Jacobian::compare() -- ERROR: the sizes of the "
"matrices do not match ('ia' indexes are different).\n";
return;
}
std::ofstream outFile;
MKL_INT ja_M = 0;
MKL_INT ja_J = 0;
outFile.open("dump_Jacobians_compare_" +
std::to_string(m_compare_file_counter++) + ".txt");
outFile << "List of differences in Jacobians for t = " << t
<< " and the tolerance tol = " << tol << ":\n";
outFile << "i" << delimiter << "j" << delimiter << "Jac_original"
<< delimiter << "Jac_reference" << delimiter << "Rel_difference"
<< '\n';
std::size_t ndiff = 0; // counts differences
for(MKL_INT j = 0; j < size; j++)
{
MKL_INT ent_M = M.ia[j + 1] - M.ia[j];
MKL_INT ent_J = J.ia[j + 1] - J.ia[j];
for(MKL_INT i = 0; i < size; i++)
{
double MA = 0.0;
double JA = 0.0;
double diff;
if((!ent_M) && (!ent_J))
break;
if((M.ja[ja_M] == i) && ent_M)
{
MA = M.A[ja_M++];
ent_M--;
}
if((J.ja[ja_J] == i) && ent_J)
{
JA = J.A[ja_J++];
ent_J--;
}
if(JA != 0.0)
{
diff = (MA - JA) / std::abs(JA);
}
else
{
diff = (MA - JA);
}
if(std::abs(diff) > tol)
{
outFile << i << delimiter << j << delimiter << MA << delimiter
<< JA << delimiter << diff << '\n';
ndiff++;
}
}
}
outFile << "Total number of differences found: " << ndiff << '\n';
std::cout << "Jacobian::compare() -- INFO: Found " << ndiff
<< " difference(s).\n";
outFile.close();
}
} // namespace daecpp_namespace_name
Boost Software License - Version 1.0 - August 17th, 2003
Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:
The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
//
// boost/assert.hpp - BOOST_ASSERT(expr)
// BOOST_ASSERT_MSG(expr, msg)
// BOOST_VERIFY(expr)
// BOOST_VERIFY_MSG(expr, msg)
// BOOST_ASSERT_IS_VOID
//
// Copyright (c) 2001, 2002 Peter Dimov and Multi Media Ltd.
// Copyright (c) 2007, 2014 Peter Dimov
// Copyright (c) Beman Dawes 2011
// Copyright (c) 2015 Ion Gaztanaga
//
// Distributed under the Boost Software License, Version 1.0.
// See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt
//
// Note: There are no include guards. This is intentional.
//
// See http://www.boost.org/libs/assert/assert.html for documentation.
//
//
// Stop inspect complaining about use of 'assert':
//
// boostinspect:naassert_macro
//
//
// BOOST_ASSERT, BOOST_ASSERT_MSG, BOOST_ASSERT_IS_VOID
//
#undef BOOST_ASSERT
#undef BOOST_ASSERT_MSG
#undef BOOST_ASSERT_IS_VOID
#if defined(BOOST_DISABLE_ASSERTS) || ( defined(BOOST_ENABLE_ASSERT_DEBUG_HANDLER) && defined(NDEBUG) )
# define BOOST_ASSERT(expr) ((void)0)
# define BOOST_ASSERT_MSG(expr, msg) ((void)0)
# define BOOST_ASSERT_IS_VOID
#elif defined(BOOST_ENABLE_ASSERT_HANDLER) || ( defined(BOOST_ENABLE_ASSERT_DEBUG_HANDLER) && !defined(NDEBUG) )
#include <boost/config.hpp> // for BOOST_LIKELY
#include <boost/current_function.hpp>
namespace boost
{
void assertion_failed(char const * expr, char const * function, char const * file, long line); // user defined
void assertion_failed_msg(char const * expr, char const * msg, char const * function, char const * file, long line); // user defined
} // namespace boost
#define BOOST_ASSERT(expr) (BOOST_LIKELY(!!(expr))? ((void)0): ::boost::assertion_failed(#expr, BOOST_CURRENT_FUNCTION, __FILE__, __LINE__))
#define BOOST_ASSERT_MSG(expr, msg) (BOOST_LIKELY(!!(expr))? ((void)0): ::boost::assertion_failed_msg(#expr, msg, BOOST_CURRENT_FUNCTION, __FILE__, __LINE__))
#else
# include <assert.h> // .h to support old libraries w/o <cassert> - effect is the same
# define BOOST_ASSERT(expr) assert(expr)
# define BOOST_ASSERT_MSG(expr, msg) assert((expr)&&(msg))
#if defined(NDEBUG)
# define BOOST_ASSERT_IS_VOID
#endif
#endif
//
// BOOST_VERIFY, BOOST_VERIFY_MSG
//
#undef BOOST_VERIFY
#undef BOOST_VERIFY_MSG
#if defined(BOOST_DISABLE_ASSERTS) || ( !defined(BOOST_ENABLE_ASSERT_HANDLER) && defined(NDEBUG) )
# define BOOST_VERIFY(expr) ((void)(expr))
# define BOOST_VERIFY_MSG(expr, msg) ((void)(expr))
#else
# define BOOST_VERIFY(expr) BOOST_ASSERT(expr)
# define BOOST_VERIFY_MSG(expr, msg) BOOST_ASSERT_MSG(expr,msg)
#endif
The MIT License (MIT)
Copyright (c) 2014 Benno Evers
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
This library does not contain any files from the matplotlib project, nor
does it make any changes to it. On the other hand, the code contained herein
is perfectly useless without a separate installation of matplotlib.
I don't know enough about US copyright law to decide whether this implies
that this library "uses" or is "based on" matplotlib.
In any case, matplotlib comes with the following license:
License agreement for matplotlib 1.4.3
1. This LICENSE AGREEMENT is between the Matplotlib Development Team (“MDT”),
and the Individual or Organization (“Licensee”) accessing and otherwise
using matplotlib software in source or binary form and its associated documentation.
2. Subject to the terms and conditions of this License Agreement, MDT hereby grants
Licensee a nonexclusive, royalty-free, world-wide license to reproduce, analyze,
test, perform and/or display publicly, prepare derivative works, distribute, and
otherwise use matplotlib 1.4.3 alone or in any derivative version, provided, however,
that MDT’s License Agreement and MDT’s notice of copyright, i.e.,
“Copyright (c) 2012-2013 Matplotlib Development Team; All Rights Reserved” are retained
in matplotlib 1.4.3 alone or in any derivative version prepared by Licensee.
3. In the event Licensee prepares a derivative work that is based on or incorporates
matplotlib 1.4.3 or any part thereof, and wants to make the derivative work available
to others as provided herein, then Licensee hereby agrees to include in any such work a
brief summary of the changes made to matplotlib 1.4.3.
4. MDT is making matplotlib 1.4.3 available to Licensee on an “AS IS” basis. MDT MAKES NO
REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. BY WAY OF EXAMPLE, BUT NOT LIMITATION,
MDT MAKES NO AND DISCLAIMS ANY REPRESENTATION OR WARRANTY OF MERCHANTABILITY OR FITNESS
FOR ANY PARTICULAR PURPOSE OR THAT THE USE OF MATPLOTLIB 1.4.3 WILL NOT INFRINGE ANY
THIRD PARTY RIGHTS.
5. MDT SHALL NOT BE LIABLE TO LICENSEE OR ANY OTHER USERS OF MATPLOTLIB 1.4.3 FOR ANY
INCIDENTAL, SPECIAL, OR CONSEQUENTIAL DAMAGES OR LOSS AS A RESULT OF MODIFYING,
DISTRIBUTING, OR OTHERWISE USING MATPLOTLIB 1.4.3, OR ANY DERIVATIVE THEREOF,
EVEN IF ADVISED OF THE POSSIBILITY THEREOF.
6. This License Agreement will automatically terminate upon a material breach of
its terms and conditions.
7. Nothing in this License Agreement shall be deemed to create any relationship of
agency, partnership, or joint venture between MDT and Licensee. This License
Agreement does not grant permission to use MDT trademarks or trade name in a
trademark sense to endorse or promote products or services of Licensee, or any
third party.
8. By copying, installing or otherwise using matplotlib 1.4.3, Licensee agrees to be
bound by the terms and conditions of this License Agreement.
/*
* The RHS class.
* This class is abstract and must be inherited.
*/
#pragma once
#include "typedefs.h"
namespace daecpp_namespace_name
{
class RHS
{
std::size_t m_dump_file_counter = 0;
public:
/*
* Takes vector x and time t and returns vector f.
* This function is pure virtual and must be overriden.
*/
virtual void operator()(const state_type &x, state_type &f,
const double t) = 0;
/*
* User-defined condition, when the solver should stop and return the
* solution at the current time step
*/
virtual bool stop_condition(const state_type &x, const double t)
{
return false;
}
/*
* Helper function to write the RHS vector to a file
*/
void dump(const state_type &x, const double t);
};
} // namespace daecpp_namespace_name
/* #undef DAE_FORTRAN_STYLE */
/* #undef DAE_SINGLE */
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment