Skip to content

GitLab

E
EMAM2Cpp
Commit 16ae85fa authored by Malte Heithoff's avatar Malte Heithoff
Browse files
Options

Remove DAECPP from native

parent 0b3d21e1
Pipeline #393174 failed with stage
    in 11 minutes and 51 seconds
    Expand all Hide whitespace changes
    Inline Side-by-side
    Showing with 0 additions and 5884 deletions
    native/dae-cpp_linux/include/RHS.h deleted 100644 → 0
    View file @ 0b3d21e1
    /*
    * 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
    native/dae-cpp_linux/include/cmake_config.h deleted 100644 → 0
    View file @ 0b3d21e1
    /* #undef DAE_FORTRAN_STYLE */
    /* #undef DAE_SINGLE */
    native/dae-cpp_linux/include/jacobian.h deleted 100644 → 0
    View file @ 0b3d21e1
    /*
    * Calculates Jacobian matrix.
    * If not overridden by a user, performs numerical differentiation of the RHS
    * with the given tolerance to estimate numerical Jacobian matrix.
    */
    #pragma once
    #include "typedefs.h"
    #include "RHS.h"
    namespace daecpp_namespace_name
    {
    class Jacobian
    {
    RHS &m_rhs;
    #ifdef DAE_SINGLE
    const double m_tol = 1.0e-3; // Default tolerance
    const double m_eps = 1.0e-6; // The order of the rounding unit
    #else
    const double m_tol = 1.0e-6; // Default tolerance
    const double m_eps = 1.0e-13; // The order of the rounding unit
    #endif
    std::size_t m_dump_file_counter = 0;
    std::size_t m_compare_file_counter = 0;
    int m_jac_type = 0; // This will be changed to 1
    // if numerical Jacobian is used
    /*
    * Sparse matrix converter from simple three-array format to Intel MKL
    * three array format.
    * Input: matrix holder M with simple three-array format
    * Output: matrix holder M with Intel MKL three-array format
    */
    void m_matrix_converter(daecpp::sparse_matrix_holder &M);
    public:
    explicit Jacobian(RHS &rhs) : m_rhs(rhs) {}
    Jacobian(RHS &rhs, const double tol) : m_rhs(rhs), m_tol(tol)
    {
    // TODO: Check user's tol parameter. Too small tol may lead to crash.
    }
    /*
    * Can be overriden to provide analytical Jacobian
    */
    virtual void operator()(sparse_matrix_holder &J, const state_type &x,
    const double t);
    /*
    * Helper function to show Jacobian structure on screen (in sparse format)
    */
    void print(const state_type &x, const double t);
    /*
    * Helper function to write Jacobian matrix to a file (in dense format)
    */
    void dump(const state_type &x, const double t);
    /*
    * Helper function to compare two Jacobians and write the differences.
    * Comparison will be made with the external Jacobian jac (usually,
    * numerical Jacobian) using vector x at time t with the given tolerance.
    */
    void compare(Jacobian jac, const state_type &x, const double t,
    const double tol);
    };
    } // namespace daecpp_namespace_name
    native/dae-cpp_linux/include/mass_matrix.h deleted 100644 → 0
    View file @ 0b3d21e1
    /*
    * Mass matrix classes
    */
    #pragma once
    #include "typedefs.h"
    namespace daecpp_namespace_name
    {
    /*
    * Parent class. This class is abstract and must be inherited.
    */
    class MassMatrix
    {
    /*
    * Sparse matrix converter from simple three-array format to Intel MKL
    * three array format.
    * Input: matrix holder M with simple three-array format
    * Output: matrix holder M with Intel MKL three-array format
    */
    void m_matrix_converter(daecpp::sparse_matrix_holder &M);
    public:
    /*
    * The matrix should be defined in sparse format,
    * see three array sparse format decription on
    * https://software.intel.com/en-us/mkl-developer-reference-c-sparse-blas-csr-matrix-storage-format
    *
    * The Mass matrix is static, i.e. it does not depend on time t and
    * the vector x
    *
    * This function is pure virtual and must be overriden.
    */
    virtual void operator()(sparse_matrix_holder &M) = 0;
    /*
    * Helper function to write the Mass matrix to a file
    */
    void dump();
    };
    /*
    * Helper child class to create identity Mass matrix of size N
    */
    class MassMatrixIdentity : public MassMatrix
    {
    const MKL_INT m_N;
    public:
    MassMatrixIdentity(const MKL_INT N) : daecpp::MassMatrix(), m_N(N) {}
    void operator()(daecpp::sparse_matrix_holder &M);
    };
    } // namespace daecpp_namespace_name
    native/dae-cpp_linux/include/solver.h deleted 100644 → 0
    View file @ 0b3d21e1
    /*
    * The main solver class definition
    */
    #pragma once
    #include "typedefs.h"
    #include "RHS.h"
    #include "jacobian.h"
    #include "mass_matrix.h"
    #include "solver_options.h"
    #include "time_integrator.h"
    namespace daecpp_namespace_name
    {
    class Solver
    {
    RHS &m_rhs; // RHS
    Jacobian &m_jac; // Jacobian matrix
    MassMatrix &m_mass; // Mass matrix
    TimeIntegrator *m_ti; // Pointer to the time integrator
    struct m_iterator_state_struct // Keeps the current time layer state
    {
    double t; // current time
    double dt[2]; // current and previous time steps
    double dt_eval; // new time step
    int current_scheme; // current BDF order
    int step_counter_local; // local time step counter
    bool final_time_step; // do final time step
    } m_iterator_state;
    MKL_INT m_size; // System size
    std::size_t m_steps = 0; // Total time iteration counter
    std::size_t m_calls = 0; // Total linear algebra solver calls counter
    // Count the number of the DAE solver calls (for output)
    std::size_t m_dae_solver_calls = 0;
    // Timers
    double m_timer_lin = 0;
    double m_timer_rhs = 0;
    double m_timer_jac = 0;
    double m_timer_tot = 0;
    // Contains a few latest successful time steps for the time integrator
    state_type_matrix m_x_prev;
    // Intel MKL PARDISO control parameters
    MKL_INT m_phase; // Current phase of the solver
    MKL_INT m_maxfct = 1; // Maximum number of numerical factorizations
    MKL_INT m_mnum = 1; // Which factorization to use
    MKL_INT m_mtype = 11; // Real unsymmetric matrix
    MKL_INT m_nrhs = 1; // Number of right hand sides
    MKL_INT m_msglvl = 0; // Print statistical information
    MKL_INT m_error = 0; // Error flag
    // Intel MKL PARDISO sparse matrix indeces
    MKL_INT *m_ia = nullptr;
    MKL_INT *m_ja = nullptr;
    // Intel MKL PARDISO vectors and sparse matrix non-zero elements
    float_type *m_mkl_a = nullptr;
    float_type *m_mkl_b = nullptr;
    float_type *m_mkl_x = nullptr;
    // Intel MKL PARDISO internal solver memory pointer pt,
    // 32-bit: int pt[64]; 64-bit: long int pt[64]
    // or void *pt[64] should be OK on both architectures
    void *m_pt[64];
    // Intel MKL PARDISO auxiliary variables
    double m_ddum; // Double dummy
    MKL_INT m_idum; // Integer dummy
    // Intel MKL PARDISO iparm parameter
    MKL_INT m_iparm[64];
    // Simple yet efficient Adaptive Time Stepping
    int m_adaptive_time_stepping(state_type &x, const state_type_matrix &x_prev,
    int iter);
    // Scrapes the current time iteration and decreases the time step
    // Return -1 in case the time step is below dt_min
    int m_reset_ti_state(state_type &x, const state_type_matrix &x_prev);
    // Updates time integrator scheme when the time step changes
    int m_reset_ti_scheme();
    // Increases the time step
    void m_increase_dt();
    // Decreases the time step
    void m_decrease_dt();
    // Checks if dt is within the interval defined in solver_options.h
    int m_check_dt();
    // Checks PARDISO solver error messages
    void m_check_pardiso_error(MKL_INT err);
    protected:
    /*
    * Expose solver options for the observer in the children classes
    */
    SolverOptions &m_opt; // Solver options
    public:
    /*
    * Receives user-defined RHS, Jacobian, Mass matrix and solver options.
    * Defined in solver.cpp
    */
    Solver(RHS &rhs, Jacobian &jac, MassMatrix &mass, SolverOptions &opt);
    /*
    * Releases memory. Defined in solver.cpp.
    */
    virtual ~Solver();
    /*
    * Integrates the system of DAEs on the interval t = [t0; t1] and returns
    * result in the array x. Parameter t0 can be overriden in the solver
    * options (t0 = 0 by default).
    * The data stored in x (initial conditions) will be overwritten.
    * Returns 0 in case of success or error code if integration is failed.
    */
    int operator()(state_type &x, double &t1);
    /*
    * Virtual Observer. Called by the solver every time step.
    * Receives current solution vector and the current time.
    * Can be overriden by a user to get intermediate results (for example,
    * for plotting).
    */
    virtual void observer(state_type &x, const double t)
    {
    return; // It does nothing by default
    }
    };
    } // namespace daecpp_namespace_name
    native/dae-cpp_linux/include/solver_options.h deleted 100644 → 0
    View file @ 0b3d21e1
    /*
    * Contains public solver parameters
    */
    #pragma once
    #include "typedefs.h"
    namespace daecpp_namespace_name
    {
    #define BDF_MAX_ORDER 6 // Max BDF scheme order currently implemented
    class SolverOptions
    {
    public:
    // You can control the parallel execution of the solver by explicitly
    // setting the MKL_NUM_THREADS environment variable. If fewer OpenMP threads
    // are available than specified, the execution may slow down instead of
    // speeding up. If MKL_NUM_THREADS is not defined, then the solver uses all
    // available processors.
    // In order to control the parallel execution of the numerical Jacobian, use
    // OMP_NUM_THREADS environment variable (uses all cores by default).
    // Perform Jacobian update, Reordering, Symbolic and Numerical Factorization
    // every Newton iteration. Changing to 'false' can increase speed but also
    // can lead to instability.
    bool fact_every_iter = true;
    // If fact_every_iter = false, update Jacobian every fact_iter Newton
    // iterations
    int fact_iter = 15;
    // Order of BDF implicit numerical integration method:
    // 1 - first order BDF, 2 - BDF-2, ..., 6 - BDF-6
    // Default is BDF-2 since it fully supports variable time stepping
    int bdf_order = 2;
    // Time stepping algorithm:
    // 1 - Stability-based Simple Adaptive Time Stepping (S-SATS),
    // 2 - Variability-based Simple Adaptive Time Stepping (V-SATS) from
    // https://www.sciencedirect.com/science/article/pii/S0377042705005534
    int time_stepping = 2;
    // Maximum number of Newton iterations. If the Newton method fails to
    // converge after max_Newton_iter iterations, the solver reduces time step
    // and tries to make the current step again.
    int max_Newton_iter = 15;
    #ifdef DAE_SINGLE
    double atol = 1.0e-3; // Absolute tolerance for the Newton algorithm
    double rtol = 1.0e-6; // Relative tolerance for the Newton algorithm
    double dt_eps_m = 1.0e-6; // The order of the rounding unit
    double value_max = 1.0e20; // Solution shouldn't be higher than this
    #else
    double atol = 1.0e-6; // Absolute tolerance for the Newton algorithm
    double rtol = 1.0e-6; // Relative tolerance for the Newton algorithm
    double dt_eps_m = 1.0e-14; // The order of the rounding unit
    double value_max = 1.0e100; // Solution shouldn't be higher than this
    #endif
    // Initial time step
    double dt_init = 0.1;
    // Initial integration time t0
    double t0 = 0.0;
    // Minimum time step
    double dt_min = dt_eps_m;
    // Maximum time step
    double dt_max = 1.0 / dt_eps_m;
    // Verbosity level of the solver:
    // 0 - silent, 1 - basic information, 2 - time stepping info, 3 - all info
    int verbosity = 2;
    // Simple Adaptive Time Stepping options
    int dt_increase_threshold = 2; // Time step amplification threshold
    // (S-SATS only)
    int dt_decrease_threshold = 7; // Time step reduction threshold
    // (S-SATS only)
    double dt_increase_factor = 2.0; // Time step amplification factor
    double dt_decrease_factor = 2.0; // Time step reduction factor
    double dt_eta_min = 0.05; // Monitor function lower threshold (V-SATS only)
    double dt_eta_max = 0.5; // Monitor function higher threshold (V-SATS only)
    // Try to roll back and reduce the time step if Newton iterations diverged.
    // Otherwise stop with error message.
    bool redo_newton = false;
    // If Newton method fails to converge within 'max_Newton_iter' iterations
    // 'newton_failed_attempts' times in a row, the solver will try to update
    // Jacobian every single iteration next time step.
    int newton_failed_attempts = 3;
    // 1 - V-SATS will use NORM_infinity to estimate solution variability,
    // 2 - V-SATS will use NORM_2 (default)
    int vsats_norm = 2;
    // Intel MKL PARDISO parameters (iparam). More about iparam:
    // https://software.intel.com/en-us/mkl-developer-reference-c-pardiso-iparm-parameter
    // iparm[3]: Controls preconditioned CGS.
    // 0 - The factorization is always computed as required.
    // 31 - LU-preconditioned CGS iteration with a stopping criterion of 1.0E-3.
    // 61 - LU-preconditioned CGS iteration with a stopping criterion of 1.0E-6.
    MKL_INT preconditioned_CGS = 0;
    // iparm[7]: Maximum number of iterative refinement steps.
    // 0 - The solver automatically performs two steps of iterative refinement.
    // >0 - Maximum number of iterative refinement steps that the solver
    // performs.
    MKL_INT refinement_steps = 2;
    // iparm[12]: Maximum weighted matching algorithm.
    // 0 - OFF
    // 1 - ON (may conflict with parallel_fact_control)
    MKL_INT matching_alg = 1;
    // iparm[23]: Parallel factorization control.
    // 0 - Intel MKL PARDISO uses the classic algorithm for factorization.
    // 1 - Two-level factorization algorithm. This algorithm generally improves
    // scalability in case of parallel factorization on many OpenMP threads
    // (more than eight).
    // 10 - Improved two-level factorization algorithm.
    MKL_INT parallel_fact_control = 10;
    // iparm[24]: Parallel forward/backward solve control.
    // 0 - Intel MKL PARDISO uses the sequential forward and backward solve.
    // 1 - Parallel algorithm for the solve step (in-core mode only).
    // 2 - If there is only one right hand side, Intel MKL PARDISO performs
    // parallelization. If there are multiple right hand sides, PARDISO
    // performs parallel forward and backward substitution via matrix.
    MKL_INT parallel_solve_control = 0;
    SolverOptions() = default;
    // Initialises Intel MKL PARDISO parameters (iparam) array
    void set_iparm_for_pardiso(MKL_INT *iparm);
    // Checks correctness of the solver parameters
    // TODO: should return error code
    void check_options();
    };
    } // namespace daecpp_namespace_name
    native/dae-cpp_linux/include/time_integrator.h deleted 100644 → 0
    View file @ 0b3d21e1
    /*
    * Numerical time integrator class definition
    */
    #pragma once
    #include <mkl_spblas.h>
    #include "typedefs.h"
    #include "RHS.h"
    #include "jacobian.h"
    #include "mass_matrix.h"
    #include "solver_options.h"
    namespace daecpp_namespace_name
    {
    class TimeIntegrator
    {
    RHS &m_rhs;
    Jacobian &m_jac;
    MassMatrix &m_mass;
    SolverOptions &m_opt;
    // clang-format off
    const double BDF_COEF[6][8] = {
    { 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0},
    { 3.0, -4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 2.0},
    { 11.0, -18.0, 9.0, -2.0, 0.0, 0.0, 0.0, 6.0},
    { 25.0, -48.0, 36.0, -16.0, 3.0, 0.0, 0.0, 12.0},
    {137.0, -300.0, 300.0, -200.0, 75.0, -12.0, 0.0, 60.0},
    {147.0, -360.0, 450.0, -400.0, 225.0, -72.0, 10.0, 60.0}};
    const double ALPHA_COEF[6] =
    {1.0, 1.5, 11.0/6.0, 25.0/12.0, 137.0/60.0, 2.45};
    // clang-format on
    // The first time step will be performed using BDF-1
    int m_scheme = 1;
    // Total time spent to estimate Jacobian, sec.
    double m_jac_time = 0.0;
    // Total time spent to calculate the RHS, sec.
    double m_rhs_time = 0.0;
    // Temporary Jacobian matrix holder
    sparse_matrix_holder m_J;
    // Mass matrix container
    sparse_matrix_holder m_M;
    // Descriptor of sparse mass matrix properties
    struct matrix_descr m_descrA;
    // Structure with sparse mass matrix stored in CSR format
    sparse_matrix_t m_csrA;
    /*
    * Sparse matrix converter from simple three-array format to Intel MKL
    * three array format.
    * Input: matrix holder M with simple three-array format
    * Output: matrix holder M with Intel MKL three-array format
    */
    void m_matrix_converter(daecpp::sparse_matrix_holder &M);
    /*
    * Sparse matrix checker
    */
    int m_matrix_checker(sparse_matrix_holder &A, MKL_INT size);
    /*
    * Performs matrix-matrix addition: C = alpha*A + B.
    * Replaces deprecated Intel MKL mkl_dcsradd() function.
    */
    void m_matrix_add(const float_type alpha, const sparse_matrix_holder &A,
    const sparse_matrix_holder &B, sparse_matrix_holder &C);
    public:
    TimeIntegrator(RHS &rhs, Jacobian &jac, MassMatrix &mass,
    SolverOptions &opt);
    ~TimeIntegrator() { mkl_sparse_destroy(m_csrA); }
    void set_scheme(int scheme) { m_scheme = scheme; }
    void reset_jac_time() { m_jac_time = 0.0; }
    void reset_rhs_time() { m_rhs_time = 0.0; }
    double get_jac_time() { return m_jac_time; }
    double get_rhs_time() { return m_rhs_time; }
    void integrate(sparse_matrix_holder &J, state_type &b, const state_type &x,
    const state_type_matrix &x_prev, const double t,
    const double dt[], const bool do_jac);
    };
    } // namespace daecpp_namespace_name
    native/dae-cpp_linux/include/typedefs.h deleted 100644 → 0
    View file @ 0b3d21e1
    /*
    * The main types and definitions of the solver
    */
    #pragma once
    #include <mkl_types.h>
    #include <cstddef>
    #include <vector>
    #include "cmake_config.h"
    #define daecpp_namespace_name daecpp
    namespace daecpp_namespace_name
    {
    #ifdef DAE_SINGLE
    typedef float float_type;
    #else
    typedef double float_type;
    #endif
    typedef std::vector<float_type> state_type;
    typedef std::vector<MKL_INT> vector_type_int;
    typedef std::vector<std::vector<float_type>> state_type_matrix;