<commit_before>/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011-2015 FlowKit Sarl
 * Route d'Oron 2
 * 1010 Lausanne, Switzerland
 * E-mail contact: <EMAIL>
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/* <NAME> contributed this code.
 */

/** \file
 * Template specializations for some computationally intensive LB 
 * functions of the header file mrtTemplates.h, for the D3Q19 grid.
 */

#ifndef MRT_TEMPLATES_3D_H
#define MRT_TEMPLATES_3D_H

#include "core/globalDefs.h"

namespace plb {

// Efficient specialization for D3Q19 lattice
template<typename T> struct mrtTemplatesImpl<T, descriptors::MRTD3Q19DescriptorBase<T> > {
    
    typedef descriptors::D3Q19DescriptorBase<T> Descriptor; 
    typedef descriptors::MRTD3Q19DescriptorBase<T> MRTDescriptor; 

    /// Computation of all equilibrium distribution (in moments space)
    static void computeEquilibriumMoments( Array<T,Descriptor::q>& momentsEq,
                                    T rhoBar, Array<T,3> const& j, T jSqr )
    {
        T invRho = Descriptor::invRho(rhoBar);
        
        momentsEq[0] = rhoBar;
        momentsEq[1] = (T)19*jSqr*invRho-(T)11*rhoBar;
        momentsEq[2] = -(T)5.5*jSqr*invRho+(T)3*rhoBar;
        momentsEq[3] = j[0];
        momentsEq[4] = -((T)2/(T)3)*j[0];
        momentsEq[5] = j[1];
        momentsEq[6] = -((T)2/3)*j[1];
        momentsEq[7] = j[2];
        momentsEq[8] = -((T)2/(T)3)*j[2];
        momentsEq[9] = ((T)2*j[0]*j[0]-j[1]*j[1]-j[2]*j[2])*invRho;
        momentsEq[10] = (-j[0]*j[0]+(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2])*invRho;
        momentsEq[11] = (j[1]*j[1]-j[2]*j[2])*invRho;
        momentsEq[12] = (-(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2])*invRho;
        momentsEq[13] = j[1]*j[0]*invRho;
        momentsEq[14] = j[2]*j[1]*invRho;
        momentsEq[15] = j[2]*j[0]*invRho;
        momentsEq[16] = T();
        momentsEq[17] = T();
        momentsEq[18] = T();
        
    }
    
    /// Computation of all equilibrium distribution (in moments space)
    static void computeSmagorinskyEquilibriumMoments( Array<T,Descriptor::q>& momentsEq,
                                               T rhoBar, Array<T,3> const& j, T jSqr, const Array<T,6> &strain, T cSmago )
    {
        typedef SymmetricTensorImpl<T,3> S;
        
        T invRho = Descriptor::invRho(rhoBar);
        T rho = Descriptor::fullRho(rhoBar);
        T rho2 = rho*rho;
        
        T sNorm = std::sqrt((T)2*SymmetricTensorImpl<T,3>::tensorNormSqr(strain));
        T smagoFactor = (T)2*cSmago * cSmago * sNorm;
        
        T ux2 = rho2 * smagoFactor * strain[S::xx];
        T uy2 = rho2 * smagoFactor * strain[S::yy];
        T uz2 = rho2 * smagoFactor * strain[S::zz];
        
        T uxuy = rho2 * smagoFactor * strain[S::xy];
        T uyuz = rho2 * smagoFactor * strain[S::yz];
        T uxuz = rho2 * smagoFactor * strain[S::xz];
        
        momentsEq[0] = rhoBar;
        momentsEq[1] = (T)19*(jSqr + ux2+uy2+uz2)*invRho-(T)11*rhoBar;
        momentsEq[2] = -(T)5.5*(jSqr + ux2+uy2+uz2)*invRho + (T)3*rhoBar;
        momentsEq[3] = j[0];
        momentsEq[4] = -((T)2/(T)3)*j[0];
        momentsEq[5] = j[1];
        momentsEq[6] = -((T)2/3)*j[1];
        momentsEq[7] = j[2];
        momentsEq[8] = -((T)2/(T)3)*j[2];
        momentsEq[9] = ((T)2*(j[0]*j[0]+ux2)-(j[1]*j[1]+uy2)-(j[2]*j[2]+uz2))*invRho;
        momentsEq[10] = (-(j[0]*j[0]+ux2)+(T)0.5*(j[1]*j[1]+uy2)+(T)0.5*(j[2]*j[2]+uz2))*invRho;
        momentsEq[11] = (j[1]*j[1]+uy2-(j[2]*j[2]+uz2))*invRho;
        momentsEq[12] = (-(T)0.5*(j[1]*j[1]+uy2)+(T)0.5*(j[2]*j[2]+uz2))*invRho;
        momentsEq[13] = (j[1]*j[0]+uxuy)*invRho;
        momentsEq[14] = (j[2]*j[1]+uyuz)*invRho;
        momentsEq[15] = (j[2]*j[0]+uxuz)*invRho;
        momentsEq[16] = T();
        momentsEq[17] = T();
        momentsEq[18] = T();
    }
    
    /// Computation of all moments (specialized for d3q19)
    static void computeMoments(Array<T,Descriptor::q>& moments, const Array<T,Descriptor::q>& f)
    {
        moments[0] = f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8]+f[9]+f[10]+f[11]+f[12]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18];
        moments[1] =
            -(T)30*f[0]-(T)11*(f[1]+f[2]+f[3]+f[10]+f[11]+f[12])
            +(T)8*(f[4]+f[5]+f[6]+f[7]+f[8]+f[9]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18]);
        moments[2] = 12*f[0]-4*f[1]-4*f[2]-4*f[3]+f[4]+f[5]+f[6]+f[7]+f[8]+f[9]-4*f[10]-4*f[11]-4*f[12]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18];
        moments[3] = -f[1]-f[4]-f[5]-f[6]-f[7]+f[10]+f[13]+f[14]+f[15]+f[16];
        moments[4] = 4*f[1]-f[4]-f[5]-f[6]-f[7]-4*f[10]+f[13]+f[14]+f[15]+f[16];
        moments[5] = -f[2]-f[4]+f[5]-f[8]-f[9]+f[11]+f[13]-f[14]+f[17]+f[18];
        moments[6] = 4*f[2]-f[4]+f[5]-f[8]-f[9]-4*f[11]+f[13]-f[14]+f[17]+f[18];
        moments[7] = -f[3]-f[6]+f[7]-f[8]+f[9]+f[12]+f[15]-f[16]+f[17]-f[18];
        moments[8] = 4*f[3]-f[6]+f[7]-f[8]+f[9]-4*f[12]+f[15]-f[16]+f[17]-f[18];
        moments[9] = 2*f[1]-f[2]-f[3]+f[4]+f[5]+f[6]+f[7]-2*f[8]-2*f[9]+2*f[10]-f[11]-f[12]+f[13]+f[14]+f[15]+f[16]-2*f[17]-2*f[18];
        moments[10] = -4*f[1]+2*f[2]+2*f[3]+f[4]+f[5]+f[6]+f[7]-2*f[8]-2*f[9]-4*f[10]+2*f[11]+2*f[12]+f[13]+f[14]+f[15]+f[16]-2*f[17]-2*f[18];
        moments[11] = f[2]-f[3]+f[4]+f[5]-f[6]-f[7]+f[11]-f[12]+f[13]+f[14]-f[15]-f[16];
        moments[12] = -2*f[2]+2*f[3]+f[4]+f[5]-f[6]-f[7]-2*f[11]+2*f[12]+f[13]+f[14]-f[15]-f[16];
        moments[13] = f[4]-f[5]+f[13]-f[14];
        moments[14] = f[8]-f[9]+f[17]-f[18];
        moments[15] = f[6]-f[7]+f[15]-f[16];
        moments[16] = -f[4]-f[5]+f[6]+f[7]+f[13]+f[14]-f[15]-f[16];
        moments[17] = f[4]-f[5]-f[8]-f[9]-f[13]+f[14]+f[17]+f[18];
        moments[18] = -f[6]+f[7]+f[8]-f[9]+f[15]-f[16]-f[17]+f[18];
        
    }
    
    static void computef_InvM_Smoments(Array<T,19>& f, const Array<T,19> &moments, const T &omega) 
    {
        T mom0 = moments[0] * MRTDescriptor::S[0];
        T mom1 = moments[1] * MRTDescriptor::S[1];
        T mom2 = moments[2] * MRTDescriptor::S[2];
        T mom3 = moments[3] * MRTDescriptor::S[3];
        T mom4 = moments[4] * MRTDescriptor::S[4];
        T mom5 = moments[5] * MRTDescriptor::S[5];
        T mom6 = moments[6] * MRTDescriptor::S[6];
        T mom7 = moments[7] * MRTDescriptor::S[7];
        T mom8 = moments[8] * MRTDescriptor::S[8];
        T mom9 = moments[9] * omega;
        T mom10 = moments[10] * MRTDescriptor::S[10];
        T mom11 = moments[11] * omega;
        T mom12 = moments[12] * MRTDescriptor::S[12];
        T mom13 = moments[13] * omega;
        T mom14 = moments[14] * omega;
        T mom15 = moments[15] * omega;
        T mom16 = moments[16] * MRTDescriptor::S[16];
        T mom17 = moments[17] * MRTDescriptor::S[17];
        T mom18 = moments[18] * MRTDescriptor::S[18];
        
        T mom0tmp = mom0 / (T)19;
        
        f[0] -= mom0tmp-5/(T)399*mom1+1/(T)21*mom2;
        
        T mom1tmp = (T)11/(T)2394*mom1;
        T mom2tmp = mom2/(T)63;
        T mom3_m4 = (T)0.1*(mom3-mom4);
        T mom9_m10 = (mom9-mom10)/(T)18;
        f[1] -= mom0tmp-mom1tmp-mom2tmp-mom3_m4+mom9_m10;
        f[10] -= mom0tmp-mom1tmp-mom2tmp+mom3_m4+mom9_m10;
        
        T mom5_m6 = (T)0.1*(mom5-mom6);
        T mom7_m8 = (T)0.1*(mom7-mom8);
        mom9_m10 *= (T)0.5;
        T mom11_m12 = (mom11-mom12)/(T)12;
        f[2] -= mom0tmp-mom1tmp-mom2tmp-mom5_m6-mom9_m10+mom11_m12;
        f[3] -= mom0tmp-mom1tmp-mom2tmp-mom7_m8-mom9_m10-mom11_m12;
        f[11] -= mom0tmp-mom1tmp-mom2tmp+mom5_m6-mom9_m10+mom11_m12;
        f[12] -= mom0tmp-mom1tmp-mom2tmp+mom7_m8-mom9_m10-mom11_m12;
        
        mom1tmp = (T)4/(T)1197*mom1;
        mom2tmp *= (T)0.25;
        T mom5_p6 = (T)0.1*mom5+(T)0.025*mom6;
        T mom7_p8 = (T)0.1*mom7+(T)0.025*mom8;
        T mom9_p10 = (mom9+mom10*(T)0.5)/(T)18;
        mom14 *= (T)0.25;
        mom17 *= (T)0.125;
        mom18 *= (T)0.125;
        f[8] -= mom0tmp+mom1tmp+mom2tmp-mom5_p6-mom7_p8-mom9_p10+mom14-mom17+mom18;<commit_msg>made different relaxation time settings available in mrtTemplates3D.h<commit_after>/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011-2015 FlowKit Sarl
 * Route d'Oron 2
 * 1010 Lausanne, Switzerland
 * E-mail contact: <EMAIL>
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/* Orestis Malaspinas contributed this code.
 */

/** \file
 * Template specializations for some computationally intensive LB 
 * functions of the header file mrtTemplates.h, for the D3Q19 grid.
 */

#ifndef MRT_TEMPLATES_3D_H
#define MRT_TEMPLATES_3D_H

#include "core/globalDefs.h"

namespace plb {

// Efficient specialization for D3Q19 lattice
template<typename T> struct mrtTemplatesImpl<T, descriptors::MRTD3Q19DescriptorBase<T> > {
    
    typedef descriptors::D3Q19DescriptorBase<T> Descriptor; 
    typedef descriptors::MRTD3Q19DescriptorBase<T> MRTDescriptor; 

    /// Computation of all equilibrium distribution (in moments space)
    static void computeEquilibriumMoments( Array<T,Descriptor::q>& momentsEq,
                                    T rhoBar, Array<T,3> const& j, T jSqr )
    {
        T invRho = Descriptor::invRho(rhoBar);
        
        momentsEq[0] = rhoBar;
        momentsEq[1] = (T)19*jSqr*invRho-(T)11*rhoBar;
        momentsEq[2] = -(T)5.5*jSqr*invRho+(T)3*rhoBar;
        momentsEq[3] = j[0];
        momentsEq[4] = -((T)2/(T)3)*j[0];
        momentsEq[5] = j[1];
        momentsEq[6] = -((T)2/3)*j[1];
        momentsEq[7] = j[2];
        momentsEq[8] = -((T)2/(T)3)*j[2];
        momentsEq[9] = ((T)2*j[0]*j[0]-j[1]*j[1]-j[2]*j[2])*invRho;
        momentsEq[10] = (-j[0]*j[0]+(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2])*invRho;
        momentsEq[11] = (j[1]*j[1]-j[2]*j[2])*invRho;
        momentsEq[12] = (-(T)0.5*j[1]*j[1]+(T)0.5*j[2]*j[2])*invRho;
        momentsEq[13] = j[1]*j[0]*invRho;
        momentsEq[14] = j[2]*j[1]*invRho;
        momentsEq[15] = j[2]*j[0]*invRho;
        momentsEq[16] = T();
        momentsEq[17] = T();
        momentsEq[18] = T();
        
    }
    
    /// Computation of all equilibrium distribution (in moments space)
    static void computeSmagorinskyEquilibriumMoments( Array<T,Descriptor::q>& momentsEq,
                                               T rhoBar, Array<T,3> const& j, T jSqr, const Array<T,6> &strain, T cSmago )
    {
        typedef SymmetricTensorImpl<T,3> S;
        
        T invRho = Descriptor::invRho(rhoBar);
        T rho = Descriptor::fullRho(rhoBar);
        T rho2 = rho*rho;
        
        T sNorm = std::sqrt((T)2*SymmetricTensorImpl<T,3>::tensorNormSqr(strain));
        T smagoFactor = (T)2*cSmago * cSmago * sNorm;
        
        T ux2 = rho2 * smagoFactor * strain[S::xx];
        T uy2 = rho2 * smagoFactor * strain[S::yy];
        T uz2 = rho2 * smagoFactor * strain[S::zz];
        
        T uxuy = rho2 * smagoFactor * strain[S::xy];
        T uyuz = rho2 * smagoFactor * strain[S::yz];
        T uxuz = rho2 * smagoFactor * strain[S::xz];
        
        momentsEq[0] = rhoBar;
        momentsEq[1] = (T)19*(jSqr + ux2+uy2+uz2)*invRho-(T)11*rhoBar;
        momentsEq[2] = -(T)5.5*(jSqr + ux2+uy2+uz2)*invRho + (T)3*rhoBar;
        momentsEq[3] = j[0];
        momentsEq[4] = -((T)2/(T)3)*j[0];
        momentsEq[5] = j[1];
        momentsEq[6] = -((T)2/3)*j[1];
        momentsEq[7] = j[2];
        momentsEq[8] = -((T)2/(T)3)*j[2];
        momentsEq[9] = ((T)2*(j[0]*j[0]+ux2)-(j[1]*j[1]+uy2)-(j[2]*j[2]+uz2))*invRho;
        momentsEq[10] = (-(j[0]*j[0]+ux2)+(T)0.5*(j[1]*j[1]+uy2)+(T)0.5*(j[2]*j[2]+uz2))*invRho;
        momentsEq[11] = (j[1]*j[1]+uy2-(j[2]*j[2]+uz2))*invRho;
        momentsEq[12] = (-(T)0.5*(j[1]*j[1]+uy2)+(T)0.5*(j[2]*j[2]+uz2))*invRho;
        momentsEq[13] = (j[1]*j[0]+uxuy)*invRho;
        momentsEq[14] = (j[2]*j[1]+uyuz)*invRho;
        momentsEq[15] = (j[2]*j[0]+uxuz)*invRho;
        momentsEq[16] = T();
        momentsEq[17] = T();
        momentsEq[18] = T();
    }
    
    /// Computation of all moments (specialized for d3q19)
    static void computeMoments(Array<T,Descriptor::q>& moments, const Array<T,Descriptor::q>& f)
    {
        moments[0] = f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8]+f[9]+f[10]+f[11]+f[12]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18];
        moments[1] =
            -(T)30*f[0]-(T)11*(f[1]+f[2]+f[3]+f[10]+f[11]+f[12])
            +(T)8*(f[4]+f[5]+f[6]+f[7]+f[8]+f[9]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18]);
        moments[2] = 12*f[0]-4*f[1]-4*f[2]-4*f[3]+f[4]+f[5]+f[6]+f[7]+f[8]+f[9]-4*f[10]-4*f[11]-4*f[12]+f[13]+f[14]+f[15]+f[16]+f[17]+f[18];
        moments[3] = -f[1]-f[4]-f[5]-f[6]-f[7]+f[10]+f[13]+f[14]+f[15]+f[16];
        moments[4] = 4*f[1]-f[4]-f[5]-f[6]-f[7]-4*f[10]+f[13]+f[14]+f[15]+f[16];
        moments[5] = -f[2]-f[4]+f[5]-f[8]-f[9]+f[11]+f[13]-f[14]+f[17]+f[18];
        moments[6] = 4*f[2]-f[4]+f[5]-f[8]-f[9]-4*f[11]+f[13]-f[14]+f[17]+f[18];
        moments[7] = -f[3]-f[6]+f[7]-f[8]+f[9]+f[12]+f[15]-f[16]+f[17]-f[18];
        moments[8] = 4*f[3]-f[6]+f[7]-f[8]+f[9]-4*f[12]+f[15]-f[16]+f[17]-f[18];
        moments[9] = 2*f[1]-f[2]-f[3]+f[4]+f[5]+f[6]+f[7]-2*f[8]-2*f[9]+2*f[10]-f[11]-f[12]+f[13]+f[14]+f[15]+f[16]-2*f[17]-2*f[18];
        moments[10] = -4*f[1]+2*f[2]+2*f[3]+f[4]+f[5]+f[6]+f[7]-2*f[8]-2*f[9]-4*f[10]+2*f[11]+2*f[12]+f[13]+f[14]+f[15]+f[16]-2*f[17]-2*f[18];
        moments[11] = f[2]-f[3]+f[4]+f[5]-f[6]-f[7]+f[11]-f[12]+f[13]+f[14]-f[15]-f[16];
        moments[12] = -2*f[2]+2*f[3]+f[4]+f[5]-f[6]-f[7]-2*f[11]+2*f[12]+f[13]+f[14]-f[15]-f[16];
        moments[13] = f[4]-f[5]+f[13]-f[14];
        moments[14] = f[8]-f[9]+f[17]-f[18];
        moments[15] = f[6]-f[7]+f[15]-f[16];
        moments[16] = -f[4]-f[5]+f[6]+f[7]+f[13]+f[14]-f[15]-f[16];
        moments[17] = f[4]-f[5]-f[8]-f[9]-f[13]+f[14]+f[17]+f[18];
        moments[18] = -f[6]+f[7]+f[8]-f[9]+f[15]-f[16]-f[17]+f[18];
        
    }
    
    static void computef_InvM_Smoments(Array<T,19>& f, const Array<T,19> &moments, const T &omega) 
    {

      // TRT relaxation rates
#ifdef MRT_USE_TRT_RELAXATION

        T omega_min = 8.*(2.-omega)/(8.-omega);

        T mom0 = moments[0];
        T mom1 = moments[1] * omega;
        T mom2 = moments[2] * omega;
        T mom3 = moments[3];
        T mom4 = moments[4] * omega_min;
        T mom5 = moments[5];
        T mom6 = moments[6] * omega_min;
        T mom7 = moments[7];
        T mom8 = moments[8] * omega_min;
        T mom9 = moments[9] * omega;
        T mom10 = moments[10] * omega;
        T mom11 = moments[11] * omega;
        T mom12 = moments[12] * omega;
        T mom13 = moments[13] * omega;
        T mom14 = moments[14] * omega;
        T mom15 = moments[15] * omega;
        T mom16 = moments[16] * omega_min;
        T mom17 = moments[17] * omega_min;
        T mom18 = moments[18] * omega_min;

        // relaxation rates from Pan et al. (2006)
#elif defined MRT_USE_PAN_RELAXATION

        T omega_min = 8.*(2.-omega)/(8.-omega);

        T mom0 = moments[0];
        T mom1 = moments[1] * omega_min;
        T mom2 = moments[2] * omega_min;
        T mom3 = moments[3];
        T mom4 = moments[4] * omega_min;
        T mom5 = moments[5];
        T mom6 = moments[6] * omega_min;
        T mom7 = moments[7];
        T mom8 = moments[8] * omega_min;
        T mom9 = moments[9] * omega;
        T mom10 = moments[10] * omega_min;
        T mom11 = moments[11] * omega;
        T mom12 = moments[12] * omega_min;
        T mom13 = moments[13] * omega;
        T mom14 = moments[14] * omega;
        T mom15 = moments[15] * omega;
        T mom16 = moments[16] * omega_min;
        T mom17 = moments[17] * omega_min;
        T mom18 = moments[18] * omega_min;

        // relaxation rates proposed by Palabos authors (and probably d'Humieres et al)
#else

        T mom0 = moments[0] * MRTDescriptor::S[0];
        T mom1 = moments[1] * MRTDescriptor::S[1];
        T mom2 = moments[2] * MRTDescriptor::S[2];
        T mom3 = moments[3] * MRTDescriptor::S[3];
        T mom4 = moments[4] * MRTDescriptor::S[4];
        T mom5 = moments[5] * MRTDescriptor::S[5];
        T mom6 = moments[6] * MRTDescriptor::S[6];
        T mom7 = moments[7] * MRTDescriptor::S[7];
        T mom8 = moments[8] * MRTDescriptor::S[8];
        T mom9 = moments[9] * omega;
        T mom10 = moments[10] * MRTDescriptor::S[10];
        T mom11 = moments[11] * omega;
        T mom12 = moments[12] * MRTDescriptor::S[12];
        T mom13 = moments[13] * omega;
        T mom14 = moments[14] * omega;
        T mom15 = moments[15] * omega;
        T mom16 = moments[16] * MRTDescriptor::S[16];
        T mom17 = moments[17] * MRTDescriptor::S[17];
        T mom18 = moments[18] * MRTDescriptor::S[18];

#endif
        
        T mom0tmp = mom0 / (T)19;
        
        f[0] -= mom0tmp-5/(T)399*mom1+1/(T)21*mom2;
        
        T mom1tmp = (T)11/(T)2394*mom1;
        T mom2tmp = mom2/(T)63;
        T mom3_m4 = (T)0.1*(mom3-mom4);
        T mom9_m10 = (mom9-mom10)/(T)18;
        f[1] -= mom0tmp-mom1tmp-mom2tmp-mom3_m4+mom9_m10;
        f[10] -= mom0tmp-mom1tmp-mom2tmp+mom3_m4+mom9_m10;
        
        T mom5_m6 = (T)0.1*(mom5-mom6);
        T mom7_m8 = (T)0.1*(mom7-mom8);
        mom9_m10 *= (T)0.5;
        T mom11_m12 = (mom11-mom12)/(T)12;
        f[2] -= mom0tmp-mom1tmp-mom2tmp-mom5_m6-mom9_m10+mom11_m12;
        f[3] -= mom0tmp-mom1tmp-mom2tmp-mom7_m8-mom9_m10-mom11_m12;
        f[11] -= mom0tmp-mom1tmp-mom2tmp+mom5_m6-mom9_m10+mom11_m12;
        f[12] -= mom0tmp-mom1tmp-mom2tmp+mom7_m8-mom9_m10-mom11_m12;
        
        mom1tmp = (T)4/(T)1197*mom1;
        mom2tmp *= (T)0.25;
        T mom5_p6 = (T)0.1*mom5+(T)0.025*mom6;
        T mom7_p8 = (T)0.1*mom7+(T)0.025*mom8;
        T mom9_p10 = (mom9+mom10*(T)0.5)/(T)18;
        mom14 *= (T)0.25;
        mom17 *= (T)0.125;
        mom18 *= (T)0.125;
        f[8] -= mom0tmp+mom1tmp+mom2tmp-mom5_p6-mom7_p8-mom9_p10+mom14-mom17+mom18;
{'fraction_non_alphanumeric': 0.21453053593685084, 'fraction_numerical': 0.10100747818861654, 'mean_word_length': 2.775142128994315, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6799627', 'n_tokens_mistral': 9347, 'n_tokens_neox': 8516, 'n_words': 1307}