<commit_before>//-----------------------------------------------------------------------bl-
//--------------------------------------------------------------------------
// 
// GRINS - General Reacting Incompressible Navier-Stokes 
//
// Copyright (C) 2010-2012 The PECOS Development Team
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the Version 2 GNU General
// Public License as published by the Free Software Foundation.
//
// 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 library; if not, write to the Free Software
// Foundation, Inc. 51 Franklin Street, Fifth Floor, Boston, MA
// 02110-1301 USA
//
//-----------------------------------------------------------------------el-
//
// $Id$
//
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------

#include "reacting_low_mach_navier_stokes.h"

namespace GRINS
{
  template<class Mixture>
  ReactingLowMachNavierStokes<Mixture>::ReactingLowMachNavierStokes(const PhysicsName& physics_name, const GetPot& input)
    : ReactingLowMachNavierStokesBase<Mixture>(physics_name,input),
      _p_pinning(input,physics_name)
  {
    this->read_input_options(input);

    return;
  }

  template<class Mixture>
  ReactingLowMachNavierStokes<Mixture>::~ReactingLowMachNavierStokes()
  {
    return;
  }

  template<class Mixture>
  void ReactingLowMachNavierStokes<Mixture>::read_input_options( const GetPot& input )
  {
    // Other quantities read in base class

    // Read pressure pinning information
    this->_pin_pressure = input("Physics/"+reacting_low_mach_navier_stokes+"/pin_pressure", false );
  
    return;
  }

  template<class Mixture>
  void ReactingLowMachNavierStokes<Mixture>::init_context( libMesh::DiffContext &context )
  {
    // First call base class
    GRINS::ReactingLowMachNavierStokesBase<Mixture>::init_context(context);

    libMesh::FEMContext &c = libmesh_cast_ref<libMesh::FEMContext&>(context);

    // We also need the side shape functions, etc.
    c.side_fe_var[this->_u_var]->get_JxW();
    c.side_fe_var[this->_u_var]->get_phi();
    c.side_fe_var[this->_u_var]->get_dphi();
    c.side_fe_var[this->_u_var]->get_xyz();

    c.side_fe_var[this->_T_var]->get_JxW();
    c.side_fe_var[this->_T_var]->get_phi();
    c.side_fe_var[this->_T_var]->get_dphi();
    c.side_fe_var[this->_T_var]->get_xyz();

    return;
  }

  template<class Mixture>
  bool ReactingLowMachNavierStokes<Mixture>::element_time_derivative( bool request_jacobian,
								      libMesh::DiffContext& context,
								      libMesh::FEMSystem* system )
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
    
    unsigned int n_qpoints = c.element_qrule->n_points();

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
	// Compute species mass fractions at quadrature points
	std::vector<Real> Y;
	Y.reserve(this->_n_species);
	for( unsigned int s = 0; s < this->_n_species; s++ )
	  {
	    Y.push_back( c.interior_value(this->_species_vars[s],qp) );
	  }

	// Build up cache at element interior quadrature point
	ReactingFlowCache cache( c.interior_value(this->_T_var, qp), 
				 this->get_p0_steady(c,qp), Y );

	this->build_reacting_flow_cache(c, cache, qp);

	this->assemble_mass_time_deriv(c, cache, qp);
	this->assemble_species_time_deriv(c, cache, qp);
	this->assemble_momentum_time_deriv(c, cache, qp);
	this->assemble_energy_time_deriv(c, cache, qp);
      }

    // Pin p = p_value at p_point
    if( this->_pin_pressure )
      {
	this->_p_pinning.pin_value( context, request_jacobian, this->_p_var );
      }

    return request_jacobian;
  }

  template<class Mixture>
  void ReactingLowMachNavierStokes<Mixture>::assemble_mass_time_deriv(libMesh::FEMContext& c, 
								      const ReactingFlowCache& cache, 
								      unsigned int qp)
  {
    // The number of local degrees of freedom in each variable.
    const unsigned int n_p_dofs = c.dof_indices_var[this->_p_var].size();
    
    // Element Jacobian * quadrature weights for interior integration.
    const std::vector<libMesh::Real> &JxW =
      c.element_fe_var[this->_u_var]->get_JxW();
    
    // The pressure shape functions at interior quadrature points.
    const std::vector<std::vector<libMesh::Real> >& p_phi =
      c.element_fe_var[this->_p_var]->get_phi();
    
    libMesh::DenseSubVector<Number> &Fp = *c.elem_subresiduals[this->_p_var]; // R_{p}
    
    libMesh::Number T = cache.T();
    const libMesh::NumberVectorValue& U = cache.U();
    libMesh::Number divU = cache.divU();
    const libMesh::Gradient& grad_T = cache.grad_T();
    
    for (unsigned int i=0; i != n_p_dofs; i++)
      {
	Fp(i) += MOLAR_MASS_TERM (-U*grad_T/T + divU)*p_phi[i][qp]*JxW[qp];<commit_msg>[GRINS]: Checkpointing before meeting (#2619).<commit_after>//-----------------------------------------------------------------------bl-
//--------------------------------------------------------------------------
// 
// GRINS - General Reacting Incompressible Navier-Stokes 
//
// Copyright (C) 2010-2012 The PECOS Development Team
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the Version 2 GNU General
// Public License as published by the Free Software Foundation.
//
// 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 library; if not, write to the Free Software
// Foundation, Inc. 51 Franklin Street, Fifth Floor, Boston, MA
// 02110-1301 USA
//
//-----------------------------------------------------------------------el-
//
// $Id$
//
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------

#include "reacting_low_mach_navier_stokes.h"

namespace GRINS
{
  template<class Mixture>
  ReactingLowMachNavierStokes<Mixture>::ReactingLowMachNavierStokes(const PhysicsName& physics_name, const GetPot& input)
    : ReactingLowMachNavierStokesBase<Mixture>(physics_name,input),
      _p_pinning(input,physics_name)
  {
    this->read_input_options(input);

    return;
  }

  template<class Mixture>
  ReactingLowMachNavierStokes<Mixture>::~ReactingLowMachNavierStokes()
  {
    return;
  }

  template<class Mixture>
  void ReactingLowMachNavierStokes<Mixture>::read_input_options( const GetPot& input )
  {
    // Other quantities read in base class

    // Read pressure pinning information
    this->_pin_pressure = input("Physics/"+reacting_low_mach_navier_stokes+"/pin_pressure", false );
  
    return;
  }

  template<class Mixture>
  void ReactingLowMachNavierStokes<Mixture>::init_context( libMesh::DiffContext &context )
  {
    // First call base class
    GRINS::ReactingLowMachNavierStokesBase<Mixture>::init_context(context);

    libMesh::FEMContext &c = libmesh_cast_ref<libMesh::FEMContext&>(context);

    // We also need the side shape functions, etc.
    c.side_fe_var[this->_u_var]->get_JxW();
    c.side_fe_var[this->_u_var]->get_phi();
    c.side_fe_var[this->_u_var]->get_dphi();
    c.side_fe_var[this->_u_var]->get_xyz();

    c.side_fe_var[this->_T_var]->get_JxW();
    c.side_fe_var[this->_T_var]->get_phi();
    c.side_fe_var[this->_T_var]->get_dphi();
    c.side_fe_var[this->_T_var]->get_xyz();

    return;
  }

  template<class Mixture>
  bool ReactingLowMachNavierStokes<Mixture>::element_time_derivative( bool request_jacobian,
								      libMesh::DiffContext& context,
								      libMesh::FEMSystem* system )
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
    
    unsigned int n_qpoints = c.element_qrule->n_points();

    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
	// Compute species mass fractions at quadrature points
	std::vector<Real> Y;
	Y.reserve(this->_n_species);
	for( unsigned int s = 0; s < this->_n_species; s++ )
	  {
	    Y.push_back( c.interior_value(this->_species_vars[s],qp) );
	  }

	// Build up cache at element interior quadrature point
	ReactingFlowCache cache( c.interior_value(this->_T_var, qp), 
				 this->get_p0_steady(c,qp), Y );

	this->build_reacting_flow_cache(c, cache, qp);

	this->assemble_mass_time_deriv(c, cache, qp);
	this->assemble_species_time_deriv(c, cache, qp);
	this->assemble_momentum_time_deriv(c, cache, qp);
	this->assemble_energy_time_deriv(c, cache, qp);
      }

    // Pin p = p_value at p_point
    if( this->_pin_pressure )
      {
	this->_p_pinning.pin_value( context, request_jacobian, this->_p_var );
      }

    return request_jacobian;
  }

  template<class Mixture>
  void ReactingLowMachNavierStokes<Mixture>::assemble_mass_time_deriv(libMesh::FEMContext& c, 
								      const ReactingFlowCache& cache, 
								      unsigned int qp)
  {
    // The number of local degrees of freedom in each variable.
    const unsigned int n_p_dofs = c.dof_indices_var[this->_p_var].size();
    
    // Element Jacobian * quadrature weights for interior integration.
    const std::vector<libMesh::Real> &JxW =
      c.element_fe_var[this->_u_var]->get_JxW();
    
    // The pressure shape functions at interior quadrature points.
    const std::vector<std::vector<libMesh::Real> >& p_phi =
      c.element_fe_var[this->_p_var]->get_phi();
    
    libMesh::DenseSubVector<Number> &Fp = *c.elem_subresiduals[this->_p_var]; // R_{p}
    
    libMesh::Number T = cache.T();
    const libMesh::NumberVectorValue& U = cache.U();
    libMesh::Number divU = cache.divU();
    const libMesh::Gradient& grad_T = cache.grad_T();
    Real M = cache.M();
    const std::vector<Real>& M_species = cache.M_species();
    const std::vector<libMesh::Gradient>& grad_w = cache.mass_fractions();

    libmesh_assert_equal_to( M_species.size(), this->_n_species );
    libmesh_assert_equal_to( grad_w.size(), this->_n_species );
    
    libMesh::Gradient mass_term(0.0,0.0,0.0);
    for(unsigned int s=0; s < this->_n_species; s++ )
      {
	mass_term += grad_w[s]/M_species[s];
      }
    mass_term *= M;
    
    for (unsigned int i=0; i != n_p_dofs; i++)
      {
	Fp(i) += (-U*(mass_term + grad_T/T) + divU)*p_phi[i][qp]*JxW[qp];
{'fraction_non_alphanumeric': 0.18606395948607335, 'fraction_numerical': 0.0055331520210072215, 'mean_word_length': 2.7896233120113716, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 126, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '365400', 'n_tokens_mistral': 3649, 'n_tokens_neox': 3419, 'n_words': 738}