SemiImplicit1D.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

#ifndef __solve_SemiImplicit1D_h__
#define __solve_SemiImplicit1D_h__

#include "signal/signal.h"
#include "datatool/datatool.h"
#include "shape/shape.h"

#define PARTICLE_DEBUG
//#define USE_FILTERS

/**

@page semiimplicit_1d_solver 1D solver

One main use case for this solver is for mechanical rotational systems.
However, the language of this solver is in simple clumped mass system
terminology.   So, in the interface you see "particle", "mass", "velocity",
and "force".
When using this for a mechanical rotational system, you can think of these
as "shaft", "moment of inertia", "angular velocity", and "torque".
The solve math is equivalent.

The main thing this solver has beyond a simple clumped mass solver to support
mechanical rotational systems is the idea of ratios.   Every "mass" (or shaft)
has a ratio for which it "spins" from the point of view of a use case, that
may differ from the point of view of the numerical solve.   The reason for this
is that the numerical solver can only solver in a single "spin space" for
all connected shafts.

Ratio groups ...

@verbatim

simple example

shaft <-> gear <-> shaft <-> clutch <-> shaft
  |         |        |         |
  |  <--- pair --->  |   <--- pair --->   |
  |                  |                    |
 mass <-----------> mass <-------------> mass
  |                     \               /
  |                      \             /
ratio group                ratio group

@endverbatim

Conversion functions when operating with ratios as a mechanical rotational
system.

Functions that convert values from ??? solver space to user space:
    - l_ext -- Linear (postion)                 rad
    - v_ext -- angular Velocity                 rad/s
    - f_ext -- torque (Force)                   Nm
    - k_ext -- stiffness (K)                    Nm/rad
    - d_ext -- Damping/Drag (viscous)           Nm/rad/s

Functions that convert values from user space to solver space:
    - m_int -- moment of inertia (Mass)         kgm


connect ratio into pairs?

ACCUM GUIDE

USER                    |   SOLVER
params ---- ext ------------>
                             force



compile
    precompile (add more particles into data)
        process dataset to create sim particles and constraints
    pairs
        process dataset to create collector sparse matrix, connectivity groups, etc, from pairs info
    compile (populate sparse matrix with entries necessary to accumulate)
        bake/optimize solver sparse matrix using collector sparse matrix
    reconfigure (run time dynamic change of configuration (parameters, etc))  in particular, ratios


step
    accumulate

validate

reconfigure
    reconfigure
        bake ratios in pairs down to solver per-shaft ratios.

ratios
    per shaft for solve
    per pair

*/

namespace fe
{
namespace ext
{

inline t_solve_real l_ext(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value / a_ratio;
}
inline t_solve_real v_ext(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value / a_ratio;
}
inline t_solve_real f_ext(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value * a_ratio;
}
inline t_solve_real k_ext(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value * a_ratio * a_ratio;
}
inline t_solve_real d_ext(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value * a_ratio * a_ratio;
}
inline t_solve_real m_int_deprecate(t_solve_real a_mass, t_solve_real a_ratio)
{
    return a_mass * a_ratio * a_ratio;
}


// these seem unnecessary, and tempt bugs
#if 0
inline t_solve_real l_int(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value * a_ratio;
}
inline t_solve_real v_int(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value * a_ratio;
}

inline t_solve_real f_int(t_solve_real a_value, t_solve_real a_ratio)
{
    return a_value / a_ratio;
}
#endif

/** Semi Implicit time integration

    @copydoc SemiImplicit1D_info
    */
class FE_DL_EXPORT SemiImplicit1D : public Handled<SemiImplicit1D>
{
    public:
                SemiImplicit1D(void);
virtual         ~SemiImplicit1D(void);

virtual void    compile(sp<RecordGroup> rg_input);
virtual void    initialize(sp<Scope> a_spScope);
virtual void    extract(sp<RecordGroup> rg_output);
virtual void    step(t_solve_real a_timestep, t_solve_real &a_totalConstraintForce);

virtual void    setRayleighDamping(bool a_flag) { m_rayleigh_damping = a_flag; }
virtual void    setRayleighDamping(t_solve_real a_stiffness, t_solve_real a_mass)
                {
                    m_rayleigh_stiffness = a_stiffness;
                    m_rayleigh_mass = -a_mass;
                }

virtual void    shiftDomain(t_solve_real a_period, unsigned int a_p);
virtual void    shiftDomain(void);

        void    reset(void);

        void    propogateRatiosFrom(unsigned int a_p);


        void    reconfigure(void);



    typedef unsigned int t_size;
    class t_pair : public Counted
    {
        public:
            friend SemiImplicit1D;
            t_pair(void)
            {
                m_ratio = 1.0;
            }
            t_pair(const t_size &a_first, const t_size &a_second)
            {
                first = a_first;
                second = a_second;
                m_ratio = 1.0;
            }
            void setGearRatio(int a_from, int a_to)
            {
                if (a_from < 1) { a_from = 1; }
                if (a_to < 1) { a_to = 1; }
                m_ratio = (t_solve_real)a_from/(t_solve_real)a_to;
            }
            /// This likely screws up any viable domain shifting
            void setIrrational(t_solve_real a_maybe_irrational)
            {
                m_ratio = a_maybe_irrational;
            }


            t_size          first;      // syntax to match std::pair
            t_size          second;     // syntax to match std::pair
        private:
            t_solve_real    m_ratio;
    };
    typedef std::vector< fe::sp<t_pair> > t_pairs;


    class Particle
    {
        public:
            t_solve_real        m_location;
            t_solve_real        m_adjustment;
            t_solve_real        m_velocity;
            t_solve_real        m_force;
            t_solve_real        m_force_weak;
            t_solve_real        m_force_external;
            t_solve_real        m_mass;
            t_solve_real        m_prev_location;
            t_solve_real        m_prev_velocity;
            t_solve_real        m_constraint_force;

            t_solve_real        m_ratio;
            // correct location upon changing ratios needs to be integrated on its own
            t_solve_real        m_location_ratio;
            bool                m_adjustable;
            bool                m_shift_root;

            t_pairs             m_pairs;

#ifdef PARTICLE_DEBUG
            Color                   m_color;
#endif
    };

        std::vector<Particle>   &particles(void) { return m_particles; }
        bool                    lookupIndex(unsigned int &a_particle, Record &r_particle);

    class FE_DL_EXPORT CompileMatrix
    {
        public:
                                            CompileMatrix(void);
                                            ~CompileMatrix(void);

            typedef t_solve_real                            **t_ppBlock;
            typedef std::vector<t_ppBlock>                  t_dfdx_array;
            typedef std::vector<t_ppBlock>                  t_dfdv_array;
            typedef std::pair<t_dfdx_array, t_dfdv_array>   t_entry;
            typedef std::map<unsigned int, t_entry>         t_row;

            std::vector<t_row>              m_rows;

            void                            clear(void);
            void                            setRows(unsigned int a_count);
            unsigned int                    rows(void);
            t_row                           &row(unsigned int a_index);
            t_entry                         &entry( unsigned int a_i,
                                                    unsigned int a_j);

            typedef std::set<unsigned int>                  t_nonzero_set;
            typedef std::vector< t_nonzero_set >            t_nonzero_pattern;

            void                            symbolicFill(void);
    };

    typedef std::vector< std::vector<unsigned int> >    t_shift;
    typedef std::vector< std::vector<unsigned int> >    t_neighbors;
    class Force : public Counted, public fe::CastableAs<Force>
    {
        public:
                    Force(void){ m_timestep = 0.0; }
        virtual     ~Force(void){}

        //virtual       void clear(void){}
        virtual     void accumulate(void)                                   = 0;
        virtual     bool validate(void){ return true; }

        virtual     void compile(   sp<RecordGroup> rg_input,
                                    std::vector<Particle> &a_particles,
                                    CompileMatrix &a_compileMatrix){}
        virtual     void precompile(sp<RecordGroup> rg_input) {}
        virtual     void pairs(     sp<RecordGroup> rg_input,
                                    t_pairs &a_pairs) { }
                    void setTimestep(t_solve_real a_timestep) { m_timestep = a_timestep; }
        virtual     void reconfigure(void) {}
        protected:
            t_solve_real    m_timestep;
    };

virtual void    addForce(sp<Force> a_force, bool a_add_damping=false);



    private:
        t_solve_real    m_int(t_solve_real a_mass, t_solve_real a_ratio)
        {
            return a_mass * a_ratio * a_ratio;
        }

        void reorder(std::vector<unsigned int> &a_order, t_pairs &a_pairs);

        std::vector<Particle>       m_particles;

        AsParticle1D            m_asParticle;
        AsSolverParticle1D      m_asSolverParticle;
        AsComponent             m_asComponent;
        AsForcePoint1D          m_asForcePoint;
        AsTemporal              m_asTemporal;
        AsValidate              m_asValidate;
        AsAccumulate            m_asAccumulate;
        AsClear                 m_asClear;
        AsUpdate                m_asUpdate;
#ifdef PARTICLE_DEBUG
        AsColored               m_asColored;
#endif
        AsForceFilter           m_asForceFilter;
        hp<SignalerI>           m_hpSignaler;

        unsigned int                    m_n;
        unsigned int                    m_n_sim;
        std::vector< sp<Force> >        m_forces_add_damping;
        std::vector< sp<Force> >        m_forces_as_is;

        t_solve_real                    m_dummy_block;

        sp< UpperTriangularBCRS1<t_solve_real> >            m_dfdx;
        sp< UpperTriangularBCRS1<t_solve_real> >            m_dfdv;
        sp< UpperTriangularBCRS1<t_solve_real> >            m_lhs;
        sp< UpperTriangularBCRS1<t_solve_real> >            m_lhs_snapshot;

        std::vector<t_solve_real>       m_rhs;
        std::vector<t_solve_real>       m_dv;
        std::vector<t_solve_real>       m_tmp;


        std::map<FE_UWORD, unsigned int>    m_recordToParticle;

        t_solve_real                    m_dv2dxRatio;
        t_solve_real                    m_dxImplicitness;
        t_solve_real                    m_dvImplicitness;

        unsigned int                    m_subdivcnt;
        unsigned int                    m_subdivsz;
        t_solve_real                    m_subdivmult;

        bool                            m_rayleigh_damping;
        t_solve_real                    m_rayleigh_stiffness;
        t_solve_real                    m_rayleigh_mass;
        std::vector<t_solve_real>       m_perturb;

        std::vector<t_solve_real>       m_adjustableState;

        t_shift                         m_shift;
};

} /* namespace */
} /* namespace */


#endif /* __solve_SemiImplicit1D_h__ */