tetra.cpp

/**
  Elementary matrix Calculation for a tetrahedron element
 */

#include <set>

#include "config.h"  // to get gamma0 constant
#include "tetra.h"
#include "time_integration.h"
#include "facette.h"

using namespace Tetra;
using namespace Nodes;

void prm::infos(void)
    {
    using std::cout;
    cout << "volume region name = " << regName;
    if(Ms > 0)
        {
        cout << ", magnetic.\n";
        cout << "alpha_LLG = " << alpha_LLG << std::endl;
        cout << "A = " << A << std::endl;
        cout << "Ms = " << Ms << std::endl;

        if (K != 0)
            {
            cout << "K*uk =" << K << "*[ " << uk << "]\n";
            }

        if (K3 != 0)
            {
            cout << "K3 = " << K3 << "; ex=[ " << ex << "], ey=[" << ey << "], ez=[" << ez << "]\n";
            }
        cout << "l_sd = " << lsd << std::endl;
        }
    else
        {
        cout << ", non magnetic metal.\n";
        cout << "l_sd = " << lsd << std::endl;
        cout << "l_sf = " << lsf << std::endl;
        cout << "polarization rate = " << P << std::endl;
        cout << "density of states at Fermi level = " << N0 << std::endl;
        cout << "sigma = " << sigma << std::endl;
        }
    }

Eigen::Matrix<double,NPI,1> Tetra::calc_alpha_eff(const double dt, const double alpha,
                                                      Eigen::Ref<Eigen::Matrix<double,NPI,1>> uHeff)
    {
    double reduced_dt = gamma0 * dt;
    Eigen::Matrix<double,NPI,1> a_eff;
    a_eff.setConstant(alpha);
    const double r = 0.1;  // what is that constant, where does it come from ?
    const double M = 2. * alpha * r / reduced_dt;

    for(int npi=0;npi<NPI;npi++)
        {
        double h = uHeff(npi);
        if (h > 0.)
            {
            if (h > M)
                a_eff(npi) = alpha + reduced_dt / 2. * M;
            else
                a_eff(npi) = alpha + reduced_dt / 2. * h;
            }
        else
            {
            if (h < -M)
                a_eff(npi) = alpha / (1. + reduced_dt / (2. * alpha) * M);
            else
                a_eff(npi) = alpha / (1. - reduced_dt / (2. * alpha) * h);
            }
        }
    return a_eff;
    }

Eigen::Matrix<double,Nodes::DIM,NPI> Tetra::calc_gradV(Tet const &tet, std::vector<double> &V)
    {
    Eigen::Matrix<double,N,1> _V;
    for (int i = 0; i < N; i++)
        { _V[i] = V[tet.ind[i]]; }
    Eigen::Matrix<double,Nodes::DIM,NPI> _gradV;
    for (int npi = 0; npi < NPI; npi++)
        {
        Eigen::Vector3d v(0,0,0);
        for (int i = 0; i < N; i++)
            { v += _V[i] * tet.da.row(i); }
        _gradV.col(npi) = v;
        }
    return _gradV;
    }

Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> Tetra::calc_Hst(Tetra::Tet const &tet,
        const double prefactor, std::vector<Eigen::Vector3d> &s)
    {
    Eigen::Matrix<double,Nodes::DIM,Tetra::N> s_nod;
    for(size_t ie=0;ie<Tetra::N;ie++)
        {
        size_t i = tet.ind[ie];
        for(size_t d=0;d<Nodes::DIM;d++)
            { s_nod(d,ie)= s[i][d]; }
        }

    Eigen::Matrix<double,Nodes::DIM,Tetra::NPI> Hst = prefactor * s_nod * Tetra::eigen_a;
    return Hst;
    }

void Tet::lumping(Eigen::Ref<Eigen::Matrix<double,NPI,1>> alpha_eff, double prefactor,
                  Eigen::Ref<Eigen::Matrix<double,3*N,3*N>> AE ) const
    {
    // contrib is alpha contribution to the diagonal of AE; to help stabilizing the scheme, alpha is
    // modified
    Eigen::Matrix<double,N,1> contrib = eigen_a * weight.cwiseProduct(alpha_eff);
    Eigen::Matrix<double,N,N> exch_block = calcDiagBlock(prefactor,contrib);
    AE.block<N,N>(0,0) += exch_block;
    AE.block<N,N>(N,N) += exch_block;
    AE.block<N,N>(2*N,2*N) += exch_block;
      
    // off diagonal blocks are filled with magnetization components weighted products
    Eigen::Matrix<double,N,1> a_w = eigen_a * weight;
    Eigen::Matrix<double,N,1> diag = a_w.cwiseProduct(calcOffDiagBlock(IDX_X));
    AE.block<N,N>(N,2*N).diagonal() -= diag;
    AE.block<N,N>(2*N,N).diagonal() += diag;
    diag = a_w.cwiseProduct(calcOffDiagBlock(IDX_Y));
    AE.block<N,N>(0,2*N).diagonal() += diag;
    AE.block<N,N>(2*N,0).diagonal() -= diag;
    diag = a_w.cwiseProduct(calcOffDiagBlock(IDX_Z));
    AE.block<N,N>(0,N).diagonal() -= diag;
    AE.block<N,N>(N,0).diagonal() += diag;
    }

Eigen::Matrix<double,N,N> Tet::calcDiagBlock(const double c, Eigen::Matrix<double,N,1> &x) const
    {
    Eigen::Matrix<double,N,N> result = da*da.transpose();
    result *= c*weight.sum();
    result.diagonal() += x;
    return result;
    }

Eigen::Matrix<double,N,1> Tet::calcOffDiagBlock(const Nodes::index idx) const
    {
    Eigen::Matrix<double,N,1> result;
    for (int i = 0; i < N; i++)
        result[i] = getNode(i).get_u(Nodes::CURRENT)[idx];
    return result;
    }

void Tet::add_drift_BE(double alpha, double s_dt, double Vdrift,
                       Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> U,
                       Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> V, 
                       Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dUd_,
                       Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dVd_,
                       Eigen::Ref<Eigen::Matrix<double,DIM,N>> BE) const
    {  // the artificial drift from eventual recentering is along x,y or z
    Eigen::Matrix<double,DIM,N> interim;

    for (int npi = 0; npi < NPI; npi++)
        {
        for (int i = 0; i < N; i++)
            {
            interim.col(i) = a[i][npi]*( alpha*dUd_.col(npi) + U.col(npi).cross(dUd_.col(npi))
                    + s_dt*(alpha*dVd_.col(npi) + U.col(npi).cross(dVd_.col(npi))
                        + V.col(npi).cross(dUd_.col(npi))) );
            }
        BE += Vdrift*weight[npi]*interim;
        }
    }

Eigen::Matrix<double,NPI,1> Tet::calc_aniso_uniax(Eigen::Ref<const Eigen::Vector3d> uk,
        const double Kbis,
        const double s_dt,
        Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> U,
        Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> V,
        Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> H_aniso) const
    {
    for(int npi = 0;npi<NPI;npi++)
        { H_aniso.col(npi) += (Kbis * uk.dot( U.col(npi) + s_dt * V.col(npi))) * uk; }
    return Kbis * (U.transpose() * uk).array().square();
    }

Eigen::Matrix<double,NPI,1> Tet::calc_aniso_cub(Eigen::Ref<const Eigen::Vector3d> ex,
                           Eigen::Ref<const Eigen::Vector3d> ey,
                           Eigen::Ref<const Eigen::Vector3d> ez,
                           const double K3bis, const double s_dt,
                           Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> U,
                           Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> V,
                           Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> H_aniso) const
    {
    Eigen::Matrix<double,NPI,1> result;
    for(int npi = 0;npi<NPI;npi++)
        {
        Eigen::Vector3d uk_u = Eigen::Vector3d(ex.dot(U.col(npi)),
                ey.dot(U.col(npi)), ez.dot(U.col(npi)));
        Eigen::Vector3d uk_v = Eigen::Vector3d(ex.dot(V.col(npi)),
                ey.dot(V.col(npi)), ez.dot(V.col(npi)));
        Eigen::Vector3d uk_uuu = uk_u.unaryExpr( [](double x){ return x*(1.0 - x*x);} );
        Eigen::Vector3d tmp = uk_v.cwiseProduct(ex);

        H_aniso.col(npi) += -K3bis * (uk_uuu(0) * ex + uk_uuu(1) * ey + uk_uuu(2) * ez
                  + s_dt * tmp.cwiseProduct( Eigen::Vector3d(1, 1, 1)
                      - 3*uk_u.cwiseProduct(uk_u) ));

        result[npi] = uk_u.dot(uk_uuu);
        }
    return -K3bis*result;
    }

void Tet::integrales(Tetra::prm const &param, timing const &prm_t,
                     std::function< Eigen::Matrix<double,DIM,NPI> (void)> calc_Hext,
                     Nodes::index idx_dir, double Vdrift)
    {
    const double alpha = param.alpha_LLG;
    const double Ms = param.Ms;
    const double Abis = 2.0 * param.A / (mu0*Ms);
    const double dt = prm_t.get_dt();
    const double s_dt = THETA * dt * gamma0;  // theta from theta scheme in config.h.in

    /*-------------------- INTERPOLATION --------------------*/
    Eigen::Matrix<double,DIM,NPI> U,dUdx,dUdy,dUdz;
    interpolation(Nodes::get_u<CURRENT>, U, dUdx, dUdy, dUdz);
    Eigen::Matrix<double,Nodes::DIM,N> v_nod;
    for (int i = 0; i < N; i++)
        { v_nod.col(i) = Nodes::get_v<CURRENT>(getNode(i)); }
    Eigen::Matrix<double,DIM,NPI> V = v_nod * eigen_a;
    Eigen::Matrix<double,DIM,NPI> Hd;
    interpolation_field(Nodes::get_phi<CURRENT>, Hd);
    Eigen::Matrix<double,DIM,NPI> Hv;
    interpolation_field(Nodes::get_phiv<CURRENT>, Hv);
    /*-------------------- END INTERPOLATION ----------------*/
    Eigen::Matrix<double,NPI,1> uHeff = -Abis *( dUdx.colwise().squaredNorm()
            + dUdy.colwise().squaredNorm() + dUdz.colwise().squaredNorm());
    Eigen::Matrix<double,DIM,NPI> H_aniso;
    H_aniso.setZero();

    if(param.K != 0)
        {
        double Kbis = 2.0 * param.K/(mu0*Ms);
        uHeff += calc_aniso_uniax(param.uk, Kbis, s_dt/gamma0, U, V, H_aniso);
        }
    if(param.K3 != 0)
        {
        double K3bis = 2.0 * param.K3/(mu0*Ms);
        uHeff += calc_aniso_cub(param.ex, param.ey, param.ez, K3bis, s_dt/gamma0, U, V, H_aniso);
        }

    Eigen::Matrix<double,DIM,NPI> Heff = Hd + calc_Hext();
    Eigen::Matrix<double,DIM,NPI> H = Heff; // we need Hd+Hext for future computations
    Eigen::Matrix<double,DIM,NPI> Hst; // Hst is an effective field for spin accumulation
    Hst.setZero();
    extraField(Hst); // carefull, extrafield do a += like operation on Hst, not a =

    Heff += Hst;
    uHeff += (U.cwiseProduct(Heff)).colwise().sum();//dot product on each col of U and Heff

    Eigen::Matrix<double,NPI,1> a_eff = calc_alpha_eff(dt, alpha, uHeff);
    
    Eigen::Matrix<double,3*N,3*N> AE;
    AE.setZero();
    lumping(a_eff, prm_t.prefactor * s_dt * Abis, AE);

    Eigen::Matrix<double,2*N,3*N> P;
    buildMatP(P);

    // Perm is a permutation that swaps the block-lines of P:
    //   Perm × P = ⎛Eqx Eqy Eqz⎞
    //              ⎝Epx Epy Epz⎠
    // This permutation makes the global linear system better conditioned.
    Eigen::PermutationMatrix<8> Perm;
    Perm.indices() = {4, 5, 6, 7, 0, 1, 2, 3};

    /*--------------------   PROJECTION: AE->Kp   --------------------*/
    Kp = Perm * P * AE * P.transpose();// with MKL installed this operation should call dgemm_direct
    
    /*********************   calculations on BE   *********************/
    Eigen::Matrix<double,DIM,N> BE;
    BE.setZero();
    if (idx_dir != IDX_UNDEF)
        {
        Eigen::Matrix<double,DIM,NPI>  dVd_dir = v_nod * (da.col(idx_dir)).replicate(1,NPI);
        if (idx_dir == IDX_Z)
            add_drift_BE(alpha, s_dt, Vdrift, U, V, dUdz, dVd_dir, BE);
        else if (idx_dir == IDX_Y)
            add_drift_BE(alpha, s_dt, Vdrift, U, V, dUdy, dVd_dir, BE);
        else if (idx_dir == IDX_X)
            add_drift_BE(alpha, s_dt, Vdrift, U, V, dUdx, dVd_dir, BE);
        }
    H += H_aniso + (s_dt / gamma0) * Hv;
    
    for (int npi = 0; npi < NPI; npi++)
        {
        const double w = weight[npi];
        double scal_Hst_u = (Hst.col(npi)).dot(U.col(npi));
        for (int i = 0; i < N; i++)
            {
            const double ai_w = w*a[i][npi];
            BE.col(i) -= w*Abis*(da(i,0)*dUdx.col(npi) + da(i,1)*dUdy.col(npi)
                    + da(i,2)*dUdz.col(npi));// exchange
            BE.col(i) += ai_w*(H.col(npi) + Hst.col(npi) ); // Hst contribution to BE
            BE.col(i) -= ai_w*scal_Hst_u*s_dt*V.col(npi); // spin acc second order contrib
            }
        }

    /*--------------------   PROJECTION: BE->Lp   --------------------*/
    Lp = Perm * P * BE.reshaped<Eigen::RowMajor>();
    }

double Tet::exchangeEnergy(Tetra::prm const &param,
                           Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dudx,
                           Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dudy,
                           Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dudz) const
    {// partial reduction on columns with colwise()
    Eigen::Matrix<double,NPI,1> dens = dudx.colwise().squaredNorm()
                                     + dudy.colwise().squaredNorm()
                                     + dudz.colwise().squaredNorm();
    return param.A * weight.dot(dens);
    }

double Tet::uniaxialAnisotropyEnergy(Tetra::prm const &param,
        Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> u) const
    {
    Eigen::Matrix<double,NPI,1> dens;

    for (int npi = 0; npi < NPI; npi++)
        { dens[npi] = sq( param.uk.dot( u.col(npi) )); }
    return -param.K*weight.dot(dens);
    }

double Tet::cubicAnisotropyEnergy(Tetra::prm const &param,
        Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> u) const
    {
    Eigen::Matrix<double,NPI,1> dens;

    for (int npi = 0; npi < NPI; npi++)
        {
        Eigen::Vector3d m = u.col(npi);
        // cosinus directeurs
        double al0 = m.dot(param.ex);
        double al1 = m.dot(param.ey);
        double al2 = m.dot(param.ez);
        dens[npi] = sq(al0 * al1) + sq(al1 * al2) + sq(al2 * al0);
        }
    return param.K3*weight.dot(dens);
    }

Eigen::Matrix<double,Tetra::NPI,1> Tet::charges(const double &Ms,
                  std::function<Eigen::Vector3d(const Nodes::Node&)> getter) const
    {
    Eigen::Matrix<double,Nodes::DIM,N> vec_nod;
    for (int i = 0; i < N; i++)
        { vec_nod.col(i) = getter(getNode(i)); }

    double dud_sum = (vec_nod.row(IDX_X)).dot( da.col(IDX_X) )
                   + (vec_nod.row(IDX_Y)).dot( da.col(IDX_Y) )
                   + (vec_nod.row(IDX_Z)).dot( da.col(IDX_Z) );
    
    return -Ms*weight*dud_sum;
    }

double Tet::demagEnergy(Tetra::prm const &param,
                       Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dudx,
                       Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dudy,
                       Eigen::Ref<Eigen::Matrix<double,DIM,NPI>> dudz,
                       Eigen::Ref<Eigen::Matrix<double,NPI,1>> phi) const
    {
    Eigen::Matrix<double,NPI,1> dens;

    for (int npi = 0; npi < NPI; npi++)
        { dens[npi] = (dudx(0,npi) + dudy(1,npi) + dudz(2,npi)) * phi[npi]; }
    return -0.5*mu0*param.Ms*weight.dot(dens);
    }

double Tet::zeemanEnergy(Tetra::prm const &param, Eigen::Ref<Eigen::Vector3d> const Hext,
                        Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> const u) const
    {
    Eigen::Matrix<double,NPI,1> dens = u.transpose() * Hext;

    return -mu0*param.Ms*weight.dot(dens);
    }

double Tet::zeemanEnergy(Tetra::prm const &param, double fieldAmp,
        std::vector<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> &spaceField,
        Eigen::Ref<Eigen::Matrix<double,Nodes::DIM,Tetra::NPI>> const u) const
    {
    Eigen::Matrix<double,DIM,NPI> Hext = spaceField[idx];
    Eigen::Matrix<double,NPI,1> dens = u.cwiseProduct(Hext).colwise().sum(); // dot product column
                                                                             // to column

    return -mu0*param.Ms*fieldAmp*weight.dot(dens);
    }

double Tet::Jacobian(Eigen::Ref<Eigen::Matrix3d> J)
    {
    Eigen::Vector3d p0p1 = getNode(1).p - getNode(0).p;
    Eigen::Vector3d p0p2 = getNode(2).p - getNode(0).p;
    Eigen::Vector3d p0p3 = getNode(3).p - getNode(0).p;
    J(0,0) = p0p1.x();
    J(0,1) = p0p2.x();
    J(0,2) = p0p3.x();
    J(1,0) = p0p1.y();
    J(1,1) = p0p2.y();
    J(1,2) = p0p3.y();
    J(2,0) = p0p1.z();
    J(2,1) = p0p2.z();
    J(2,2) = p0p3.z();
    return J.determinant();
    }

double Tet::calc_vol(void) const
    {
    Eigen::Vector3d p0p1 = getNode(1).p - getNode(0).p;
    Eigen::Vector3d p0p2 = getNode(2).p - getNode(0).p;
    Eigen::Vector3d p0p3 = getNode(3).p - getNode(0).p;

    return p0p1.dot(p0p2.cross(p0p3))/6.0;
    }