<commit_before>#pragma once
#include <exception>
#include "dg/algorithm.h"
#include "parameters.h"
#include "dg/polarization.h"

namespace poet
{


template<class Geometry, class Matrix, class container>
struct Implicit
{
    Implicit( const Geometry& g, double nu):
        nu_(nu),     
        temp( dg::evaluate(dg::zero, g)),
        LaplacianM_perp( g, dg::normed, dg::centered){ }
    void operator()(double t, const std::vector<container>& x, std::vector<container>& y)
    {
        /* x[0] := N_e - 1
         * x[2] := N_i - 1 
         */
        for( unsigned i=0; i<x.size(); i++)
        {
//             dg::blas2::gemv( LaplacianM_perp, x[i], temp);
//             dg::blas2::gemv( LaplacianM_perp, temp, y[i]);
//             dg::blas1::scal( y[i], -nu_);
            dg::blas2::gemv( -nu_, LaplacianM_perp, x[i], 0., y[i]);
        }
    }
    dg::Elliptic<Geometry, Matrix, container>& laplacianM() {return LaplacianM_perp;}
    const container& weights(){return LaplacianM_perp.weights();}
    const container& inv_weights(){return LaplacianM_perp.inv_weights();}
    const container& precond(){return LaplacianM_perp.precond();}

  private:
    double nu_;
    container temp;
    dg::Elliptic<Geometry, Matrix, container> LaplacianM_perp;
};

template< class Geometry, class Matrix, class DiaMatrix, class CooMatrix, class container >
struct Explicit
{
    /**
     * @brief Construct a Explicit solver object
     *
     * @param g The grid on which to operate
     * @param p The parameters
     */
    Explicit( const Geometry& g, const Parameters& p );


    /**
     * @brief Returns phi and psi that belong to the last y in operator()
     *
     * In a multistep scheme this belongs to the point HEAD-1
     * @return phi[0] is the electron and phi[1] the generalized ion potential
     */
    const std::vector<container>& potential( ) const { return phi;}

    /**
     * @brief Return the normalized negative laplacian used by this object
     *
     * @return cusp matrix
     */
    dg::Elliptic<Geometry, Matrix, container>& laplacianM( ) { return laplaceM;}

    void gamma1_y( const container& y, container& yp)
    {
        if (equations == "global-arbpolO2") {
            dg::blas1::scal(yp,0.0);
            sqrtinvert(yp, y);
        }
        else {
            std::vector<unsigned> number = multigrid.direct_solve( multi_gamma1, yp, y, eps_gamma);
            if(  number[0] == multigrid.max_iter())
                throw dg::Fail( eps_gamma);
        }
    }
    void gamma1inv_y( const container& y, container& yp)
    {
        if (equations == "global-arbpolO2")
        {
            sqrtsolve(y, yp);
        }
        else
        {
        dg::blas2::symv( multi_gamma1[0],y,chi); //invG ne-1
        dg::blas2::symv( v2d, chi, yp);
        }
    }
    /**
     * @brief Compute the right-hand side of the poet equations
     *
     * y[0] = N_e - 1, 
     * y[1] = N_i - 1 || y[1] = Omega
     * @param y input vector
     * @param yp the rhs yp = f(y)
     */
    void operator()( double t, const std::vector<container>& y, std::vector<container>& yp);

    /**
     * @brief Return the mass of the last field in operator() in a global computation
     *
     * @return int exp(y[0]) dA
     * @note undefined for a local computation
     */
    double mass( ) {return mass_;}
    /**
     * @brief Return the last integrated mass diffusion of operator() in a global computation
     *
     * @return int \nu \Delta (exp(y[0])-1)
     * @note undefined for a local computation
     */
    double mass_diffusion( ) {return diff_;}
    /**
     * @brief Return the energy of the last field in operator() in a global computation
     *
     * @return integrated total energy in {ne, ni}
     * @note undefined for a local computation
     */
    double energy( ) {return energy_;}
    /**
     * @brief Return the integrated energy diffusion of the last field in operator() in a global computation
     *
     * @return integrated total energy diffusion
     * @note undefined for a local computation
     */
    double energy_diffusion( ){ return ediff_;}

  private:
    //use chi and omega as helpers to compute square velocity in omega
    const container& compute_psi( double t, const container& potential);
    const container& polarisation( double t, const std::vector<container>& y);

    container chi, omega, iota, gamma_n, gamma_phi;
    const container binv; //magnetic field
    container n_old, gamma_n_old, phi_old, gamma_phi_old;
    std::vector<container> phi, dyphi, ype;
    std::vector<container> dyy, lny, lapy;
    
    
    //matrices and solvers
    dg::Elliptic<Geometry, Matrix, container>  laplaceM; //contains normalized laplacian
    std::vector<dg::Elliptic<Geometry, Matrix, container> > multi_elliptic;
    std::vector<dg::TensorElliptic<Geometry, Matrix, container> > multi_tensorelliptic;
    std::vector<dg::Helmholtz<Geometry,  Matrix, container> > multi_gamma1, multi_gamma0;
    
    dg::KrylovSqrtCauchySolve< Geometry, Matrix, DiaMatrix, CooMatrix, container, dg::DVec> sqrtsolve;
    dg::KrylovSqrtCauchyinvert<Geometry, Matrix, DiaMatrix, CooMatrix, container, dg::DVec> sqrtinvert;
    
    dg::ArakawaX< Geometry, Matrix, container> arakawa;
    dg::Advection<Geometry, Matrix, container> m_adv;
    
    dg::MultigridCG2d<Geometry, Matrix, container> multigrid;
    dg::Extrapolation<container> old_phi, old_gamma_phi, old_psi, old_gamma_n;
    std::vector<container> multi_chi, multi_iota;
    
    Matrix m_centered[2];
 
    
    const container w2d,v2d, one;
    const std::vector<double> eps_pol;
    const double eps_gamma;
    const double kappa,  nu, tau;
    const std::string equations;
    bool boussinesq;

    double mass_, energy_, diff_, ediff_;


};

template< class Geometry, class M, class DM, class CM, class container>
Explicit< Geometry, M, DM, CM, container>::Explicit( const Geometry& grid, const Parameters& p ):
    chi( evaluate( dg::zero, grid)), omega(chi), iota(chi), gamma_n(chi), gamma_phi(chi),
    binv( evaluate( dg::LinearX( p.kappa, 1.-p.kappa*p.posX*p.lx), grid)),
    n_old(chi),    gamma_n_old(chi), phi_old(chi),    gamma_phi_old(chi),
    phi( 2, chi), dyphi( phi), ype(phi),
    dyy(2,chi), lny( dyy), lapy(dyy), 
    laplaceM( grid, dg::normed, dg::centered),
    multigrid( grid, 3),
    old_phi( 2, chi),  old_gamma_phi(2, chi), old_psi( 2, chi), old_gamma_n( 2, chi),
    w2d( dg::create::volume(grid)), v2d( dg::create::inv_weights(grid)), one( dg::evaluate(dg::one, grid)),
    eps_pol(p.eps_pol), eps_gamma( p.eps_gamma), kappa(p.kappa), nu(p.nu), tau( p.tau), equations( p.equations),  boussinesq(p.boussinesq)
{
    multi_chi= multigrid.project( chi);
    multi_iota= multigrid.project( chi);
    multi_elliptic.resize(3);
    multi_tensorelliptic.resize(3);
    multi_dfpolcharge.resize(3);
    multi_gamma1.resize(3);
    multi_gamma0.resize(3);
    arakawa.construct(grid);
    m_adv.construct(grid);
    m_centered[0] = dg::create::dx( grid, grid.bcx(), dg::centered);
    m_centered[1] = dg::create::dy( grid, grid.bcy(), dg::centered);
    for( unsigned u=0; u<3; u++)
    {
        multi_elliptic[u].construct( multigrid.grid(u), dg::not_normed, dg::centered, p.jfactor);
        multi_tensorelliptic[u].construct( multigrid.grid(u),  dg::not_normed, dg::centered, p.jfactor); //only centered implemented

        
        multi_gamma0[u].construct( multigrid.grid(u), -p.tau, dg::centered, p.jfactor);
        if( equations == "global-arbpolO2" ) {
            multi_gamma1[u].construct( multigrid.grid(u), -p.tau, dg::centered, p.jfactor);
            sqrtsolve.construct( multi_gamma1[0], multigrid.grid(0), chi,  p.eps_time, 200, 20,  p.eps_gamma);
            sqrtinvert.construct(multi_gamma1[0], multigrid.grid(0), chi,  p.eps_time, 200, 20,  p.eps_gamma);
        }
        else {<commit_msg>cleaned poet code<commit_after>#pragma once
#include <exception>
#include "dg/algorithm.h"
#include "parameters.h"
#include "dg/polarization.h"

namespace poet
{


template<class Geometry, class Matrix, class container>
struct Implicit
{
    Implicit( const Geometry& g, double nu):
        nu_(nu),     
        temp( dg::evaluate(dg::zero, g)),
        LaplacianM_perp( g, dg::normed, dg::centered){ }
    void operator()(double t, const std::vector<container>& x, std::vector<container>& y)
    {
        /* x[0] := N_e - 1
         * x[2] := N_i - 1 
         */
        for( unsigned i=0; i<x.size(); i++)
        {
//             dg::blas2::gemv( LaplacianM_perp, x[i], temp);
//             dg::blas2::gemv( LaplacianM_perp, temp, y[i]);
//             dg::blas1::scal( y[i], -nu_);
            dg::blas2::gemv( -nu_, LaplacianM_perp, x[i], 0., y[i]);
        }
    }
    dg::Elliptic<Geometry, Matrix, container>& laplacianM() {return LaplacianM_perp;}
    const container& weights(){return LaplacianM_perp.weights();}
    const container& inv_weights(){return LaplacianM_perp.inv_weights();}
    const container& precond(){return LaplacianM_perp.precond();}

  private:
    double nu_;
    container temp;
    dg::Elliptic<Geometry, Matrix, container> LaplacianM_perp;
};

template< class Geometry, class Matrix, class DiaMatrix, class CooMatrix, class container >
struct Explicit
{
    /**
     * @brief Construct a Explicit solver object
     *
     * @param g The grid on which to operate
     * @param p The parameters
     */
    Explicit( const Geometry& g, const Parameters& p );


    /**
     * @brief Returns phi and psi that belong to the last y in operator()
     *
     * In a multistep scheme this belongs to the point HEAD-1
     * @return phi[0] is the electron and phi[1] the generalized ion potential
     */
    const std::vector<container>& potential( ) const { return phi;}

    /**
     * @brief Return the normalized negative laplacian used by this object
     *
     * @return cusp matrix
     */
    dg::Elliptic<Geometry, Matrix, container>& laplacianM( ) { return laplaceM;}

    void gamma1_y( const container& y, container& yp)
    {
        if (equations == "global-arbpolO2") {
            dg::blas1::scal(yp,0.0);
            sqrtinvert(yp, y);
        }
        else {
            std::vector<unsigned> number = multigrid.direct_solve( multi_gamma1, yp, y, eps_gamma);
            if(  number[0] == multigrid.max_iter())
                throw dg::Fail( eps_gamma);
        }
    }
    void gamma1inv_y( const container& y, container& yp)
    {
        if (equations == "global-arbpolO2")
        {
            sqrtsolve(y, yp);
        }
        else
        {
        dg::blas2::symv( multi_gamma1[0],y,chi); //invG ne-1
        dg::blas2::symv( v2d, chi, yp);
        }
    }
    /**
     * @brief Compute the right-hand side of the poet equations
     *
     * y[0] = N_e - 1, 
     * y[1] = N_i - 1 || y[1] = Omega
     * @param y input vector
     * @param yp the rhs yp = f(y)
     */
    void operator()( double t, const std::vector<container>& y, std::vector<container>& yp);

    /**
     * @brief Return the mass of the last field in operator() in a global computation
     *
     * @return int exp(y[0]) dA
     * @note undefined for a local computation
     */
    double mass( ) {return mass_;}
    /**
     * @brief Return the last integrated mass diffusion of operator() in a global computation
     *
     * @return int \nu \Delta (exp(y[0])-1)
     * @note undefined for a local computation
     */
    double mass_diffusion( ) {return diff_;}
    /**
     * @brief Return the energy of the last field in operator() in a global computation
     *
     * @return integrated total energy in {ne, ni}
     * @note undefined for a local computation
     */
    double energy( ) {return energy_;}
    /**
     * @brief Return the integrated energy diffusion of the last field in operator() in a global computation
     *
     * @return integrated total energy diffusion
     * @note undefined for a local computation
     */
    double energy_diffusion( ){ return ediff_;}

  private:
    //use chi and omega as helpers to compute square velocity in omega
    const container& compute_psi( double t, const container& potential);
    const container& polarisation( double t, const std::vector<container>& y);

    container chi, omega, iota, gamma_n, gamma_phi;
    const container binv; //magnetic field
    container n_old, gamma_n_old, phi_old, gamma_phi_old;
    std::vector<container> phi, dyphi, ype;
    std::vector<container> dyy, lny, lapy;
    
    
    //matrices and solvers
    dg::Elliptic<Geometry, Matrix, container>  laplaceM; //contains normalized laplacian
    std::vector<dg::Elliptic<Geometry, Matrix, container> > multi_elliptic;
    std::vector<dg::TensorElliptic<Geometry, Matrix, container> > multi_tensorelliptic;
    std::vector<dg::Helmholtz<Geometry,  Matrix, container> > multi_gamma1, multi_gamma0;
    
    dg::KrylovSqrtCauchySolve< Geometry, Matrix, DiaMatrix, CooMatrix, container, dg::DVec> sqrtsolve;
    dg::KrylovSqrtCauchyinvert<Geometry, Matrix, DiaMatrix, CooMatrix, container, dg::DVec> sqrtinvert;
    
    dg::ArakawaX< Geometry, Matrix, container> arakawa;
    dg::Advection<Geometry, Matrix, container> m_adv;
    
    dg::MultigridCG2d<Geometry, Matrix, container> multigrid;
    dg::Extrapolation<container> old_phi, old_gamma_phi, old_psi, old_gamma_n;
    std::vector<container> multi_chi, multi_iota;
    
    Matrix m_centered[2];
 
    
    const container w2d,v2d, one;
    const std::vector<double> eps_pol;
    const double eps_gamma;
    const double kappa,  nu, tau;
    const std::string equations;
    bool boussinesq;

    double mass_, energy_, diff_, ediff_;


};

template< class Geometry, class M, class DM, class CM, class container>
Explicit< Geometry, M, DM, CM, container>::Explicit( const Geometry& grid, const Parameters& p ):
    chi( evaluate( dg::zero, grid)), omega(chi), iota(chi), gamma_n(chi), gamma_phi(chi),
    binv( evaluate( dg::LinearX( p.kappa, 1.-p.kappa*p.posX*p.lx), grid)),
    n_old(chi),    gamma_n_old(chi), phi_old(chi),    gamma_phi_old(chi),
    phi( 2, chi), dyphi( phi), ype(phi),
    dyy(2,chi), lny( dyy), lapy(dyy), 
    laplaceM( grid, dg::normed, dg::centered),
    multigrid( grid, 3),
    old_phi( 2, chi),  old_gamma_phi(2, chi), old_psi( 2, chi), old_gamma_n( 2, chi),
    w2d( dg::create::volume(grid)), v2d( dg::create::inv_weights(grid)), one( dg::evaluate(dg::one, grid)),
    eps_pol(p.eps_pol), eps_gamma( p.eps_gamma), kappa(p.kappa), nu(p.nu), tau( p.tau), equations( p.equations),  boussinesq(p.boussinesq)
{
    multi_chi= multigrid.project( chi);
    multi_iota= multigrid.project( chi);
    multi_elliptic.resize(3);
    multi_tensorelliptic.resize(3);
    multi_gamma1.resize(3);
    multi_gamma0.resize(3);
    arakawa.construct(grid);
    m_adv.construct(grid);
    m_centered[0] = dg::create::dx( grid, grid.bcx(), dg::centered);
    m_centered[1] = dg::create::dy( grid, grid.bcy(), dg::centered);
    for( unsigned u=0; u<3; u++)
    {
        multi_elliptic[u].construct( multigrid.grid(u), dg::not_normed, dg::centered, p.jfactor);
        multi_tensorelliptic[u].construct( multigrid.grid(u),  dg::not_normed, dg::centered, p.jfactor); //only centered implemented

        
        multi_gamma0[u].construct( multigrid.grid(u), -p.tau, dg::centered, p.jfactor);
        if( equations == "global-arbpolO2" ) {
            multi_gamma1[u].construct( multigrid.grid(u), -p.tau, dg::centered, p.jfactor);
            sqrtsolve.construct( multi_gamma1[0], multigrid.grid(0), chi,  p.eps_time, 200, 20,  p.eps_gamma);
            sqrtinvert.construct(multi_gamma1[0], multigrid.grid(0), chi,  p.eps_time, 200, 20,  p.eps_gamma);
        }
        else {
{'fraction_non_alphanumeric': 0.13630857507172459, 'fraction_numerical': 0.01039209435766656, 'mean_word_length': 3.071113418115754, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5748715', 'n_tokens_mistral': 5566, 'n_tokens_neox': 5243, 'n_words': 1383}