Fullupwind about conservative like-advection

[wangzhaohao]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Dear MOOSE Developers.
I am solving this problem:
$\frac{\partial p}{\partial t}+\nabla•[(1-p)p\overrightarrow{v}-\nu\nabla p]=0$
in which $p$ is variable, $\overrightarrow{v}$ is velocity about temperature, $\nu$ is a diffusion coefficient.
I would like to creat Fulupwind kernel to compute $\nabla•(1-p)p\overrightarrow{v}$, and there is ConservativeAdvection and a disscusion.
To be more comprehensive, I paste all the code below

#pragma once
#include "Kernel.h"

/*
adcection-diffusion of pores
this is advection part
$\bigtriangledown \dot [(1 - p) p \vec{v}]$
*/

class PoreAdvection : public Kernel
{
public:
  PoreAdvection(const InputParameters & parameters);

  static InputParameters validParams();

protected:
  virtual Real computeQpResidual() override;
  virtual Real computeQpJacobian() override;
  virtual void computeResidual() override;
  virtual void computeJacobian() override;
  RealVectorValue velocity() const; // m/s

  const enum class UpwindType { none, full } _upwinding;
  enum class JacRes
  {
    CALCULATE_RESIDUAL = 0,
    CALCULATE_JACOBIAN = 1
  };
  const VariableValue & _u_nodal;
  std::vector<bool> _upwind_node;
  std::vector<Real> _dtotal_mass_out;

private:
  const VariableValue & _temperature;
  const VariableGradient & _temperature_grad;
  const VariableValue & _plutonium_concentration;
  Real negSpeedQp() const;
  void fullUpwind(JacRes res_or_jac);
};

#include "PoreAdvection.h"

registerMooseObject("MiracleApp", PoreAdvection);

InputParameters
PoreAdvection::validParams()
{
  InputParameters params = Kernel::validParams();
  params.addRequiredCoupledVar("temperature", "the MOX temperature in K");
  params.addCoupledVar("plutonium_concentration", "the plutonium concentration in range 0 - 1");
  MooseEnum upwinding("none full", "none");
  params.addParam<MooseEnum>("upwinding", upwinding, "the upwinding type");

  return params;
}

PoreAdvection::PoreAdvection(const InputParameters & parameters)
  : Kernel(parameters),
    _temperature(coupledValue("temperature")),
    _temperature_grad(coupledGradient("temperature")),
    _plutonium_concentration(coupledValue("plutonium_concentration")),
    _upwinding(getParam<MooseEnum>("upwinding").getEnum<UpwindType>()),
    _u_nodal(_var.dofValues()),
    _upwind_node(0),
    _dtotal_mass_out(0)
{
}

RealVectorValue
PoreAdvection::velocity() const
{
  RealVectorValue pore_grad = 0; //[K/m]
  if (_plutonium_concentration[_qp] == 0)
    pore_grad = 4 * _temperature_grad[_qp];
  else // MOX
    pore_grad = 5 * _temperature_grad[_qp];

  Real molecules_volume = 0;              //[m^3]
  if (_plutonium_concentration[_qp] == 0) // from THE KINETICS OF PORE MOVEMENT IN UOz FUEL RODS
  {
    molecules_volume =
        6.023e23 * 270e-3 /
        (1e3 * 10.97 /
         (0.998 + 6.395e-6 * _temperature[_qp] + 3.543e-9 * std::pow(_temperature[_qp], 2) +
          3e-12 * std::pow(_temperature[_qp], 3)));
  }
  else // MOX
  {
    Real puo2_mole_fraction =
        _plutonium_concentration[_qp] / 274. /
        (_plutonium_concentration[_qp] / 274. + (1 - _plutonium_concentration[_qp] / 270.));
    Real uo2_mole_fraction =
        (1 - _plutonium_concentration[_qp]) / 270. /
        (_plutonium_concentration[_qp] / 274. + (1 - _plutonium_concentration[_qp]) / 270.);
    molecules_volume =
        0.25 * std::pow(puo2_mole_fraction * 5.396e-10 + uo2_mole_fraction * 5.470e-10, 3);
  }

  Real fuel_number_density = 5e17;   //[m^-3]
  Real gas_number_density = 2.46e25; //[m^-3]
  Real equivalent_mass = 1 / ((1 / 28.e-3 + 1 / 270e-3) * 6.023e23);//[kg]
  Real diffusion_coeff =
      3. / 8. * 1. /
      (libMesh::pi * (fuel_number_density + gas_number_density) * std::pow(4.33e-10, 2)) *
      std::sqrt(libMesh::pi * Boltzmann_constant * _temperature[_qp] / (2 *
 equivalent_mass));//[m^2/s]

  RealVectorValue velocity = molecules_volume * diffusion_coeff * molecules_volume_dT * pore_grad;
  Real universal_gas_constant = 8.31446261815324; //[J/mol/K]
  Real Boltzmann_constant = 1.380649e-23;         //[J/K]
  Real Avogadro_constant = 6.02214179e23;         //[molecules/mol]
  Real corros_section_avg = 3e-10;                //[m^2]
  Real uo2_molecules = 5e17;                      //[molecules/m^3]
  Real He_molecules = 2.46e25;                    //[molecules/m^3]
  Real He_molar = 4.0026020e-3;                   //[kg/mol]
  Real uo2_molar = 270.03e-3;                     //[kg/mol]
  Real diffusion_coeff =
      std::sqrt(8 * universal_gas_constant * _temperature[_qp] / libMesh::pi) /
      (3 * libMesh::pi * std::pow(corros_section_avg, 2) * (He_molecules + uo2_molecules)) *
      std::sqrt((He_molar + uo2_molar) / (He_molar * uo2_molar)); //[m^2/s]
  Real pre_factor_pressure = 0;
  Real vaporization = 0;

  if (_temperature[_qp] < 2000)
  {
    pre_factor_pressure = 0.1 * std::exp(33.785);
    vaporization = 5.98e5;
  }
  else if (_temperature[_qp] < 2200)
  {
    pre_factor_pressure = 0.1 * std::exp(34.107);
    vaporization = 6.17e5;
  }
  else
  {
    pre_factor_pressure = 0.1 * std::exp(35.813);
    vaporization = 6.35e5;
  }
  Real vapour_pressure =
      pre_factor_pressure * std::exp(-vaporization / universal_gas_constant / _temperature[_qp]);

  // transport velocity in SI [m/s]
  Real uo2_density =
      1e3 * 10.97 /
      (0.988 + 6.395e-6 * _temperature[_qp] + 3.543e-9 * std::pow(_temperature[_qp], 2) +
       3e-12 * std::pow(_temperature[_qp], 3)); //[kg/m^3]
  // material_flow = -D/kT *dp  velocity = -material_flow * \omega
  RealVectorValue material_flow = -diffusion_coeff / (Boltzmann_constant * _temperature[_qp]) *
                                  vapour_pressure * vaporization /
                                  (universal_gas_constant * std::pow(_temperature[_qp], 2)) * 4 *
                                  _temperature_grad[_qp]; //[1/m^2s])

  Real molecular_volume_uo2 = uo2_molar / (uo2_density * Avogadro_constant); //[m^3]
  RealVectorValue velocity = -material_flow * molecular_volume_uo2;

  RealVectorValue velocity_sens =
      5.006e-12 *
      (0.988 + 6.395e-6 * _temperature[_qp] + 3.543e-9 * std::pow(_temperature[_qp], 2) +
       3e-12 * std::pow(_temperature[_qp], 3)) *
      std::pow(_temperature[_qp], -2.5) * vaporization * vapour_pressure * _temperature_grad[_qp];
  return velocity_sens;
}

Real
PoreAdvection::computeQpResidual()
{
  return -_grad_test[_i][_qp] * (_u[_qp] - _u[_qp] * _u[_qp]) * velocity();
}

Real
PoreAdvection::computeQpJacobian()
{
  return -_grad_test[_i][_qp] * (_phi[_j][_qp] - 2 * _u[_qp] * _phi[_j][_qp]) * velocity();
}

Real
PoreAdvection::negSpeedQp() const
{
  return -_grad_test[_i][_qp] * velocity();
}

void
PoreAdvection::computeResidual()
{
  switch (_upwinding)
  {
    case UpwindType::none:
      Kernel::computeResidual();
      break;
    case UpwindType::full:
      fullUpwind(JacRes::CALCULATE_RESIDUAL);
      break;
  }
}

void
PoreAdvection::computeJacobian()
{
  switch (_upwinding)
  {
    case UpwindType::none:
      Kernel::computeJacobian();
      break;

    case UpwindType::full:
      fullUpwind(JacRes::CALCULATE_JACOBIAN);
      break;
  }
}

void
PoreAdvection::fullUpwind(JacRes res_or_jac)
{
  // The number of nodes in the element
  const unsigned int num_nodes = _test.size();

  // Even if we are computing the Jacobian we still need to compute the outflow from each node to
  // see which nodes are upwind and which are downwind
  prepareVectorTag(_assembly, _var.number());

  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
    prepareMatrixTag(_assembly, _var.number(), _var.number());

  // Compute the outflux from each node and store in _local_re
  // If _local_re is positive at the node, mass (or whatever the Variable represents) is flowing out
  // of the node
  _upwind_node.resize(num_nodes);
  for (_i = 0; _i < num_nodes; ++_i)
  {
    for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      _local_re(_i) += _JxW[_qp] * _coord[_qp] * negSpeedQp();
    _upwind_node[_i] = (_local_re(_i) >= 0.0);
  }

  // Variables used to ensure mass conservation
  Real total_mass_out = 0.0;
  Real total_in = 0.0;
  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
    _dtotal_mass_out.assign(num_nodes, 0.0);

  for (unsigned int n = 0; n < num_nodes; ++n)
  {
    if (_upwind_node[n])
    {
      if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
      {
        if (_test.size() == _phi.size())
          /* u at node=n depends only on the u at node=n, by construction.  For
           * linear-lagrange variables, this means that Jacobian entries involving the derivative
           * will only be nonzero for derivatives wrt variable at node=n.  Hence the
           * (n, n) in the line below.  The above "if" statement catches other variable types
           * (eg constant monomials)
           */
          _local_ke(n, n) += _local_re(n) * (1.0 - 2.0 * _u_nodal[n]);

        _dtotal_mass_out[n] += _local_ke(n, n);
      }
      _local_re(n) *= _u_nodal[n] * (1.0 - _u_nodal[n]);
      total_mass_out += _local_re(n);
    }
    else                        // downwind node
      total_in -= _local_re(n); // note the -= means the result is positive
  }

  // Conserve mass over all phases by proportioning the total_mass_out mass to the inflow nodes,
  // weighted by their local_re values
  for (unsigned int n = 0; n < num_nodes; ++n)
  {
    if (!_upwind_node[n]) // downwind node
    {
      if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
        for (_j = 0; _j < _phi.size(); _j++)
        {
          Real factor = 1.0;
          if (_j < num_nodes && _upwind_node[_j])
            factor = (1.0 - 2.0 * _u_nodal[_j]);
          _local_ke(n, _j) += _local_re(n) * _dtotal_mass_out[_j] * factor / total_in;
        }
      _local_re(n) *= total_mass_out / total_in;
    }
  }

  // Add the result to the residual and jacobian
  if (res_or_jac == JacRes::CALCULATE_RESIDUAL)
  {
    accumulateTaggedLocalResidual();

    if (_has_save_in)
    {
      Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
      for (const auto & var : _save_in)
        var->sys().solution().add_vector(_local_re, var->dofIndices());
    }
  }

  if (res_or_jac == JacRes::CALCULATE_JACOBIAN)
  {
    accumulateTaggedLocalMatrix();

    if (_has_diag_save_in)
    {
      unsigned int rows = _local_ke.m();
      DenseVector<Number> diag(rows);
      for (unsigned int i = 0; i < rows; i++)
        diag(i) = _local_ke(i, i);

      Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
      for (const auto & var : _diag_save_in)
        var->sys().solution().add_vector(diag, var->dofIndices());
    }
  }
}

when I run this, I find the results are more oscillation and instability.

How did I fix this problem? I am not very familiar with upwind finite element.

Thank you very much in advance for your help.
Zhaohao Wang

[GiudGiud]

Hello

I am not that familiar either. The kernel you listed is a good reference.
I would attempt to separate this kernel in 3 parts and only upwind one of them at a time to see which one is causing the oscillations

Guillaume

[wangzhaohao]

Yes, I also thought this method. Three part is good. I also read the Ref in the ConservativeAdvection, Umm, hard to understand the method. So does anyone know of a way to do this? I don't want my code to be too verbose.