text:
/*
 *  iaf_cond_exp_sfa_rr.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef IAF_COND_EXP_SFA_RR_H
#define IAF_COND_EXP_SFA_RR_H

// Generated includes:
#include "config.h"

#ifdef HAVE_GSL

// External includes:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"

/* BeginDocumentation
Name: iaf_cond_exp_sfa_rr - Simple conductance based leaky integrate-and-fire
                            neuron model.

Description:
iaf_cond_exp_sfa_rr is an iaf_cond_exp_sfa_rr i.e. an implementation of a
spiking neuron using IAF dynamics with conductance-based synapses,
with additional spike-frequency adaptation and relative refractory
mechanisms as described in Dayan+Abbott, 2001, page 166.

As for the iaf_cond_exp_sfa_rr, Incoming spike events induce a post-synaptic
change  of  conductance  modelled  by an  exponential  function.  The
exponential function is  normalised such that an event  of weight 1.0
results in a peak current of 1 nS.

Outgoing spike events induce a change of the adaptation and relative
refractory conductances by q_sfa and q_rr, respectively.  Otherwise
these conductances decay exponentially with time constants tau_sfa
and tau_rr, respectively.

Parameters:
The following parameters can be set in the status dictionary.

V_m        double - Membrane potential in mV
E_L        double - Leak reversal potential in mV.
C_m        double - Capacity of the membrane in pF
t_ref      double - Duration of refractory period in ms.
V_th       double - Spike threshold in mV.
V_reset    double - Reset potential of the membrane in mV.
E_ex       double - Excitatory reversal potential in mV.
E_in       double - Inhibitory reversal potential in mV.
g_L        double - Leak conductance in nS;
tau_syn_ex double - Time constant of the excitatory synaptic exponential
                    function in ms.
tau_syn_in double - Time constant of the inhibitory synaptic exponential
                    function in ms.
q_sfa      double - Outgoing spike activated quantal spike-frequency adaptation
                    conductance increase in nS.
q_rr       double - Outgoing spike activated quantal relative refractory
                    conductance increase in nS.
tau_sfa    double - Time constant of spike-frequency adaptation in ms.
tau_rr     double - Time constant of the relative refractory mechanism in ms.
E_sfa      double - spike-frequency adaptation conductance reversal potential in
                    mV.
E_rr       double - relative refractory mechanism conductance reversal potential
                    in mV.
I_e        double - an external stimulus current in pA.

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

References:

Meffin, H., Burkitt, A. N., & Grayden, D. B. (2004). An analytical
model for the large, fluctuating synaptic conductance state typical of
neocortical neurons in vivo. J.  Comput. Neurosci., 16, 159-175.

Dayan, P. and Abbott, L. F. (2001). Theoretical Neuroscience, MIT Press (p166)

Author: Sven Schrader, Eilif Muller

SeeAlso: iaf_cond_exp_sfa_rr, aeif_cond_alpha, iaf_psc_delta, iaf_psc_exp,
iaf_cond_alpha
*/

namespace nest
{
/**
 * Function computing right-hand side of ODE for GSL solver.
 * @note Must be declared here so we can befriend it in class.
 * @note Must have C-linkage for passing to GSL. Internally, it is
 *       a first-class C++ function, but cannot be a member function
 *       because of the C-linkage.
 * @note No point in declaring it inline, since it is called
 *       through a function pointer.
 * @param void* Pointer to model neuron instance.
 */
extern "C" int
iaf_cond_exp_sfa_rr_dynamics( double, const double*, double*, void* );

class iaf_cond_exp_sfa_rr : public Archiving_Node
{

public:
  iaf_cond_exp_sfa_rr();
  iaf_cond_exp_sfa_rr( const iaf_cond_exp_sfa_rr& );
  ~iaf_cond_exp_sfa_rr();

  /**
   * Import sets of overloaded virtual functions.
   * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
   * Hiding
   */
  using Node::handle;
  using Node::handles_test_event;

  port send_test_event( Node&, rport, synindex, bool );

  void handle( SpikeEvent& );
  void handle( CurrentEvent& );
  void handle( DataLoggingRequest& );

  port handles_test_event( SpikeEvent&, rport );
  port handles_test_event( CurrentEvent&, rport );
  port handles_test_event( DataLoggingRequest&, rport );

  void get_status( DictionaryDatum& ) const;
  void set_status( const DictionaryDatum& );

private:
  void init_state_( const Node& proto );
  void init_buffers_();
  void calibrate();
  void update( Time const&, const long, const long );

  // END Boilerplate function declarations ----------------------------

  // Friends --------------------------------------------------------

  // make dynamics function quasi-member
  friend int
  iaf_cond_exp_sfa_rr_dynamics( double, const double*, double*, void* );

  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< iaf_cond_exp_sfa_rr >;
  friend class UniversalDataLogger< iaf_cond_exp_sfa_rr >;

private:
  // ----------------------------------------------------------------

  //! Independent parameters
  struct Parameters_
  {
    double V_th_;    //!< Threshold Potential in mV
    double V_reset_; //!< Reset Potential in mV
    double t_ref_;   //!< Refractory period in ms
    double g_L;      //!< Leak Conductance in nS
    double C_m;      //!< Membrane Capacitance in pF
    double E_ex;     //!< Excitatory reversal Potential in mV
    double E_in;     //!< Inhibitory reversal Potential in mV
    double E_L;      //!< Leak reversal Potential (aka resting potential) in mV
    double tau_synE; //!< Synaptic Time Constant Excitatory Synapse in ms
    double tau_synI; //!< Synaptic Time Constant for Inhibitory Synapse in ms
    double I_e;      //!< Constant Current in pA
    double tau_sfa;  //!< spike-frequency adaptation (sfa) time constant
    double tau_rr;   //!< relative refractory (rr) time constant
    double E_sfa;    //!< spike-frequency adaptation (sfa) reversal Potential
                     //!< in mV
    double E_rr;     //!<  relative refractory (rr) reversal Potential in mV
    double q_sfa;    //!< spike-frequency adaptation (sfa) quantal conductance
                     //!< increase in nS
    double q_rr;     //!< relative refractory (rr) quantal conductance increase
                     //!< in nS

    Parameters_(); //!< Sets default parameter values

    void get( DictionaryDatum& ) const; //!< Store current values in dictionary
    void set( const DictionaryDatum& ); //!< Set values from dicitonary
  };

public:
  // ----------------------------------------------------------------

  /**
   * State variables of the model.
   * @note Copy constructor and assignment operator required because
   *       of C-style array.
   */
  struct State_
  {

    //! Symbolic indices to the elements of the state vector y
    enum StateVecElems
    {
      V_M = 0,
      G_EXC,
      G_INH,
      G_SFA,
      G_RR,
      STATE_VEC_SIZE
    };

    //! neuron state, must be C-array for GSL solver
    double y_[ STATE_VEC_SIZE ];
    int r_; //!< number of refractory steps remaining

    State_( const Parameters_& ); //!< Default initialization
    State_( const State_& );
    State_& operator=( const State_& );

    void get( DictionaryDatum& ) const;
    void set( const DictionaryDatum&, const Parameters_& );
  };

private:
  // ----------------------------------------------------------------

  /**
   * Buffers of the model.
   */
  struct Buffers_
  {
    Buffers_( iaf_cond_exp_sfa_rr& ); //!<Sets buffer pointers to 0
    //! Sets buffer pointers to 0
    Buffers_( const Buffers_&, iaf_cond_exp_sfa_rr& );

    //! Logger for all analog data
    UniversalDataLogger< iaf_cond_exp_sfa_rr > logger_;

    /** buffers and sums up incoming spikes/currents */
    RingBuffer spike_exc_;
    RingBuffer spike_inh_;
    RingBuffer currents_;

    /** GSL ODE stuff */
    gsl_odeiv_step* s_;    //!< stepping function
    gsl_odeiv_control* c_; //!< adaptive stepsize control function
    gsl_odeiv_evolve* e_;  //!< evolution function
    gsl_odeiv_system sys_; //!< struct describing system

    // IntergrationStep_ should be reset with the neuron on ResetNetwork,
    // but remain unchanged during calibration. Since it is initialized with
    // step_, and the resolution cannot change after nodes have been created,
    // it is safe to place both here.
    double step_;            //!< step size in ms
    double IntegrationStep_; //!< current integration time step, updated by GSL

    /**
     * Input current injected by CurrentEvent.
     * This variable is used to transport the current applied into the
     * _dynamics function computing the derivative of the state vector.
     * It must be a part of Buffers_, since it is initialized once before
     * the first simulation, but not modified before later Simulate calls.
     */
    double I_stim_;
  };

  // ----------------------------------------------------------------

  /**
   * Internal variables of the model.
   */
  struct Variables_
  {
    int RefractoryCounts_;
  };

  // Access functions for UniversalDataLogger -------------------------------

  //! Read out state vector elements, used by UniversalDataLogger
  template < State_::StateVecElems elem >
  double
  get_y_elem_() const
  {
    return S_.y_[ elem ];
  }

  // ----------------------------------------------------------------

  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;

  //! Mapping of recordables names to access functions
  static RecordablesMap< iaf_cond_exp_sfa_rr > recordablesMap_;
};

inline port
nest::iaf_cond_exp_sfa_rr::send_test_event( Node& target,
  rport receptor_type,
  synindex,
  bool )
{
  SpikeEvent e;
  e.set_sender( *this );
  return target.handles_test_event( e, receptor_type );
}

inline port
iaf_cond_exp_sfa_rr::handles_test_event( SpikeEvent&, rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}

inline port
iaf_cond_exp_sfa_rr::handles_test_event( CurrentEvent&, rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}

inline port
iaf_cond_exp_sfa_rr::handles_test_event( DataLoggingRequest& dlr,
  rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}

inline void
iaf_cond_exp_sfa_rr::get_status( DictionaryDatum& d ) const
{
  P_.get( d );
  S_.get( d );
  Archiving_Node::get_status( d );

  ( *d )[ names::recordables ] = recordablesMap_.get_list();
}

inline void
iaf_cond_exp_sfa_rr::set_status( const DictionaryDatum& d )
{
  Parameters_ ptmp = P_; // temporary copy in case of errors
  ptmp.set( d );         // throws if BadProperty
  State_ stmp = S_;      // temporary copy in case of errors
  stmp.set( d, ptmp );   // throws if BadProperty

  // We now know that (ptmp, stmp) are consistent. We do not
  // write them back to (P_, S_) before we are also sure that
  // the properties to be set in the parent class are internally
  // consistent.
  Archiving_Node::set_status( d );

  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}

} // namespace

#endif // HAVE_GSL
#endif // IAF_COND_EXP_SFA_RR_H

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


text:
#ifndef S3D_VIDEO_FILE_PARSER_FFMPEG_SEEKER_H
#define S3D_VIDEO_FILE_PARSER_FFMPEG_SEEKER_H

#include <chrono>

#include <gsl/gsl>

struct AVFormatContext;

namespace s3d {

class Seeker {
 public:
  Seeker(gsl::not_null<AVFormatContext*> formatContext, int streamIndex);

  void seekTo(std::chrono::microseconds timestamp);

 private:
  AVFormatContext* formatContext_;
  int streamIndex_;
};

}  // namespace s3d

#endif  // S3D_VIDEO_FILE_PARSER_FFMPEG_SEEKER_H

ext: h


text:
#pragma once

#include <cstddef>
#include <memory>
#include <iostream>
#include <array>
#include <gsl\gsl>
#include "../utils/CRSP.h"
#include "./Pool.h"
#include "../math/OtherMath.h"
#include "../types/def.h"

namespace config {
    /**
    Default pool size for our Small Memory Allocator pools.
    */
    constexpr size_t sma_pool_size = 1024 * 1024;
	constexpr size_t sma_pool_amount = 10;
}  // namespace config

/**
Our allocator class that will be the general interface for any future
allocation. Inherits from the Curiously Recurring Singleton Pattern in order to
have just one global one.

@see CRSP
*/
class SmallMemoryAllocator : public CRSP<SmallMemoryAllocator> {
private:
    friend class CRSP<SmallMemoryAllocator>;

    SmallMemoryAllocator();
    ~SmallMemoryAllocator();

public:
    Pool * get_pool_for_size(size_t size);

    bool fits(size_t size) const;

    void* alloc(size_t size);

    template <typename T>
    T* alloc() {
        return reinterpret_cast<T*>(alloc(sizeof(T)));
    }

    void dealloc(void* elem);

    void print_status(std::ostream& stream);

    void print_status(void* elem, std::ostream& stream) const;

private:
	size_t get_pool_index(size_t size) const;

    std::array<std::unique_ptr<Pool>, config::sma_pool_amount> m_pool_array{};
	std::array<size_t, config::sma_pool_amount + 1> m_extra_requested_elements{};
    std::array<size_t, config::sma_pool_amount> m_max_allocated_elements{};

    byte* m_general_pool{ nullptr };
};

max_stars_repo_name: osor-io/NGene

max_issues_repo_name: osor-io/NGene

max_forks_repo_name: osor-io/NGene


text:
#pragma once

#if __has_include(<span>) && (__cplusplus > 201703L)
#include <span>
#else
#include <gsl/span>
namespace std {
template <class T, std::size_t Extent = gsl::dynamic_extent>
using span = gsl::span<T, Extent>;
}
#endif

max_stars_repo_name: VITObelgium/cpp-infra

max_issues_repo_name: VITObelgium/cpp-infra

max_forks_repo_name: VITObelgium/cpp-infra


text:
/*
   The following program solves the second-order nonlinear
   Van der Pol oscillator equation (see Landau/Paez 14.12, part 1),
   x"(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0
   This can be converted into a first order system suitable for
   use with the library by introducing a separate variable for
   the velocity, v = x'(t).  We assign x --> y[0] and v --> y[1].
   So the equations are:
   x' = v                  ==>  dy[0]/dt = f[0] = y[1]
   v' = -x + \mu v (1-x^2) ==>  dy[1]/dt = f[1] = -y[0] + mu*y[1]*(1-y[0]*y[0])
   x"(t) = - \mu x'(t) (x(t)^2 - 1) - x(t)

   x''(t) + x'(t) + x(t) = 0
   x''(t) = -x'(t) - x(t)

   x' = v ===> dy[0]/dt = f[0] = y[1]
   v' = -v - x
*/

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

/* function prototypes */
int rhs (double t, const double y[], double f[], void *params_ptr);
//int jacobian (double t, const double y[], double *dfdy, double dfdt[], void *params_ptr);

/*************************** main program ****************************/

int *jac;

int
main (void)
{
    int dimension = 2;		/* number of differential equations */

    double eps_abs = 1.e-8;	/* absolute error requested */
    double eps_rel = 1.e-10;	/* relative error requested */

    /* define the type of routine for making steps: */
    const gsl_odeiv_step_type *type_ptr = gsl_odeiv_step_rkf45;
    /* some other possibilities (see GSL manual):
       = gsl_odeiv_step_rk4;
       = gsl_odeiv_step_rkck;
       = gsl_odeiv_step_rk8pd;
       = gsl_odeiv_step_rk4imp;
       = gsl_odeiv_step_bsimp;
       = gsl_odeiv_step_gear1;
       = gsl_odeiv_step_gear2;
    */

    /*
       allocate/initialize the stepper, the control function, and the
       evolution function.
    */
	// 룬게쿠타 방식의 2계 미분 방정식으로 stepping을 진행합니다.
    gsl_odeiv_step *step_ptr = gsl_odeiv_step_alloc (type_ptr, dimension);
	// 허용 오차
    gsl_odeiv_control *control_ptr = gsl_odeiv_control_y_new (eps_abs, eps_rel);
	// 계산하는 차원
    gsl_odeiv_evolve *evolve_ptr = gsl_odeiv_evolve_alloc (dimension);

    gsl_odeiv_system my_system;	/* structure with the rhs function, etc. */

    double mu = 10;		/* parameter for the diffeq */
    double y[2];			/* current solution vector */

    double t, t_next;		/* current and next independent variable */
    double tmin, tmax, delta_t;	/* range of t and step size for output */

    double h = 1e-6;		/* starting step size for ode solver */

    /* load values into the my_system structure */
    my_system.function = rhs;	/* the right-hand-side functions dy[i]/dt */
    //my_system.jacobian = jacobian;	/* the Jacobian df[i]/dy[j] */
    my_system.jacobian = NULL;
    my_system.dimension = dimension;	/* number of diffeq's */
    //my_system.params = μ	/* parameters to pass to rhs and jacobian */
    my_system.params = NULL;

    tmin = 0.;			/* starting t value */
    tmax = 10.;			/* final t value */
    delta_t = 1.;

	// 주의점이라면 2계 미분 방정식은 최소 2개의 초기값이 필요하다는 부분을 주의!
    y[0] = 1.;			/* initial x value */
    y[1] = 0.;			/* initial v value */

    t = tmin;             /* initialize t */
    printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);	/* initial values */

    /* step to tmax from tmin */
    for (t_next = tmin + delta_t; t_next <= tmax; t_next += delta_t)
    {
        while (t < t_next)	/* evolve from t to t_next */
        {
            gsl_odeiv_evolve_apply (evolve_ptr, control_ptr, step_ptr,
                                    &my_system, &t, t_next, &h, y);
        }
		// y[0], y[1] 결과 확인이 필요함
        printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); /* print at t=t_next */
    }

    /* all done; free up the gsl_odeiv stuff */
    gsl_odeiv_evolve_free (evolve_ptr);
    gsl_odeiv_control_free (control_ptr);
    gsl_odeiv_step_free (step_ptr);

    return 0;
}

/*************************** rhs ****************************/
/*
   x''(t) + x'(t) + x(t) = 0
   x''(t) = -x'(t) - x(t)

   x' = v            ===> dy[0]/dt = f[0] = y[1]
   v' = -v - x       ===> dy[1]/dt = f[1] = -y[1] - y[0]
*/

// https://www.wolframalpha.com/input/?i=y%27%27+%2B+y%27+%2B+y+%3D+0%2C+y%280%29+%3D+1%2C+y%27%280%29+%3D+0, y(0) = 1
// https://www.wolframalpha.com/input/?i=%28sqrt%283%29+*+sin%28sqrt%283%29+%2F+2%29+%2B+3+*+cos%28sqrt%283%29+%2F+2%29%29+%2F+%283+*+sqrt%28e%29%29, y(1)
// (sqrt(3) * sin(sqrt(3) / 2) + 3 * cos(sqrt(3) / 2)) / (3 * sqrt(e)), y(1) = 0.6597001...
int
rhs (double t, const double v[], double f[], void *params_ptr)
//rhs (double t, const double y[], double f[], void *params_ptr)
{
	/*
    f[0] = y[1];
    f[1] = -y[0] + mu * y[1] * (1. - y[0] * y[0]);
	*/
	// * 코드 작성 팁: 미분항은 f에 배치합니다.
	// 1계와 마찬가지로 반드시 표준형이어야 합니다!
	// 첫번째 1차 미분이 무엇(치환)입니다 ~
	f[0] = v[1];
	// 두번째 2차 미분을 치환을 기반으로 어떻게 표현해야 하는지를 기록합니다.
	f[1] = -v[1] - v[0];
	/* Tip: 배열의 인덱스는 결국 미분의 횟수라고 생각해도 좋습니다 */
	// 현재 초기값이 있으니까 두 개의 근이 나오지 않고 한 개의 근이 나옴!

    return GSL_SUCCESS;		/* GSL_SUCCESS defined in gsl/errno.h as 0 */
}

/*************************** Jacobian ****************************/
/*
  Define the Jacobian matrix using GSL matrix routines.
  (see the GSL manual under "Ordinary Differential Equations")

  * params is a void pointer that is used in many GSL routines
  to pass parameters to a function
*/

#if 0
int
jacobian (double t, const double y[], double *dfdy,
          double dfdt[], void *params_ptr)
{
    /* get parameter(s) from params_ptr; here, just a double */
    double mu = *(double *) params_ptr;

    gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);

    gsl_matrix *m_ptr = &dfdy_mat.matrix;	/* m_ptr points to the matrix */

    /* fill the Jacobian matrix as shown */
    gsl_matrix_set (m_ptr, 0, 0, 0.0);	/* df[0]/dy[0] = 0 */
    gsl_matrix_set (m_ptr, 0, 1, 1.0);	/* df[0]/dy[1] = 1 */
    gsl_matrix_set (m_ptr, 1, 0, -2.0 * mu * y[0] * y[1] - 1.0); /* df[1]/dy[0] */
    gsl_matrix_set (m_ptr, 1, 1, -mu * (y[0] * y[0] - 1.0));     /* df[1]/dy[1] */

    /* set explicit t dependence of f[i] */
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;

    return GSL_SUCCESS;		/* GSL_SUCCESS defined in gsl/errno.h as 0 */
}
#endif

ext: c


text:
#ifdef __cplusplus
extern "C" {
#endif

#include "src/allvars.h"
#include "src/proto.h"
#include "src/tags.h"
#include <gsl/gsl_rng.h>

typedef struct {
    double mass;                                        /// mass
    double x, y, z;                                     /// position
    double vx, vy, vz;                                  /// velocity
} dynamics_state;

typedef struct {
    double mass;                                        /// mass
    double x, y, z;                                     /// position
    double vx, vy, vz;                                  /// velocity
    double u;                                           /// entropy
#ifdef MORRIS97VISC
    double alpha, dalphadt;				///viscosity
#endif
} sph_state;

void   begrun(void);
double second(void);

int found_particle(int index_of_the_particle, int *local_index);
void update_particle_map(void);

void hydro_state_at_point(FLOAT pos[3], FLOAT vel[3], FLOAT *h_out,
  FLOAT *ngb_out, FLOAT *dhsml_out, FLOAT *rho_out, FLOAT *rhov_out,
  FLOAT *rhov2_out, FLOAT *rhoe_out);

#ifdef __cplusplus
}
#endif

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
#include <petsc.h>
#include <petscmath.h>
#include "compressibleFlow.h"
#include "mesh.h"
#include "petscdmplex.h"
#include "petscts.h"

//MMS from Verification of a Compressible CFD Code using the Method of Manufactured Solutions, Christopher J. Roy,† Thomas M. Smith,‡ and Curtis C. Ober§

// Define
#define Pi PETSC_PI
#define Sin PetscSinReal
#define Cos PetscCosReal
#define Power PetscPowReal

typedef struct {
    PetscReal phiO;
    PetscReal phiX;
    PetscReal phiY;
    PetscReal phiZ;
    PetscReal aPhiX;
    PetscReal aPhiY;
    PetscReal aPhiZ;
} PhiConstants;

typedef struct {
    PetscInt dim;
    PhiConstants rho;
    PhiConstants u;
    PhiConstants v;
    PhiConstants w;
    PhiConstants p;
    PetscReal L;
    PetscReal gamma;
    PetscReal R;
    PetscReal mu;
} Constants;

typedef struct {
    Constants constants;
    FlowData flowData;
} ProblemSetup;

static PetscErrorCode EulerExact(PetscInt dim, PetscReal time, const PetscReal xyz[], PetscInt Nf, PetscScalar *u, void *ctx) {
    PetscFunctionBeginUser;

    Constants *constants = (Constants *)ctx;
    PetscReal L = constants->L;
    PetscReal gamma = constants->gamma;

    PetscReal rhoO = constants->rho.phiO;
    PetscReal rhoX = constants->rho.phiX;
    PetscReal rhoY = constants->rho.phiY;
    PetscReal rhoZ = constants->rho.phiZ;
    PetscReal aRhoX = constants->rho.aPhiX;
    PetscReal aRhoY = constants->rho.aPhiY;
    PetscReal aRhoZ = constants->rho.aPhiZ;

    PetscReal uO = constants->u.phiO;
    PetscReal uX = constants->u.phiX;
    PetscReal uY = constants->u.phiY;
    PetscReal uZ = constants->u.phiZ;
    PetscReal aUX = constants->u.aPhiX;
    PetscReal aUY = constants->u.aPhiY;
    PetscReal aUZ = constants->u.aPhiZ;

    PetscReal vO = constants->v.phiO;
    PetscReal vX = constants->v.phiX;
    PetscReal vY = constants->v.phiY;
    PetscReal vZ = constants->v.phiZ;
    PetscReal aVX = constants->v.aPhiX;
    PetscReal aVY = constants->v.aPhiY;
    PetscReal aVZ = constants->v.aPhiZ;

    PetscReal wO = constants->w.phiO;
    PetscReal wX = constants->w.phiX;
    PetscReal wY = constants->w.phiY;
    PetscReal wZ = constants->w.phiZ;
    PetscReal aWX = constants->w.aPhiX;
    PetscReal aWY = constants->w.aPhiY;
    PetscReal aWZ = constants->w.aPhiZ;

    PetscReal pO = constants->p.phiO;
    PetscReal pX = constants->p.phiX;
    PetscReal pY = constants->p.phiY;
    PetscReal pZ = constants->p.phiZ;
    PetscReal aPX = constants->p.aPhiX;
    PetscReal aPY = constants->p.aPhiY;
    PetscReal aPZ = constants->p.aPhiZ;

    PetscReal x = xyz[0];
    PetscReal y = xyz[1];
    PetscReal z = dim > 2? xyz[2] : 0.0;

    u[RHO] = rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L);
    u[RHOE] = (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*((pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))/((-1. + gamma)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                (Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2) + Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2) + Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L),2))/2.);

    u[RHOU + 0] = (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*
           (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L));
    u[RHOU + 1] = (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
           (vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L));

    if(dim > 2){
        u[RHOU + 2] = (wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*
               (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L));
    }

    PetscFunctionReturn(0);
}

static PetscErrorCode EulerExactTimeDerivative(PetscInt dim, PetscReal time, const PetscReal xyz[], PetscInt Nf, PetscScalar *u, void *ctx) {
    PetscFunctionBeginUser;
    u[RHO] = 0.0;
    u[RHOE] = 0.0;
    u[RHOU + 0] = 0.0;
    u[RHOU + 1] = 0.0;
    if(dim > 2) {
        u[RHOU + 2] = 0.0;
    }

    PetscFunctionReturn(0);
}

static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx) {
    PetscFunctionBeginUser;
    PetscErrorCode     ierr;

    // Get the DM
    DM dm;
    ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);

    ierr = VecViewFromOptions(u, NULL, "-sol_view");CHKERRQ(ierr);
    ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "TS at %f\n", time);CHKERRQ(ierr);

    // Compute the error
    void            *exactCtxs[1];
    PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
    PetscDS          ds;
    ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);

    // Get the exact solution
    ierr = PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &exactCtxs[0]);CHKERRQ(ierr);

    // Create an vector to hold the exact solution
    Vec exactVec;
    ierr = VecDuplicate(u, &exactVec);CHKERRQ(ierr);
    ierr = DMProjectFunction(dm,time,exactFuncs,exactCtxs,INSERT_ALL_VALUES,exactVec);CHKERRQ(ierr);

    // For each component, compute the l2 norms
    ierr = VecAXPY(exactVec, -1.0, u);CHKERRQ(ierr);

    PetscReal ferrors[4];
    ierr = VecSetBlockSize(exactVec, 4);CHKERRQ(ierr);

    // Compute the L2 Difference
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t\n", (int) step, (double) time);CHKERRQ(ierr);
    ierr =  VecStrideNormAll(exactVec, NORM_2,ferrors);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_2 Error: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double) ferrors[0], (double) ferrors[1], (double) ferrors[2], (double) ferrors[3]);CHKERRQ(ierr);
    ierr =  VecStrideNormAll(exactVec, NORM_INFINITY,ferrors);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_Inf Error: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double) ferrors[0], (double) ferrors[1], (double) ferrors[2], (double) ferrors[3]);CHKERRQ(ierr);

    ierr = VecDestroy(&exactVec);CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

static PetscErrorCode SourceRho(PetscInt dim, PetscReal time, const PetscReal xyz[], PetscInt Nf, PetscScalar *u, void *ctx) {
    PetscFunctionBeginUser;

    Constants *constants = (Constants *)ctx;
    PetscReal L = constants->L;
    PetscReal gamma = constants->gamma;

    PetscReal rhoO = constants->rho.phiO;
    PetscReal rhoX = constants->rho.phiX;
    PetscReal rhoY = constants->rho.phiY;
    PetscReal rhoZ = constants->rho.phiZ;
    PetscReal aRhoX = constants->rho.aPhiX;
    PetscReal aRhoY = constants->rho.aPhiY;
    PetscReal aRhoZ = constants->rho.aPhiZ;

    PetscReal uO = constants->u.phiO;
    PetscReal uX = constants->u.phiX;
    PetscReal uY = constants->u.phiY;
    PetscReal uZ = constants->u.phiZ;
    PetscReal aUX = constants->u.aPhiX;
    PetscReal aUY = constants->u.aPhiY;
    PetscReal aUZ = constants->u.aPhiZ;

    PetscReal vO = constants->v.phiO;
    PetscReal vX = constants->v.phiX;
    PetscReal vY = constants->v.phiY;
    PetscReal vZ = constants->v.phiZ;
    PetscReal aVX = constants->v.aPhiX;
    PetscReal aVY = constants->v.aPhiY;
    PetscReal aVZ = constants->v.aPhiZ;

    PetscReal wO = constants->w.phiO;
    PetscReal wX = constants->w.phiX;
    PetscReal wY = constants->w.phiY;
    PetscReal wZ = constants->w.phiZ;
    PetscReal aWX = constants->w.aPhiX;
    PetscReal aWY = constants->w.aPhiY;
    PetscReal aWZ = constants->w.aPhiZ;

    PetscReal pO = constants->p.phiO;
    PetscReal pX = constants->p.phiX;
    PetscReal pY = constants->p.phiY;
    PetscReal pZ = constants->p.phiZ;
    PetscReal aPX = constants->p.aPhiX;
    PetscReal aPY = constants->p.aPhiY;
    PetscReal aPZ = constants->p.aPhiZ;

    PetscReal x = xyz[0];
    PetscReal y = xyz[1];
    PetscReal z = dim > 2? xyz[2] : 0.0;

    u[0] = (aRhoX*Pi*rhoX*Cos((aRhoX*Pi*x)/L)*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) +
                                               uX*Sin((aUX*Pi*x)/L)))/L + (aRhoZ*Pi*rhoZ*Cos((aRhoZ*Pi*z)/L)*
                                                                           (wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L)))/L +
           (aUX*Pi*uX*Cos((aUX*Pi*x)/L)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) +
                                         rhoZ*Sin((aRhoZ*Pi*z)/L)))/L + (aVY*Pi*vY*Cos((aVY*Pi*y)/L)*
                                                                         (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L)))/L -
           (aRhoY*Pi*rhoY*Sin((aRhoY*Pi*y)/L)*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) +
                                               vZ*Sin((aVZ*Pi*z)/L)))/L - (aWZ*Pi*wZ*
                                                                           (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
                                                                           Sin((aWZ*Pi*z)/L))/L;
    PetscFunctionReturn(0);
}

static PetscErrorCode SourceRhoU(PetscInt dim, PetscReal time, const PetscReal xyz[], PetscInt Nf, PetscScalar *u, void *ctx) {
    PetscFunctionBeginUser;

    Constants *constants = (Constants *)ctx;
    PetscReal L = constants->L;
    PetscReal gamma = constants->gamma;

    PetscReal rhoO = constants->rho.phiO;
    PetscReal rhoX = constants->rho.phiX;
    PetscReal rhoY = constants->rho.phiY;
    PetscReal rhoZ = constants->rho.phiZ;
    PetscReal aRhoX = constants->rho.aPhiX;
    PetscReal aRhoY = constants->rho.aPhiY;
    PetscReal aRhoZ = constants->rho.aPhiZ;

    PetscReal uO = constants->u.phiO;
    PetscReal uX = constants->u.phiX;
    PetscReal uY = constants->u.phiY;
    PetscReal uZ = constants->u.phiZ;
    PetscReal aUX = constants->u.aPhiX;
    PetscReal aUY = constants->u.aPhiY;
    PetscReal aUZ = constants->u.aPhiZ;

    PetscReal vO = constants->v.phiO;
    PetscReal vX = constants->v.phiX;
    PetscReal vY = constants->v.phiY;
    PetscReal vZ = constants->v.phiZ;
    PetscReal aVX = constants->v.aPhiX;
    PetscReal aVY = constants->v.aPhiY;
    PetscReal aVZ = constants->v.aPhiZ;

    PetscReal wO = constants->w.phiO;
    PetscReal wX = constants->w.phiX;
    PetscReal wY = constants->w.phiY;
    PetscReal wZ = constants->w.phiZ;
    PetscReal aWX = constants->w.aPhiX;
    PetscReal aWY = constants->w.aPhiY;
    PetscReal aWZ = constants->w.aPhiZ;

    PetscReal pO = constants->p.phiO;
    PetscReal pX = constants->p.phiX;
    PetscReal pY = constants->p.phiY;
    PetscReal pZ = constants->p.phiZ;
    PetscReal aPX = constants->p.aPhiX;
    PetscReal aPY = constants->p.aPhiY;
    PetscReal aPZ = constants->p.aPhiZ;

    PetscReal x = xyz[0];
    PetscReal y = xyz[1];
    PetscReal z = dim > 2? xyz[2] : 0.0;

    u[0] = -((aPX*Pi*pX*Sin((aPX*Pi*x)/L))/L) + (aRhoX*Pi*rhoX*Cos((aRhoX*Pi*x)/L)*
                                                 Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2))/L +
           (aRhoZ*Pi*rhoZ*Cos((aRhoZ*Pi*z)/L)*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) +
                                               uX*Sin((aUX*Pi*x)/L))*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L)))/
           L + (2*aUX*Pi*uX*Cos((aUX*Pi*x)/L)*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) +
                                               uX*Sin((aUX*Pi*x)/L))*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) +
                                                                      rhoZ*Sin((aRhoZ*Pi*z)/L)))/L + (aVY*Pi*vY*Cos((aVY*Pi*y)/L)*
                                                                                (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*
                                                                                (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L)))/L -
           (aUZ*Pi*uZ*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*
            (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
            Sin((aUZ*Pi*z)/L))/L - (aRhoY*Pi*rhoY*
                                    (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*Sin((aRhoY*Pi*y)/L)*
                                    (vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/L -
           (aUY*Pi*uY*Sin((aUY*Pi*y)/L)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) +
                                         rhoZ*Sin((aRhoZ*Pi*z)/L))*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) +
                                                                    vZ*Sin((aVZ*Pi*z)/L)))/L - (aWZ*Pi*wZ*
                                                                                (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*
                                                                                (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
                                                                                Sin((aWZ*Pi*z)/L))/L;

    u[1] = (aPY*Pi*pY*Cos((aPY*Pi*y)/L))/L - (aVX*Pi*vX*
                                              (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*Sin((aVX*Pi*x)/L)*
                                              (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L)))/L +
           (aVZ*Pi*vZ*Cos((aVZ*Pi*z)/L)*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) +
                                         wY*Sin((aWY*Pi*y)/L))*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) +
                                                                rhoZ*Sin((aRhoZ*Pi*z)/L)))/L + (aRhoX*Pi*rhoX*Cos((aRhoX*Pi*x)/L)*
                                                                                (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*
                                                                                (vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/L +
           (aRhoZ*Pi*rhoZ*Cos((aRhoZ*Pi*z)/L)*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) +
                                               wY*Sin((aWY*Pi*y)/L))*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/
           L + (aUX*Pi*uX*Cos((aUX*Pi*x)/L)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) +
                                             rhoZ*Sin((aRhoZ*Pi*z)/L))*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) +
                                                                        vZ*Sin((aVZ*Pi*z)/L)))/L + (2*aVY*Pi*vY*Cos((aVY*Pi*y)/L)*
                                                                                (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
                                                                                (vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/L -
           (aRhoY*Pi*rhoY*Sin((aRhoY*Pi*y)/L)*Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) +
                                                    vZ*Sin((aVZ*Pi*z)/L),2))/L - (aWZ*Pi*wZ*
                                                                                (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
                                                                                (vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L))*Sin((aWZ*Pi*z)/L))/L;

    if(dim >2){
        u[2] = (aRhoX*Pi*rhoX*Cos((aRhoX*Pi*x)/L)*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) +
                                                   uX*Sin((aUX*Pi*x)/L))*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L)))/
               L + (aRhoZ*Pi*rhoZ*Cos((aRhoZ*Pi*z)/L)*
                    Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2))/L -
               (aPZ*Pi*pZ*Sin((aPZ*Pi*z)/L))/L + (aWX*Pi*wX*Cos((aWX*Pi*x)/L)*
                                                  (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*
                                                  (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L)))/L +
               (aUX*Pi*uX*Cos((aUX*Pi*x)/L)*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) +
                                             wY*Sin((aWY*Pi*y)/L))*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) +
                                                                    rhoZ*Sin((aRhoZ*Pi*z)/L)))/L + (aVY*Pi*vY*Cos((aVY*Pi*y)/L)*
                                                                                (wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*
                                                                                (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L)))/L -
               (aRhoY*Pi*rhoY*Sin((aRhoY*Pi*y)/L)*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) +
                                                   wY*Sin((aWY*Pi*y)/L))*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/
               L + (aWY*Pi*wY*Cos((aWY*Pi*y)/L)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) +
                                                 rhoZ*Sin((aRhoZ*Pi*z)/L))*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) +
                                                                            vZ*Sin((aVZ*Pi*z)/L)))/L - (2*aWZ*Pi*wZ*
                                                                                (wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*
                                                                                (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
                                                                                Sin((aWZ*Pi*z)/L))/L;
    }
    PetscFunctionReturn(0);
}

static PetscErrorCode SourceRhoE(PetscInt dim, PetscReal time, const PetscReal xyz[], PetscInt Nf, PetscScalar *u, void *ctx) {
    PetscFunctionBeginUser;

    Constants *constants = (Constants *)ctx;
    PetscReal L = constants->L;
    PetscReal gamma = constants->gamma;

    PetscReal rhoO = constants->rho.phiO;
    PetscReal rhoX = constants->rho.phiX;
    PetscReal rhoY = constants->rho.phiY;
    PetscReal rhoZ = constants->rho.phiZ;
    PetscReal aRhoX = constants->rho.aPhiX;
    PetscReal aRhoY = constants->rho.aPhiY;
    PetscReal aRhoZ = constants->rho.aPhiZ;

    PetscReal uO = constants->u.phiO;
    PetscReal uX = constants->u.phiX;
    PetscReal uY = constants->u.phiY;
    PetscReal uZ = constants->u.phiZ;
    PetscReal aUX = constants->u.aPhiX;
    PetscReal aUY = constants->u.aPhiY;
    PetscReal aUZ = constants->u.aPhiZ;

    PetscReal vO = constants->v.phiO;
    PetscReal vX = constants->v.phiX;
    PetscReal vY = constants->v.phiY;
    PetscReal vZ = constants->v.phiZ;
    PetscReal aVX = constants->v.aPhiX;
    PetscReal aVY = constants->v.aPhiY;
    PetscReal aVZ = constants->v.aPhiZ;

    PetscReal wO = constants->w.phiO;
    PetscReal wX = constants->w.phiX;
    PetscReal wY = constants->w.phiY;
    PetscReal wZ = constants->w.phiZ;
    PetscReal aWX = constants->w.aPhiX;
    PetscReal aWY = constants->w.aPhiY;
    PetscReal aWZ = constants->w.aPhiZ;

    PetscReal pO = constants->p.phiO;
    PetscReal pX = constants->p.phiX;
    PetscReal pY = constants->p.phiY;
    PetscReal pZ = constants->p.phiZ;
    PetscReal aPX = constants->p.aPhiX;
    PetscReal aPY = constants->p.aPhiY;
    PetscReal aPZ = constants->p.aPhiZ;

    PetscReal x = xyz[0];
    PetscReal y = xyz[1];
    PetscReal z = dim > 2? xyz[2] : 0.0;

    u[0] = -((aPX*Pi*pX*Sin((aPX*Pi*x)/L)*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L)))/L) + (aUX*Pi*uX*Cos((aUX*Pi*x)/L)*(pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L)))/L +
           (aVY*Pi*vY*Cos((aVY*Pi*y)/L)*(pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L)))/L - (aPZ*Pi*pZ*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*Sin((aPZ*Pi*z)/L))/L +
           (aPY*Pi*pY*Cos((aPY*Pi*y)/L)*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/L + (rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L))*
                                                                                ((aRhoY*Pi*rhoY*(pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))*Sin((aRhoY*Pi*y)/L))/((-1. + gamma)*L*Power(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L),2)) +
                                                                                (aPY*Pi*pY*Cos((aPY*Pi*y)/L))/((-1. + gamma)*L*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                ((-2*aUY*Pi*uY*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*Sin((aUY*Pi*y)/L))/L + (2*aWY*Pi*wY*Cos((aWY*Pi*y)/L)*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L)))/L +
                                                                                (2*aVY*Pi*vY*Cos((aVY*Pi*y)/L)*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/L)/2.) + (uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
                                                                                (-((aRhoX*Pi*rhoX*Cos((aRhoX*Pi*x)/L)*(pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L)))/((-1. + gamma)*L*Power(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L),2))) -
                                                                                (aPX*Pi*pX*Sin((aPX*Pi*x)/L))/((-1. + gamma)*L*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                ((2*aUX*Pi*uX*Cos((aUX*Pi*x)/L)*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L)))/L + (2*aWX*Pi*wX*Cos((aWX*Pi*x)/L)*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L)))/L -
                                                                                (2*aVX*Pi*vX*Sin((aVX*Pi*x)/L)*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/L)/2.) + (aRhoX*Pi*rhoX*Cos((aRhoX*Pi*x)/L)*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*
                                                                                ((pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))/((-1. + gamma)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                (Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2) + Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2) + Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L),2))/2.))/L +
           (aRhoZ*Pi*rhoZ*Cos((aRhoZ*Pi*z)/L)*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*((pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))/((-1. + gamma)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                (Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2) + Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2) + Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L),2))/2.))/L +
           (aUX*Pi*uX*Cos((aUX*Pi*x)/L)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*((pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))/((-1. + gamma)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                (Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2) + Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2) + Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L),2))/2.))/L +
           (aVY*Pi*vY*Cos((aVY*Pi*y)/L)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*((pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))/((-1. + gamma)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                (Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2) + Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2) + Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L),2))/2.))/L -
           (aRhoY*Pi*rhoY*Sin((aRhoY*Pi*y)/L)*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L))*((pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))/((-1. + gamma)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                (Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2) + Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2) + Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L),2))/2.))/L -
           (aWZ*Pi*wZ*(pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))*Sin((aWZ*Pi*z)/L))/L - (aWZ*Pi*wZ*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
                                                                                ((pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L))/((-1. + gamma)*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
                                                                                (Power(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L),2) + Power(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L),2) + Power(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L),2))/2.)*Sin((aWZ*Pi*z)/L))/L +
           (wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))*
           (-((aRhoZ*Pi*rhoZ*Cos((aRhoZ*Pi*z)/L)*(pO + pX*Cos((aPX*Pi*x)/L) + pZ*Cos((aPZ*Pi*z)/L) + pY*Sin((aPY*Pi*y)/L)))/((-1. + gamma)*L*Power(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L),2))) -
            (aPZ*Pi*pZ*Sin((aPZ*Pi*z)/L))/((-1. + gamma)*L*(rhoO + rhoY*Cos((aRhoY*Pi*y)/L) + rhoX*Sin((aRhoX*Pi*x)/L) + rhoZ*Sin((aRhoZ*Pi*z)/L))) +
            ((-2*aUZ*Pi*uZ*(uO + uY*Cos((aUY*Pi*y)/L) + uZ*Cos((aUZ*Pi*z)/L) + uX*Sin((aUX*Pi*x)/L))*Sin((aUZ*Pi*z)/L))/L + (2*aVZ*Pi*vZ*Cos((aVZ*Pi*z)/L)*(vO + vX*Cos((aVX*Pi*x)/L) + vY*Sin((aVY*Pi*y)/L) + vZ*Sin((aVZ*Pi*z)/L)))/L -
             (2*aWZ*Pi*wZ*(wO + wZ*Cos((aWZ*Pi*z)/L) + wX*Sin((aWX*Pi*x)/L) + wY*Sin((aWY*Pi*y)/L))*Sin((aWZ*Pi*z)/L))/L)/2.);
    PetscFunctionReturn(0);
}

static PetscErrorCode PhysicsBoundary_Euler(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx) {
    PetscFunctionBeginUser;
    Constants *constants = (Constants *)ctx;

    // Offset the calc assuming the cells are square
    PetscReal x[3];

    for(PetscInt i =0; i < constants->dim; i++){
        x[i] = c[i] + n[i]*0.5;
    }

    EulerExact(constants->dim, time, x, 0, a_xG, ctx);
    PetscFunctionReturn(0);
}

static PetscErrorCode ComputeRHSWithSourceTerms(DM dm, PetscReal time, Vec locXVec, Vec globFVec, void *ctx){
    PetscFunctionBeginUser;
    PetscErrorCode ierr;
    ProblemSetup *setup = (ProblemSetup *)ctx;

    // Call the flux calculation
    ierr = DMPlexTSComputeRHSFunctionFVM(dm, time, locXVec, globFVec, setup->flowData);CHKERRQ(ierr);

    // Convert the dm to a plex
    DM plex;
    DMConvert(dm, DMPLEX, &plex);

    // Extract the cell geometry, and the dm that holds the information
    Vec cellgeom;
    DM dmCell;
    const PetscScalar *cgeom;
    ierr = DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL);CHKERRQ(ierr);
    ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
    ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);

    // Get the cell start and end for the fv cells
    PetscInt cStart, cEnd;
    ierr = DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd);CHKERRQ(ierr);

    // create a local f vector
    Vec locFVec;
    PetscScalar  *locFArray;
    ierr = DMGetLocalVector(dm, &locFVec);CHKERRQ(ierr);
    ierr = VecZeroEntries(locFVec);CHKERRQ(ierr);
    ierr = VecGetArray(locFVec, &locFArray);CHKERRQ(ierr);

    // get the current values
    const PetscScalar      *locXArray;
    ierr = VecGetArrayRead(locXVec, &locXArray);CHKERRQ(ierr);

    PetscInt rank;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

    // March over each cell volume
    for (PetscInt c = cStart; c < cEnd; ++c) {
        PetscFVCellGeom       *cg;
        const PetscReal           *xc;
        PetscReal           *fc;

        ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
        ierr = DMPlexPointLocalFieldRead(plex, c, 0, locXArray, &xc);CHKERRQ(ierr);
        ierr = DMPlexPointGlobalFieldRef(plex, c, 0, locFArray, &fc);CHKERRQ(ierr);

        if(fc) {  // must be real cell and not ghost
            SourceRho(setup->constants.dim, time, cg->centroid, 0, fc + RHO, &setup->constants);
            SourceRhoE(setup->constants.dim, time, cg->centroid, 0, fc + RHOE, &setup->constants);
            SourceRhoU(setup->constants.dim, time, cg->centroid, 0, fc + RHOU, &setup->constants);
        }
    }

    // restore the cell geometry
    ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(locXVec, &locXArray);CHKERRQ(ierr);
    ierr = VecRestoreArray(locFVec, &locFArray);CHKERRQ(ierr);

    ierr = DMLocalToGlobalBegin(dm, locFVec, ADD_VALUES, globFVec);CHKERRQ(ierr);
    ierr = DMLocalToGlobalEnd(dm, locFVec, ADD_VALUES, globFVec);CHKERRQ(ierr);
    ierr = DMRestoreLocalVector(dm, &locFVec);CHKERRQ(ierr);

//    {// check rhs,
//        // temp read current residual
//        const PetscScalar *currentFArray;
//        ierr = VecGetArrayRead(globFVec, &currentFArray);CHKERRQ(ierr);
//        ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
//
//        // March over each cell volume
//        for (PetscInt c = cStart; c < cEnd; ++c) {
//            PetscFVCellGeom       *cg;
//            const PetscReal           *fcCurrent;
//
//            ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
//            ierr = DMPlexPointGlobalFieldRead(plex, c, 0, currentFArray, &fcCurrent);CHKERRQ(ierr);
//
//            if(fcCurrent) {  // must be real cell and not ghost
//                if(PetscAbsReal(cg->centroid[0] - .5 ) < 1E-8 && PetscAbsReal(cg->centroid[1] - .5 )  < 1E-8){
//                    printf("Residual(%f, %f): %f %f %f %f\n", cg->centroid[0], cg->centroid[1], fcCurrent[0], fcCurrent[1], fcCurrent[2], fcCurrent[3]);
//                }
//            }
//        }
//
//        // temp return current residual
//        ierr = VecRestoreArrayRead(globFVec, &currentFArray);CHKERRQ(ierr);
//        ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
//    }

    PetscFunctionReturn(0);
}

PetscErrorCode ComputeRHS(TS ts, DM dm, PetscReal t, Vec u, PetscInt blockSize, PetscReal residualNorm2[], PetscReal residualNormInf[], PetscReal start[], PetscReal end[])
{
    MPI_Comm       comm;
    Vec            r;
    PetscErrorCode ierr;

    PetscFunctionBeginUser;
    Vec sol;
    ierr = VecDuplicate(u, &sol);CHKERRQ(ierr);
    ierr = VecCopy(u, sol);CHKERRQ(ierr);

    ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr);
    ierr = DMComputeExactSolution(dm, t, sol, NULL);CHKERRQ(ierr);
    ierr = VecDuplicate(u, &r);CHKERRQ(ierr);
    ierr = TSComputeRHSFunction(ts, t, sol, r);CHKERRQ(ierr);

    // zero out the norms
    for(PetscInt b =0; b < blockSize; b++){
        residualNorm2[b] = 0.0;
        residualNormInf[b] = 0.0;
    }

    {//March over each cell
        // Extract the cell geometry, and the dm that holds the information
        Vec cellgeom;
        DM dmCell;
        PetscInt dim;
        const PetscScalar *cgeom;
        ierr = DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL);CHKERRQ(ierr);
        ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
        ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);

        // Get the cell start and end for the fv cells
        PetscInt cStart, cEnd;
        ierr = DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd);CHKERRQ(ierr);

        // temp read current residual
        const PetscScalar *currentRHS;
        ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
        ierr = VecGetArrayRead(r, &currentRHS);CHKERRQ(ierr);

        // Count up the cells
        PetscInt count = 0;

        // March over each cell volume
        for (PetscInt c = cStart; c < cEnd; ++c) {
            PetscFVCellGeom       *cg;
            const PetscReal           *rhsCurrent;

            ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);CHKERRQ(ierr);
            ierr = DMPlexPointGlobalFieldRead(dm, c, 0, currentRHS, &rhsCurrent);CHKERRQ(ierr);

            PetscBool countCell = PETSC_TRUE;
            for(PetscInt d =0; d < dim; d++){
                countCell = cg->centroid[d] < start[d] || cg->centroid[d] > end[d] ? PETSC_FALSE : countCell;
            }

            if(rhsCurrent && countCell ) {  // must be real cell and not ghost
                for(PetscInt b =0; b < blockSize; b++){
                    residualNorm2[b] += PetscSqr(rhsCurrent[b]);
                    residualNormInf[b] = PetscMax(residualNormInf[b], PetscAbs(rhsCurrent[b]));
                }
                count++;
            }
        }

        // normalize the norm2
        for(PetscInt b =0; b < blockSize; b++){
            residualNorm2[b] = PetscSqrtReal(residualNorm2[b]/count);
        }

        // temp return current residual
        ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
        ierr = VecRestoreArrayRead(r, &currentRHS);CHKERRQ(ierr);
    }

    ierr = VecDestroy(&sol);CHKERRQ(ierr);
    ierr = VecDestroy(&r);CHKERRQ(ierr);
    PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{

    // Setup the problem
    Constants constants;

    // sub sonic
    constants.dim = 2;
    constants.L = 1.0;
    constants.gamma = 1.4;
    constants.R = 287.0;
    constants.mu = 10;

    constants.rho.phiO = 1.0;
    constants.rho.phiX = 0.15;
    constants.rho.phiY = -0.1;
    constants.rho.phiZ = 0.0;
    constants.rho.aPhiX = 1.0;
    constants.rho.aPhiY = 0.5;
    constants.rho.aPhiZ = 0.0;

    constants.u.phiO = 70;
    constants.u.phiX = 5;
    constants.u.phiY = -7;
    constants.u.phiZ = 0.0;
    constants.u.aPhiX = 1.5;
    constants.u.aPhiY = 0.6;
    constants.u.aPhiZ = 0.0;

    constants.v.phiO = 90;
    constants.v.phiX = -15;
    constants.v.phiY = -8.5;
    constants.v.phiZ = 0.0;
    constants.v.aPhiX = 0.5;
    constants.v.aPhiY = 2.0/3.0;
    constants.v.aPhiZ = 0.0;

    constants.w.phiO = 0.0;
    constants.w.phiX = 0.0;
    constants.w.phiY = 0.0;
    constants.w.phiZ = 0.0;
    constants.w.aPhiX = 0.0;
    constants.w.aPhiY = 0.0;
    constants.w.aPhiZ = 0.0;

    constants.p.phiO = 1E5;
    constants.p.phiX = 0.2E5;
    constants.p.phiY = 0.5E5;
    constants.p.phiZ = 0.0;
    constants.p.aPhiX = 2.0;
    constants.p.aPhiY = 1.0;
    constants.p.aPhiZ = 0.0;

    PetscErrorCode ierr;
    // create the mesh
    // setup the ts
    DM dm;                 /* problem definition */
    TS ts;                 /* timestepper */

    // initialize petsc and mpi
    PetscInitialize(&argc, &argv, NULL, "HELP");

    // Create a ts
    ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
    ierr = TSSetProblemType(ts, TS_NONLINEAR);CHKERRQ(ierr);
    ierr = TSSetType(ts, TSEULER);CHKERRQ(ierr);
    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);

    // Create a mesh
    // hard code the problem setup
    PetscReal start[] = {0.0, 0.0};
    PetscReal end[] = {constants.L, constants.L};
    PetscInt nx[] = {4, 4};
    DMBoundaryType bcType[] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, constants.dim, PETSC_FALSE, nx, start, end, bcType, PETSC_TRUE, &dm);CHKERRQ(ierr);

    // Setup the flow data
    FlowData flowData;     /* store some of the flow data*/
    ierr = FlowCreate(&flowData);CHKERRQ(ierr);

    // Combine the flow data
    ProblemSetup problemSetup;
    problemSetup.flowData = flowData;
    problemSetup.constants = constants;

    //Setup
    CompressibleFlow_SetupDiscretization(flowData, &dm);

    // Add in the flow parameters
    PetscScalar params[TOTAL_COMPRESSIBLE_FLOW_PARAMETERS];
    params[CFL] = 0.5;
    params[GAMMA] = constants.gamma;

    // set up the finite volume fluxes
    CompressibleFlow_StartProblemSetup(flowData, TOTAL_COMPRESSIBLE_FLOW_PARAMETERS, params);
    DMView(flowData->dm, PETSC_VIEWER_STDERR_SELF);
    // Add in any boundary conditions
    PetscDS prob;
    ierr = DMGetDS(flowData->dm, &prob);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    const PetscInt idsLeft[]= {1, 2, 3, 4};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall left", "Face Sets", 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler, NULL, 4, idsLeft, &constants);CHKERRQ(ierr);

    // Complete the problem setup
    ierr = CompressibleFlow_CompleteProblemSetup(flowData, ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Override the flow calc for now
    ierr = DMTSSetRHSFunctionLocal(flowData->dm, ComputeRHSWithSourceTerms, &problemSetup);CHKERRQ(ierr);

    // Name the flow field
    ierr = PetscObjectSetName(((PetscObject)flowData->flowField), "Numerical Solution");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup the TS
    ierr = TSSetFromOptions(ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSMonitorSet(ts, MonitorError, &constants, NULL);CHKERRQ(ierr);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetMaxTime(ts, 0.01);CHKERRQ(ierr);

    // set the initial conditions
    PetscErrorCode     (*func[2]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {EulerExact};
    void* ctxs[1] ={&constants};
    ierr    = DMProjectFunction(flowData->dm,0.0,func,ctxs,INSERT_ALL_VALUES,flowData->flowField);CHKERRQ(ierr);

    // for the mms, add the exact solution
    ierr = PetscDSSetExactSolution(prob, 0, EulerExact, &constants);CHKERRQ(ierr);
    ierr = PetscDSSetExactSolutionTimeDerivative(prob, 0, EulerExactTimeDerivative, &constants);CHKERRQ(ierr);

    // Output the mesh
    ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);

    TSSetMaxSteps(ts, 1);
    ierr = TSSolve(ts,flowData->flowField);CHKERRQ(ierr);

    // Check the current residual
    PetscReal l2Residual[4];
    PetscReal infResidual[4];

    // Only take the residual over the central 1/3
    PetscReal resStart[2] = {1.0/3.0, 1.0/3.0};
    PetscReal resEnd[2] = {2.0/3.0, 2.0/3.0};

    ierr = ComputeRHS(ts, flowData->dm, 0.0, flowData->flowField, 4, l2Residual, infResidual, resStart, resEnd);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_2 Residual: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double) l2Residual[0], (double) l2Residual[1], (double) l2Residual[2], (double) l2Residual[3]);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "\tL_Inf Residual: [%2.3g, %2.3g, %2.3g, %2.3g]\n", (double) infResidual[0], (double) infResidual[1], (double) infResidual[2], (double) infResidual[3]);CHKERRQ(ierr);

    PetscReal time = 0.0;

//    {
//        Vec sol;
//        VecDuplicate(flowData->flowField, &sol);
//        VecCopy(flowData->flowField, sol);
//
//        SNES snes;
//        ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr);
//        ierr = DMSNESCheckDiscretization(snes, flowData->dm, time, sol, -1.0, NULL);CHKERRQ(ierr);
//
//        Vec u_t;
//        ierr = DMGetGlobalVector(flowData->dm, &u_t);CHKERRQ(ierr);
//        ierr = DMTSCheckResidual(ts, flowData->dm, time, sol, u_t, -1.0, NULL);CHKERRQ(ierr);
//        ierr = DMRestoreGlobalVector(flowData->dm, &u_t);CHKERRQ(ierr);
//
//        VecDestroy(&sol);
//    }

//    {
//        Vec sol;
//        VecDuplicate(flowData->flowField, &sol);
//        VecCopy(flowData->flowField, sol);
//
//        SNES snes;
//        ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr);
//        ierr = DMSNESCheckDiscretization(snes, flowData->dm, time, sol, -1.0, NULL);CHKERRQ(ierr);
//
//        Vec u_t;
//        ierr = DMGetGlobalVector(flowData->dm, &u_t);CHKERRQ(ierr);
//        ierr = DMTSCheckResidual(ts, flowData->dm, time, sol, u_t, -1.0, NULL);CHKERRQ(ierr);
//        ierr = DMRestoreGlobalVector(flowData->dm, &u_t);CHKERRQ(ierr);
//
//        VecDestroy(&sol);
//    }

    return PetscFinalize();

}

max_stars_repo_name: mmcgurn/MattFlowCases

max_issues_repo_name: mmcgurn/MattFlowCases

max_forks_repo_name: mmcgurn/MattFlowCases


text:
#include <gsl/gsl_sf_bessel.h>
#include "ccl.h"
#include "ctest.h"
#include "ccl_params.h"
#include "ccl_core.h"

CTEST(spacing_tests, linear_spacing_simple) {
  double * m = ccl_linear_spacing(0.0, 1.0, 5);
  ASSERT_NOT_NULL(m);
  ASSERT_DBL_NEAR_TOL(0.0, m[0], 1e-10);
  ASSERT_DBL_NEAR_TOL(0.25, m[1], 1e-10);
  ASSERT_DBL_NEAR_TOL(0.5, m[2], 1e-10);
  ASSERT_DBL_NEAR_TOL(0.75, m[3], 1e-10);
  ASSERT_DBL_NEAR_TOL(1.0, m[4], 1e-10);
  free(m);
}

CTEST(spacing_tests, log_spacing_simple) {
  double * m = ccl_log_spacing(1.0, 16.0, 5);
  ASSERT_NOT_NULL(m);
  ASSERT_DBL_NEAR_TOL(1.0, m[0], 1e-10);
  ASSERT_DBL_NEAR_TOL(2.0, m[1], 1e-10);
  ASSERT_DBL_NEAR_TOL(4.0, m[2], 1e-10);
  ASSERT_DBL_NEAR_TOL(8.0, m[3], 1e-10);
  ASSERT_DBL_NEAR_TOL(16.0, m[4], 1e-10);
  free(m);
}

CTEST(spacing_tests, linlog_spacing_simple) {
  double * m = ccl_linlog_spacing(1.0, 16.0, 32.0, 5, 9);
  ASSERT_NOT_NULL(m);
  ASSERT_DBL_NEAR_TOL(1.0, m[0], 1e-10);
  ASSERT_DBL_NEAR_TOL(2.0, m[1], 1e-10);
  ASSERT_DBL_NEAR_TOL(4.0, m[2], 1e-10);
  ASSERT_DBL_NEAR_TOL(8.0, m[3], 1e-10);
  ASSERT_DBL_NEAR_TOL(16.0, m[4], 1e-10);
  ASSERT_DBL_NEAR_TOL(18.0, m[5], 1e-10);
  ASSERT_DBL_NEAR_TOL(20.0, m[6], 1e-10);
  ASSERT_DBL_NEAR_TOL(22.0, m[7], 1e-10);
  ASSERT_DBL_NEAR_TOL(24.0, m[8], 1e-10);
  ASSERT_DBL_NEAR_TOL(26.0, m[9], 1e-10);
  ASSERT_DBL_NEAR_TOL(28.0, m[10], 1e-10);
  ASSERT_DBL_NEAR_TOL(30.0, m[11], 1e-10);
  ASSERT_DBL_NEAR_TOL(32.0, m[12], 1e-10);
  free(m);
}

CTEST(spherical_bessel_tests, compare_gsl) {
  int l, i;
  double xmin = 0.0;
  double xmax = 10.0;
  int Nx = 10000;
  double x, dx;

  dx = (xmax - xmin) / (Nx - 1);

  for (l=0; l < 15; ++l) {
    for (i=0; i < Nx; ++i) {
      x = dx * i + xmin;
      ASSERT_DBL_NEAR_TOL(
        ccl_j_bessel(l, x),
        gsl_sf_bessel_jl(l, x),
        1e-4);
    }
  }
}

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
/**
 * This file is part of the "libterminal" project
 *   Copyright (c) 2019-2020 Christian Parpart <christian@parpart.family>
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once

#include <terminal/GraphicsAttributes.h>
#include <terminal/primitives.h>

#include <crispy/Comparison.h>
#include <crispy/assert.h>

#include <gsl/span>
#include <gsl/span_ext>

#include <iterator>
#include <sstream>
#include <string>
#include <vector>

namespace terminal
{

enum class LineFlags : uint8_t
{
    None = 0x0000,
    Wrappable = 0x0001,
    Wrapped = 0x0002,
    Marked = 0x0004,
    // TODO: DoubleWidth  = 0x0010,
    // TODO: DoubleHeight = 0x0020,
};

template <typename, bool>
struct OptionalProperty;
template <typename T>
struct OptionalProperty<T, false>
{
};
template <typename T>
struct OptionalProperty<T, true>
{
    T value;
};

struct SimpleLineBuffer
{
    GraphicsAttributes attributes;
    std::string text;  // TODO: Try std::string_view later to avoid scattered copies.
    ColumnCount width; // page display width
};

template <typename Cell>
using InflatedLineBuffer = std::vector<Cell>;

template <typename Cell>
InflatedLineBuffer<Cell> inflate(SimpleLineBuffer const& input);

template <typename Cell>
using LineStorage = std::variant<SimpleLineBuffer, InflatedLineBuffer<Cell>>;

/**
 * Line<Cell> API.
 *
 * TODO: Use custom allocator for ensuring cache locality of Cells to sibling lines.
 * TODO: Make the line optimization work.
 */
template <typename Cell, bool Optimize = false>
class Line
{
  public:
    Line() = default;
    Line(Line const&) = default;
    Line(Line&&) noexcept = default;
    Line& operator=(Line const&) = default;
    Line& operator=(Line&&) noexcept = default;

    using InflatedBuffer = InflatedLineBuffer<Cell>;
    using Storage = LineStorage<Cell>;
    using value_type = Cell;
    using iterator = typename InflatedBuffer::iterator;
    using reverse_iterator = typename InflatedBuffer::reverse_iterator;
    using const_iterator = typename InflatedBuffer::const_iterator;

    Line(ColumnCount _width, LineFlags _flags, Cell _template = {}):
        buffer_(_width.as<size_t>(), _template /*, _allocator*/), flags_ { static_cast<unsigned>(_flags) }
    {
    }

    Line(ColumnCount _width, InflatedBuffer _buffer, LineFlags _flags):
        buffer_ { std::move(_buffer) }, flags_ { static_cast<unsigned>(_flags) }
    {
        buffer_.resize(unbox<size_t>(_width));
    }

    Line(InflatedBuffer _buffer, LineFlags _flags):
        buffer_ { std::move(_buffer) }, flags_ { static_cast<unsigned>(_flags) }
    {
    }

    constexpr static inline bool ColumnOptimized = Optimize;

    // This is experimental (aka. buggy) and going to be replaced with another optimization idea soon.
    //#define LINE_AVOID_CELL_RESET 1

    void reset(LineFlags _flags,
               GraphicsAttributes _attributes) noexcept // TODO: optimize by having no need to O(n) iterate
                                                        // through all buffer cells.
    {
        flags_ = static_cast<unsigned>(_flags);

        // if constexpr (ColumnOptimized)
        // {
        //     #if !defined(LINE_AVOID_CELL_RESET)
        //     if (buffer_.back().backgroundColor() != _attributes.backgroundColor)
        //         // TODO: also styles and UL color
        //         markUsedFirst(ColumnCount::cast_from(buffer_.size()));
        //
        //     for (auto i = 0; i < *columnsUsed(); ++i)
        //         buffer_[i].reset(_attributes);
        //     #endif
        //     markUsedFirst(ColumnCount(0));
        // }
        // else
        {
            for (Cell& cell: buffer_)
                cell.reset(_attributes);
        }
    }

    void markUsedFirst(ColumnCount /*_n*/) noexcept
    {
        // if constexpr (ColumnOptimized)
        //     usedColumns_.value = _n;
    }

    void reset(LineFlags _flags,
               GraphicsAttributes const& _attributes,
               char32_t _codepoint,
               uint8_t _width) noexcept
    {
        flags_ = static_cast<unsigned>(_flags);
        markUsedFirst(size());
        for (Cell& cell: buffer_)
        {
            cell.reset();
            cell.write(_attributes, _codepoint, _width);
        }
    }

    /**
     * Fills this line with the given content.
     *
     * @p _start offset into this line of the first charater
     * @p _sgr graphics rendition for the line starting at @c _start until the end
     * @p _ascii the US-ASCII characters to fill with
     */
    void fill(ColumnOffset _start, GraphicsAttributes const& _sgr, std::string_view _ascii)
    {
        auto& buffer = editable();

        assert(unbox<size_t>(_start) + _ascii.size() <= buffer.size());

        auto constexpr ASCII_Width = 1;
        auto const* s = _ascii.data();

        Cell* i = &buffer_[unbox<size_t>(_start)];
        Cell* e = i + _ascii.size();
        while (i != e)
            (i++)->write(_sgr, static_cast<char32_t>(*s++), ASCII_Width);

        // if constexpr (ColumnOptimized)
        // {
        //     #if !defined(LINE_AVOID_CELL_RESET)
        //     auto const e2 = buffer.data() + unbox<long>(columnsUsed());
        //     while (i != e2)
        //         (i++)->reset();
        //     #endif
        //     usedColumns_.value = boxed_cast<ColumnCount>(_start) + ColumnCount::cast_from(_ascii.size());
        // }
        // else
        {
            auto const e2 = buffer.data() + buffer.size();
            while (i != e2)
                (i++)->reset();
        }
    }

    ColumnCount size() const noexcept { return ColumnCount::cast_from(buffer_.size()); }

    ColumnCount columnsUsed() const noexcept
    {
        // if constexpr (ColumnOptimized)
        //     return usedColumns_.value;
        // else
        return size();
    }

    void resize(ColumnCount _count);

    gsl::span<Cell const> trim_blank_right() const noexcept;

    gsl::span<Cell const> cells() const noexcept { return buffer_; }

    gsl::span<Cell> useRange(ColumnOffset _start, ColumnCount _count) noexcept
    {
        markUsedFirst(std::max(columnsUsed(), boxed_cast<ColumnCount>(_start) + _count));
#if defined(__FreeBSD__)
        auto const bufferSpan = gsl::span(buffer_);
        return bufferSpan.subspan(unbox<size_t>(_start), unbox<size_t>(_count));
#else
        // NOTE: On FreeBSD the line below does not compile.
        return gsl::span(buffer_).subspan(unbox<size_t>(_start), unbox<size_t>(_count));
#endif
    }

    iterator begin() noexcept { return buffer_.begin(); }
    iterator end() noexcept { return std::next(std::begin(buffer_), unbox<int>(columnsUsed())); }

    const_iterator begin() const noexcept { return buffer_.begin(); }
    const_iterator end() const noexcept { return std::next(buffer_.begin(), unbox<int>(columnsUsed())); }

    reverse_iterator rbegin() noexcept { return buffer_.rbegin(); }
    reverse_iterator rend() noexcept { return buffer_.rend(); }

    Cell& front() noexcept { return buffer_.front(); }
    Cell const& front() const noexcept { return buffer_.front(); }

    Cell& back() noexcept { return *std::next(buffer_.begin(), unbox<int>(columnsUsed() - 1)); }
    Cell const& back() const noexcept { return *std::next(buffer_.begin(), unbox<int>(columnsUsed() - 1)); }

    Cell& useCellAt(ColumnOffset _column) noexcept
    {
        Require(ColumnOffset(0) <= _column);
        Require(_column <= ColumnOffset::cast_from(buffer_.size())); // Allow off-by-one for sentinel.
        return editable()[unbox<size_t>(_column)];
    }

    Cell const& at(ColumnOffset _column) const noexcept
    {
        Require(ColumnOffset(0) <= _column);
        Require(_column <= ColumnOffset::cast_from(buffer_.size())); // Allow off-by-one for sentinel.
        return buffer_[unbox<size_t>(_column)];
    }

    LineFlags flags() const noexcept { return static_cast<LineFlags>(flags_); }

    bool marked() const noexcept { return isFlagEnabled(LineFlags::Marked); }
    void setMarked(bool _enable) { setFlag(LineFlags::Marked, _enable); }

    bool wrapped() const noexcept { return isFlagEnabled(LineFlags::Wrapped); }
    void setWrapped(bool _enable) { setFlag(LineFlags::Wrapped, _enable); }

    bool wrappable() const noexcept { return isFlagEnabled(LineFlags::Wrappable); }
    void setWrappable(bool _enable) { setFlag(LineFlags::Wrappable, _enable); }

    LineFlags wrappableFlag() const noexcept { return wrappable() ? LineFlags::Wrappable : LineFlags::None; }
    LineFlags wrappedFlag() const noexcept { return marked() ? LineFlags::Wrapped : LineFlags::None; }
    LineFlags markedFlag() const noexcept { return marked() ? LineFlags::Marked : LineFlags::None; }

    LineFlags inheritableFlags() const noexcept
    {
        auto constexpr Inheritables = unsigned(LineFlags::Wrappable) | unsigned(LineFlags::Marked);
        return static_cast<LineFlags>(flags_ & Inheritables);
    }

    void setFlag(LineFlags _flag, bool _enable) noexcept
    {
        if (_enable)
            flags_ |= static_cast<unsigned>(_flag);
        else
            flags_ &= ~static_cast<unsigned>(_flag);
    }

    bool isFlagEnabled(LineFlags _flag) const noexcept
    {
        return (flags_ & static_cast<unsigned>(_flag)) != 0;
    }

    InflatedBuffer reflow(ColumnCount _newColumnCount);
    std::string toUtf8() const;
    std::string toUtf8Trimmed() const;

    // Returns a reference to this mutable grid-line buffer.
    //
    // If this line has been stored in an optimized state, then
    // the line will be first unpacked into a vector of grid cells.
    InflatedBuffer& editable();

  private:
    InflatedBuffer buffer_;
    Storage storage_;
    unsigned flags_ = 0;
    // OptionalProperty<ColumnCount, ColumnOptimized> usedColumns_;
};

constexpr LineFlags operator|(LineFlags a, LineFlags b) noexcept
{
    return LineFlags(unsigned(a) | unsigned(b));
}

constexpr LineFlags operator~(LineFlags a) noexcept
{
    return LineFlags(~unsigned(a));
}

constexpr LineFlags operator&(LineFlags a, LineFlags b) noexcept
{
    return LineFlags(unsigned(a) & unsigned(b));
}

template <typename Cell, bool Optimize>
inline typename Line<Cell, Optimize>::InflatedBuffer& Line<Cell, Optimize>::editable()
{
    // TODO: when we impement the line text buffer optimization,
    // then this is the place where we want to *promote* a possibly
    // optimized text buffer to a full grid line buffer.
#if 0
    if (std::holds_alternative<SimpleLineBuffer>(storage_))
        storage_ = inflate<Cell>(std::get<SimpleLineBuffer>(storage_));

    return std::get<InflatedBuffer>(storage_);
#else
    return buffer_;
#endif
}

} // namespace terminal

namespace fmt // {{{
{
template <>
struct formatter<terminal::LineFlags>
{
    template <typename ParseContext>
    auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(const terminal::LineFlags _flags, FormatContext& ctx) const
    {
        static const std::array<std::pair<terminal::LineFlags, std::string_view>, 3> nameMap = {
            std::pair { terminal::LineFlags::Wrappable, std::string_view("Wrappable") },
            std::pair { terminal::LineFlags::Wrapped, std::string_view("Wrapped") },
            std::pair { terminal::LineFlags::Marked, std::string_view("Marked") },
        };
        std::string s;
        for (auto const& mapping: nameMap)
        {
            if ((mapping.first & _flags) != terminal::LineFlags::None)
            {
                if (!s.empty())
                    s += ",";
                s += mapping.second;
            }
        }
        return format_to(ctx.out(), s);
    }
};
} // namespace fmt

max_stars_repo_name: sebastianrakel/contour

max_issues_repo_name: sebastianrakel/contour

max_forks_repo_name: sebastianrakel/contour


text:
#ifndef _DSGLD_LDA_MODEL_H__
#define _DSGLD_LDA_MODEL_H__

#include <El.hpp>
#include <gsl/gsl_rng.h>

#include "sgld_model.h"

using std::string;
using std::vector;

namespace dsgld {

class LDAModel : public SGLDModel<double, int> {
 public:
   LDAModel(
       const El::Matrix<int>& X,
       const int K,
       const double alpha,
       const double beta);

  ~LDAModel() {};

  El::Matrix<double> sgldEstimate(const El::Matrix<double>& theta) override;

  El::Matrix<double> nablaLogPrior(const El::Matrix<double>& theta) const override;

  void writePerplexities(const string& filename);

  int NumGibbsSteps() const;
  LDAModel* NumGibbsSteps(const int);

 protected:
  void gibbsSample(
      const El::Matrix<int>& doc,
      const El::Matrix<double>& theta,
      El::Matrix<int>& index_to_topic,
      El::Matrix<int>& topic_counts) const;

  double estimatePerplexity(
      const El::Matrix<double>& theta,
      const El::Matrix<double>& theta_sum_over_w,
      const int num_words_in_doc,
      const El::Matrix<int>& topic_counts) const;

 private:
  const double alpha_;
  const double beta_;
  const int W;
  const int K;
  int numGibbsSteps_;
  vector<double> perplexities_;
  const gsl_rng* rng;
};

}  // namespace dsgld

#endif  // _DSGLD_LDA_MODEL_H__

max_stars_repo_name: feynmanliang/travelling-cluster-chain

max_issues_repo_name: feynmanliang/travelling-cluster-chain

max_forks_repo_name: feynmanliang/travelling-cluster-chain


text:
/* spgetset.c
 *
 * Copyright (C) 2012 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spmatrix.h>

#include "avl.c"

static void *tree_find(const gsl_spmatrix *m, const size_t i, const size_t j);

double
gsl_spmatrix_get(const gsl_spmatrix *m, const size_t i, const size_t j)
{
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0.0);
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0.0);
    }
  else if (m->nz == 0)
    {
      /* no non-zero elements added to matrix */
      return (0.0);
    }
  else
    {
      if (GSL_SPMATRIX_ISTRIPLET(m))
        {
          /* traverse binary tree to search for (i,j) element */
          void *ptr = tree_find(m, i, j);
          double x = ptr ? *(double *) ptr : 0.0;

          return x;
        }
      else if (GSL_SPMATRIX_ISCCS(m))
        {
          const size_t *mi = m->i;
          const size_t *mp = m->p;
          size_t p;

          /* loop over column j and search for row index i */
          for (p = mp[j]; p < mp[j + 1]; ++p)
            {
              if (mi[p] == i)
                return m->data[p];
            }
        }
      else if (GSL_SPMATRIX_ISCRS(m))
        {
          const size_t *mj = m->i;
          const size_t *mp = m->p;
          size_t p;

          /* loop over row i and search for column index j */
          for (p = mp[i]; p < mp[i + 1]; ++p)
            {
              if (mj[p] == j)
                return m->data[p];
            }
        }
      else
        {
          GSL_ERROR_VAL("unknown sparse matrix type", GSL_EINVAL, 0.0);
        }

      /* element not found; return 0 */
      return 0.0;
    }
} /* gsl_spmatrix_get() */

/*
gsl_spmatrix_set()
  Add an element to a matrix in triplet form

Inputs: m - spmatrix
        i - row index
        j - column index
        x - matrix value
*/

int
gsl_spmatrix_set(gsl_spmatrix *m, const size_t i, const size_t j,
                 const double x)
{
  if (!GSL_SPMATRIX_ISTRIPLET(m))
    {
      GSL_ERROR("matrix not in triplet representation", GSL_EINVAL);
    }
  else if (x == 0.0)
    {
      /* traverse binary tree to search for (i,j) element */
      void *ptr = tree_find(m, i, j);

      /*
       * just set the data element to 0; it would be easy to
       * delete the node from the tree with avl_delete(), but
       * we'd also have to delete it from the data arrays which
       * is less simple
       */
      if (ptr != NULL)
        *(double *) ptr = 0.0;

      return GSL_SUCCESS;
    }
  else
    {
      int s = GSL_SUCCESS;
      void *ptr;

      /* check if matrix needs to be realloced */
      if (m->nz >= m->nzmax)
        {
          s = gsl_spmatrix_realloc(2 * m->nzmax, m);
          if (s)
            return s;
        }

      /* store the triplet (i, j, x) */
      m->i[m->nz] = i;
      m->p[m->nz] = j;
      m->data[m->nz] = x;

      ptr = avl_insert(m->tree_data->tree, &m->data[m->nz]);
      if (ptr != NULL)
        {
          /* found duplicate entry (i,j), replace with new x */
          *((double *) ptr) = x;
        }
      else
        {
          /* no duplicate (i,j) found, update indices as needed */

          /* increase matrix dimensions if needed */
          m->size1 = GSL_MAX(m->size1, i + 1);
          m->size2 = GSL_MAX(m->size2, j + 1);

          ++(m->nz);
        }

      return s;
    }
} /* gsl_spmatrix_set() */

double *
gsl_spmatrix_ptr(gsl_spmatrix *m, const size_t i, const size_t j)
{
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL);
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL);
    }
  else
    {
      if (GSL_SPMATRIX_ISTRIPLET(m))
        {
          /* traverse binary tree to search for (i,j) element */
          void *ptr = tree_find(m, i, j);
          return (double *) ptr;
        }
      else if (GSL_SPMATRIX_ISCCS(m))
        {
          const size_t *mi = m->i;
          const size_t *mp = m->p;
          size_t p;

          /* loop over column j and search for row index i */
          for (p = mp[j]; p < mp[j + 1]; ++p)
            {
              if (mi[p] == i)
                return &(m->data[p]);
            }
        }
      else if (GSL_SPMATRIX_ISCRS(m))
        {
          const size_t *mj = m->i;
          const size_t *mp = m->p;
          size_t p;

          /* loop over row i and search for column index j */
          for (p = mp[i]; p < mp[i + 1]; ++p)
            {
              if (mj[p] == j)
                return &(m->data[p]);
            }
        }
      else
        {
          GSL_ERROR_NULL("unknown sparse matrix type", GSL_EINVAL);
        }

      /* element not found; return 0 */
      return NULL;
    }
} /* gsl_spmatrix_ptr() */

/*
tree_find()
  Find node in tree corresponding to matrix entry (i,j). Adapted
from avl_find()

Inputs: m - spmatrix
        i - row index
        j - column index

Return: pointer to tree node data if found, NULL if not found
*/

static void *
tree_find(const gsl_spmatrix *m, const size_t i, const size_t j)
{
  const struct avl_table *tree = (struct avl_table *) m->tree_data->tree;
  const struct avl_node *p;

  for (p = tree->avl_root; p != NULL; )
    {
      size_t n = (double *) p->avl_data - m->data;
      size_t pi = m->i[n];
      size_t pj = m->p[n];
      int cmp = gsl_spmatrix_compare_idx(i, j, pi, pj);

      if (cmp < 0)
        p = p->avl_link[0];
      else if (cmp > 0)
        p = p->avl_link[1];
      else /* |cmp == 0| */
        return p->avl_data;
    }

  return NULL;
} /* tree_find() */

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl/span>

#include "aas/aas_math.h"

namespace aas {

	struct AnimPlayerStream;
	class AnimatedModel;
	struct SkelAnim;
	struct SkelAnimEvent;
	struct SkelBoneState;
	class IAnimEventHandler;
	class Skeleton;
	struct AnimEvent;

	class AnimPlayer {
	public:
		AnimatedModel * ownerAnim;
		char field_D;
		char field_E;
		char field_F;
		float weight; // 1.0 = Fully weighted, otherwise weighted with primary anims
		float fadingSpeed; // Velocity with which weight is being changed
		int eventHandlingDepth;
		AnimPlayer * nextRunningAnim;
		AnimPlayer * prevRunningAnim;
		const SkelAnim * animation;
		int streamCount;
		AnimPlayerStream * streams[4];
		float streamFps[4]; // "frames per drive unit"
		int8_t streamVariationIndices[4]; // indices into streamVariationIds
		int8_t streamVariationIds[4];
		int8_t variationId[4];
		float currentTimeEvents;
		float currentTime;
		float duration;
		float frameRate;
		float distancePerSecond;
		std::vector<AnimEvent> events;
		IAnimEventHandler* eventHandler;
	public:

		AnimPlayer();
		~AnimPlayer();

		void GetDistPerSec(float* distPerSec);
		void GetRotationPerSec(float* rotationPerSec);
		void AdvanceEvents(float timeChanged, float distanceChanged, float rotationChanged);
		void FadeInOrOut(float timeChanged);
		void method6(gsl::span<SkelBoneState> boneStateOut, float timeChanged, float distanceChanged, float rotationChanged);
		void SetTime(float time);
		float GetCurrentFrame();
		float method9();
		float method10();
		void EnterEventHandling();
		void LeaveEventHandling();
		int GetEventHandlingDepth();

		void Attach(AnimatedModel *owner, int animIdx, IAnimEventHandler *eventHandler);
		void Setup2(float a);

	private:
		void SetEvents(const AnimatedModel *owner, const SkelAnim *anim);

	};

	AnimPlayerStream * CreateAnimPlayerStream(Skeleton* skeleton);

}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#include "chmath.h"
#include "err.h"

static char* me = "chmath";

static double get(const gsl_matrix* m, int i, int j) {
  return gsl_matrix_get(m, i, j);
}

int chmath_eig_values(const double A[6], /**/ double VAL[3]) {
  enum { XX, XY, XZ, YY, YZ, ZZ };
  enum { YX = XY, ZX = XZ, ZY = YZ };
  enum { X, Y, Z };

  double B[3 * 3];
  gsl_matrix_view m;
  gsl_vector* val;
  gsl_matrix* vec;
  gsl_eigen_symmv_workspace* w;
  int i, status;

  val = gsl_vector_alloc(3);
  vec = gsl_matrix_alloc(3, 3);
  w = gsl_eigen_symmv_alloc(3);

  i = 0;
  B[i++] = A[XX];
  B[i++] = A[XY];
  B[i++] = A[XZ];
  B[i++] = A[YX];
  B[i++] = A[YY];
  B[i++] = A[YZ];
  B[i++] = A[ZX];
  B[i++] = A[ZY];
  B[i++] = A[ZZ];

  m = gsl_matrix_view_array(B, 3, 3);
  status = gsl_eigen_symmv(&m.matrix, val, vec, w);
  if (status != GSL_SUCCESS) ERR(("gsl_eigen_symmv failed"));
  gsl_eigen_symmv_sort(val, vec, GSL_EIGEN_SORT_ABS_ASC);

  i = 0;
  VAL[i] = gsl_vector_get(val, i);
  i++;
  VAL[i] = gsl_vector_get(val, i);
  i++;
  VAL[i] = gsl_vector_get(val, i);

  gsl_vector_free(val);
  gsl_matrix_free(vec);
  gsl_eigen_symmv_free(w);
  return 0;
}

int chmath_eig_vectors(
    const double A[6], /**/ double a[3], double b[3], double c[3]) {
  enum { XX, XY, XZ, YY, YZ, ZZ };
  enum { YX = XY, ZX = XZ, ZY = YZ };
  enum { X, Y, Z };

  double B[3 * 3];
  gsl_matrix_view m;
  gsl_vector* val;
  gsl_matrix* vec;
  gsl_eigen_symmv_workspace* w;
  int i, status;

  val = gsl_vector_alloc(3);
  vec = gsl_matrix_alloc(3, 3);
  w = gsl_eigen_symmv_alloc(3);

  i = 0;
  B[i++] = A[XX];
  B[i++] = A[XY];
  B[i++] = A[XZ];
  B[i++] = A[YX];
  B[i++] = A[YY];
  B[i++] = A[YZ];
  B[i++] = A[ZX];
  B[i++] = A[ZY];
  B[i++] = A[ZZ];

  m = gsl_matrix_view_array(B, 3, 3);
  status = gsl_eigen_symmv(&m.matrix, val, vec, w);
  if (status != GSL_SUCCESS) ERR(("gsl_eigen_symmv failed"));
  gsl_eigen_symmv_sort(val, vec, GSL_EIGEN_SORT_ABS_ASC);

  i = 0;
  a[X] = get(vec, X, X);
  a[Y] = get(vec, Y, X);
  a[Z] = get(vec, Z, X);
  b[X] = get(vec, X, Y);
  b[Y] = get(vec, Y, Y);
  b[Z] = get(vec, Z, Y);
  c[X] = get(vec, X, Z);
  c[Y] = get(vec, Y, Z);
  c[Z] = get(vec, Z, Z);
  gsl_vector_free(val);
  gsl_matrix_free(vec);
  gsl_eigen_symmv_free(w);
  return 0;
}


text:
#include <math.h>
#include <gmp.h>
#if BLAS_IMPL == openblas
#include <blas/cblas.h>
#else
#include <cblas.h>
#endif
#include "attrib.h"
#include "math.h"

#define CLEANUP \
    return err; \
    error: \
    if (ctx->err == 0) \
    	strcpy(ctx->err_msg, "Unexpected math."); \
    ctx->err = err; \
    return err;

#define CHECK(expr) \
    if ((err = (expr))) goto error;

VIO_CONST
uint32_t vmaxui(uint32_t a, uint32_t b) {
    return a > b ? a : b;
}

#define OP_IMPL(name) \
    vio_err_t generic_##name(vio_ctx *ctx, vio_val *a, vio_val *b) { \
        vio_err_t err = 0; \
        vio_val *va, *vb; \
        if (a->what == vv_str || b->what == vv_str || a->what == vv_quot || b->what == vv_quot || \
            a->what == vv_tagword || b->what == vv_tagword) { \
            err = vio_raise(ctx, VE_WRONG_TYPE, \
                            "Function '" #name "' expects numeric types or vectors; received '%s' and '%s' operands.", \
                            vio_val_type_name(a->what), vio_val_type_name(b->what));\
            goto error; \
        }

#define END_OP(_name) \
    CLEANUP \
}

#define GIVEN_VECF \
    if (a->what == vv_vecf && b->what == vv_vecf && a->vlen == b->vlen) {

#define DONE_VECF }

#define GIVEN_NUMERIC \
    else if (vio_is_numeric(a) && vio_is_numeric(b)) { \
        CHECK(vio_coerce(ctx, b, &vb, a->what)) \
        if (vb == NULL) { \
            CHECK(vio_coerce(ctx, a, &va, b->what)) \
            if (va == NULL) { \
                err = VE_NUMERIC_CONVERSION_FAIL; \
                goto error; \
            } \
            a = b; \
            b = va; \
        } \
        else \
            b = vb; \
        switch (a->what) {

#define DONE_NUMERIC \
        default: err = VE_WRONG_TYPE; goto error; \
        } \
    }

#define STANDARD_NUMERIC(op, fop) \
    GIVEN_NUMERIC \
        case vv_int: CHECK(vio_push_int(ctx, b->i32 op a->i32)) break; \
        case vv_float: CHECK(vio_push_float(ctx, b->f32 op a->f32)) break; \
        case vv_num: { mpf_t c; mpf_init(c); mpf_##fop(c, b->n, a->n); CHECK(vio_push_num(ctx, c)) break; } \
    DONE_NUMERIC

#define GIVEN_VEC \
    else CHECK(vio_coerce(ctx, a, &va, vv_vec)) \
    else CHECK(vio_coerce(ctx, b, &vb, vv_vec)) \
    else if (va && vb) { \

#define DONE_VEC }

#define AUTO_VECTORIZE(name) \
    GIVEN_VEC \
        uint32_t maxl = vmaxui(va->vlen, vb->vlen); \
        for (uint32_t i = 0; i < maxl; ++i) \
            CHECK(generic_##name(ctx, va->vv[i % va->vlen], vb->vv[i % vb->vlen])) \
        CHECK(vio_push_vec(ctx, maxl)) \
    DONE_VEC

OP_IMPL(add)
    GIVEN_VECF
        cblas_daxpy(b->vlen, 1, a->vf32, 1, b->vf32, 1);
        CHECK(vio_push_vecf32(ctx, b->vlen, b->vf32))
    DONE_VECF
    STANDARD_NUMERIC(+, add)
    AUTO_VECTORIZE(add)
END_OP(add)

OP_IMPL(sub)
    GIVEN_VECF
        cblas_daxpy(b->vlen, -1.0, a->vf32, 1, b->vf32, 1);
        CHECK(vio_push_vecf32(ctx, b->vlen, b->vf32))
    DONE_VECF
    STANDARD_NUMERIC(-, sub)
    AUTO_VECTORIZE(sub)
END_OP(sub)

OP_IMPL(mul)
    if (b->what == vv_vecf && a->what == vv_float) {
        va = a;
        a = b;
        b = va;
    }
    if (a->what == vv_vecf && b->what == vv_float) {
        cblas_dscal(a->vlen, b->f32, a->vf32, 1);
        CHECK(vio_push_vecf32(ctx, a->vlen, a->vf32))
    }
    else GIVEN_VECF
        for (uint32_t i = 0; i < a->vlen; ++i)
            a->vf32[i] *= b->vf32[i];
        CHECK(vio_push_vecf32(ctx, a->vlen, a->vf32))
    DONE_VECF
    STANDARD_NUMERIC(*, mul)
    AUTO_VECTORIZE(mul)
END_OP(mul)

OP_IMPL(div)
    if (b->what == vv_vecf && a->what == vv_float) {
        va = a;
        a = b;
        b = va;
    }
    if (a->what == vv_vecf && b->what == vv_float) {
        cblas_dscal(a->vlen, 1.0 / b->f32, a->vf32, 1);
        CHECK(vio_push_vecf32(ctx, a->vlen, a->vf32))
    }
    else GIVEN_VECF
        for (uint32_t i = 0; i < a->vlen; ++i)
            a->vf32[i] /= b->vf32[i];
        CHECK(vio_push_vecf32(ctx, a->vlen, a->vf32))
    DONE_VECF
    STANDARD_NUMERIC(/, div)
    AUTO_VECTORIZE(div)
END_OP(div)

#define BIN_OP(name) \
    vio_err_t vio_##name(vio_ctx *ctx) { \
        vio_err_t err = 0; \
        VIO__ENSURE_ATLEAST(2) \
\
        vio_val *a = ctx->stack[--ctx->sp], *b; \
        b = ctx->stack[--ctx->sp]; \
        err = generic_##name(ctx, a, b); \
\
        CLEANUP \
    }

BIN_OP(add)
BIN_OP(sub)
BIN_OP(mul)
BIN_OP(div)

#define UN_IMPL(f) \
    vio_err_t vio_##f(vio_ctx *ctx) { \
        vio_err_t err = 0; \
        vio_float x, y; \
        VIO__CHECK(vio_pop_float(ctx, &x)); \
        y = f(x); \
        VIO__CHECK(vio_push_float(ctx, y)); \
        return 0; \
        error: \
        return err; \
    }

LIST_MATH_UNARY(UN_IMPL)