Fixing jacobian bugs in kernels; removing n_gas references from all k…

…ernels
lcpp-org · Jun 5, 2019 · 66a7613 · 66a7613
1 parent 9ba8bad
commit 66a7613
Show file tree

Hide file tree

Showing 10 changed files with 299 additions and 307 deletions.
diff --git a/include/kernels/ProductFirstOrder.h b/include/kernels/ProductFirstOrder.h
@@ -31,11 +31,10 @@ class ProductFirstOrder : public Kernel
 protected:
   Real computeQpResidual() override;
   Real computeQpJacobian() override;
-  Real computeQpOffDiagJacobian(unsigned int) override;
+  Real computeQpOffDiagJacobian(unsigned int jvar) override;
 
   const VariableValue & _v;
   unsigned int _v_id;
-  // const MaterialProperty<Real> & _n_gas;
 
   // The reaction coefficient
   const MaterialProperty<Real> & _reaction_coeff;

diff --git a/include/kernels/ProductSecondOrder.h b/include/kernels/ProductSecondOrder.h
@@ -39,7 +39,6 @@ class ProductSecondOrder : public Kernel
   const VariableValue & _w;
   unsigned int _v_id;
   unsigned int _w_id;
-  const MaterialProperty<Real> & _n_gas;
 
   // The reaction coefficient
   const MaterialProperty<Real> & _reaction_coeff;

diff --git a/include/kernels/ProductThirdOrder.h b/include/kernels/ProductThirdOrder.h
@@ -41,7 +41,6 @@ class ProductThirdOrder : public Kernel
   unsigned int _v_id;
   unsigned int _w_id;
   unsigned int _x_id;
-  const MaterialProperty<Real> & _n_gas;
 
   // The reaction coefficient
   const MaterialProperty<Real> & _reaction_coeff;

diff --git a/include/kernels/ReactantSecondOrder.h b/include/kernels/ReactantSecondOrder.h
@@ -34,11 +34,9 @@ class ReactantSecondOrder : public Kernel
   virtual Real computeQpOffDiagJacobian(unsigned int jvar);
 
   // The reaction coefficient
-  // MooseVariable & _coupled_var_A;
   const MaterialProperty<Real> & _reaction_coeff;
   const VariableValue & _v;
   unsigned int _v_id;
-  const MaterialProperty<Real> & _n_gas;
   Real _stoichiometric_coeff;
   bool _v_eq_u;
 

diff --git a/include/kernels/ReactantThirdOrder.h b/include/kernels/ReactantThirdOrder.h
@@ -40,7 +40,6 @@ class ReactantThirdOrder : public Kernel
   const VariableValue & _w;
   unsigned int _v_id;
   unsigned int _w_id;
-  const MaterialProperty<Real> & _n_gas;
   Real _stoichiometric_coeff;
   bool _v_eq_u;
   bool _w_eq_u;

diff --git a/src/kernels/ProductFirstOrder.C b/src/kernels/ProductFirstOrder.C
@@ -10,7 +10,7 @@ InputParameters
 validParams<ProductFirstOrder>()
 {
   InputParameters params = validParams<Kernel>();
-  params.addCoupledVar("v", "The first variable that is reacting to create u.");
+  params.addRequiredCoupledVar("v", "The first variable that is reacting to create u.");
   params.addRequiredParam<std::string>("reaction", "The full reaction equation.");
   params.addRequiredParam<Real>("coefficient", "The stoichiometric coeffient.");
   params.addParam<bool>("_v_eq_u", false, "Whether or not v and u are the same variable.");
@@ -19,8 +19,8 @@ validParams<ProductFirstOrder>()
 
 ProductFirstOrder::ProductFirstOrder(const InputParameters & parameters)
   : Kernel(parameters),
-    _v(isCoupled("v") ? coupledValue("v") : _zero),
-    _v_id(isCoupled("v") ? coupled("v") : 0),
+    _v(coupledValue("v")),
+    _v_id(coupled("v")),
     _reaction_coeff(getMaterialProperty<Real>("k_" + getParam<std::string>("reaction"))),
     _stoichiometric_coeff(getParam<Real>("coefficient")),
     _v_eq_u(getParam<bool>("_v_eq_u"))
@@ -30,44 +30,23 @@ ProductFirstOrder::ProductFirstOrder(const InputParameters & parameters)
 Real
 ProductFirstOrder::computeQpResidual()
 {
-  if (isCoupled("v"))
-  {
-    return -_test[_i][_qp] * (_stoichiometric_coeff)*_reaction_coeff[_qp] * _v[_qp];
-  }
-  else
-  {
-    return -_test[_i][_qp] * (_stoichiometric_coeff)*_reaction_coeff[_qp] *
-           getMaterialProperty<Real>("n_gas")[_qp];
-  }
+  return -_test[_i][_qp] * (_stoichiometric_coeff)*_reaction_coeff[_qp] * _v[_qp];
 }
 
 Real
 ProductFirstOrder::computeQpJacobian()
 {
-  // return 0.0;
-  if (isCoupled("v"))
-  {
-    if (_v_eq_u)
-      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _phi[_j][_qp];
-    else
-      return 0.0;
-  }
+  if (_v_eq_u)
+    return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _phi[_j][_qp];
   else
     return 0.0;
 }
 
 Real
-ProductFirstOrder::computeQpOffDiagJacobian(unsigned int)
+ProductFirstOrder::computeQpOffDiagJacobian(unsigned int jvar)
 {
-  if (isCoupled("v"))
-  {
-    if (!_v_eq_u)
-      return -_test[_i][_qp] * (_stoichiometric_coeff)*_reaction_coeff[_qp] * _phi[_j][_qp];
-    else
-      return 0.0;
-  }
+  if (!_v_eq_u && jvar == _v_id)
+    return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _phi[_j][_qp];
   else
-  {
-    return 0;
-  }
+    return 0.0;
 }
diff --git a/src/kernels/ProductSecondOrder.C b/src/kernels/ProductSecondOrder.C
@@ -10,8 +10,8 @@ InputParameters
 validParams<ProductSecondOrder>()
 {
   InputParameters params = validParams<Kernel>();
-  params.addCoupledVar("v", "The first variable that is reacting to create u.");
-  params.addCoupledVar("w", "The second variable that is reacting to create u.");
+  params.addRequiredCoupledVar("v", "The first variable that is reacting to create u.");
+  params.addRequiredCoupledVar("w", "The second variable that is reacting to create u.");
   params.addRequiredParam<std::string>("reaction", "The full reaction equation.");
   params.addRequiredParam<Real>("coefficient", "The stoichiometric coeffient.");
   params.addParam<bool>("_v_eq_u", false, "Whether or not v and u are the same variable.");
@@ -21,11 +21,10 @@ validParams<ProductSecondOrder>()
 
 ProductSecondOrder::ProductSecondOrder(const InputParameters & parameters)
   : Kernel(parameters),
-    _v(isCoupled("v") ? coupledValue("v") : _zero),
-    _w(isCoupled("w") ? coupledValue("w") : _zero),
-    _v_id(isCoupled("v") ? coupled("v") : 0),
-    _w_id(isCoupled("w") ? coupled("w") : 0),
-    _n_gas(getMaterialProperty<Real>("n_gas")),
+    _v(coupledValue("v")),
+    _w(coupledValue("w")),
+    _v_id(coupled("v")),
+    _w_id(coupled("w")),
     _reaction_coeff(getMaterialProperty<Real>("k_" + getParam<std::string>("reaction"))),
     _stoichiometric_coeff(getParam<Real>("coefficient")),
     _v_eq_u(getParam<bool>("_v_eq_u")),
@@ -36,98 +35,141 @@ ProductSecondOrder::ProductSecondOrder(const InputParameters & parameters)
 Real
 ProductSecondOrder::computeQpResidual()
 {
-  Real mult1, mult2;
+  return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _v[_qp] * _w[_qp];
+}
 
-  if (isCoupled("v"))
-    mult1 = _v[_qp];
-  else
+Real
+ProductSecondOrder::computeQpJacobian()
+{
+  if (_v_eq_u && _w_eq_u)
   {
-    mult1 = _n_gas[_qp];
+    return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * 2.0 * _v[_qp] *
+           _phi[_j][_qp];
+  }
+  else if (_v_eq_u && !_w_eq_u)
+  {
+    return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _w[_qp] * _phi[_j][_qp];
+  }
+  else if (!_v_eq_u && _w_eq_u)
+  {
+    return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _v[_qp] * _phi[_j][_qp];
   }
-
-  if (isCoupled("w"))
-    mult2 = _w[_qp];
   else
   {
-    mult2 = _n_gas[_qp];
+    return 0.0;
   }
-  return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] *  mult1 * mult2;
-}
-
-Real
-ProductSecondOrder::computeQpJacobian()
-{
-  Real mult1, mult2, power, eq_u_mult, gas_mult;
+  /*
+  Real power, eq_u_mult, gas_mult;
   power = 0.0;
   eq_u_mult = 1.0;
   gas_mult = 1.0;
 
-  if (isCoupled("v"))
-    mult1 = _v[_qp];
-  else
-    mult1 = _n_gas[_qp];
-
-  if (isCoupled("w"))
-    mult2 = _w[_qp];
-  else
-    mult2 = _n_gas[_qp];
-
-  if (isCoupled("v") && _v_eq_u)
+  if (_v_eq_u)
   {
     power += 1.0;
-    eq_u_mult *= mult1;
+    eq_u_mult *= _v[_qp];
   }
   else
-    gas_mult *= mult1;
+    gas_mult *= _v[_qp];
 
-  if (isCoupled("w") && _w_eq_u)
+  if (_w_eq_u)
   {
     power += 1.0;
-    eq_u_mult *= mult2;
+    eq_u_mult *= _w[_qp];
   }
   else
-    gas_mult *= mult2;
+    gas_mult *= _w[_qp];
 
   if (!_v_eq_u && !_w_eq_u)
   {
     return 0.0;
   }
   else
   {
-    return -_test[_i][_qp] * _stoichiometric_coeff * power * _reaction_coeff[_qp] *  std::pow(eq_u_mult, power-1) * gas_mult * _phi[_j][_qp];
+    return -_test[_i][_qp] * _stoichiometric_coeff * power * _reaction_coeff[_qp] *
+           std::pow(eq_u_mult, power - 1) * gas_mult * _phi[_j][_qp];
   }
+  */
 }
 
 Real
 ProductSecondOrder::computeQpOffDiagJacobian(unsigned int jvar)
 {
-  Real mult1, mult2;
-
-  if (isCoupled("v"))
-    mult1 = _v[_qp];
+  if (_v_eq_u && !_w_eq_u)
+  {
+    if (jvar == _w_id)
+      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _v[_qp] *
+             _phi[_j][_qp];
+    else
+      return 0.0;
+  }
+  else if (!_v_eq_u && _w_eq_u)
+  {
+    if (jvar == _w_id)
+      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _w[_qp] *
+             _phi[_j][_qp];
+    else
+      return 0.0;
+  }
+  else if (!_v_eq_u && !_w_eq_u)
+  {
+    if (jvar == _v_id && jvar == _w_id)
+      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * 2.0 * _v[_qp] *
+             _phi[_j][_qp];
+    else if (jvar == _v_id && jvar != _w_id)
+      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _w[_qp] *
+             _phi[_j][_qp];
+    else if (jvar != _v_id && jvar == _w_id)
+      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * _v[_qp] *
+             _phi[_j][_qp];
+    else
+      return 0.0;
+  }
   else
-    mult1 = _n_gas[_qp];
+    return 0.0;
+  /*
+  Real eq_u_mult, j_mult, power;
+  eq_u_mult = 1.0;
+  j_mult = 1.0;
+  power = 0.0;
 
-  if (isCoupled("w"))
-    mult2 = _w[_qp];
-  else
-    mult2 = _n_gas[_qp];
+  if (_v_eq_u)
+  {
+    eq_u_mult *= _v[_qp];
+  }
+  else if (jvar == _v_id)
+  {
+    power += 1.0;
+    j_mult *= _v[_qp];
+  }
+
+  if (_w_eq_u)
+  {
+    eq_u_mult *= _w[_qp];
+  }
+  else if (jvar == _w_id)
+  {
+    power += 1.0;
+    j_mult *= _w[_qp];
+  }
 
-  if ((_v_eq_u && !isCoupled("w")) || (_w_eq_u && !isCoupled("v")) || (!isCoupled("v") && !isCoupled("w")) || (_v_eq_u && _w_eq_u))
+  return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * std::pow(j_mult, power) *
+         eq_u_mult * _phi[_j][_qp];
+
+  */
+
+  /*if (_v_eq_u && _w_eq_u)
   {
     return 0.0;
   }
-  else if ((isCoupled("v") && !_v_eq_u) && !isCoupled("w"))
+  else if (!_v_eq_u && jvar == _v_id)
   {
-    if (jvar == _v_id)
-      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] *  mult2 * _phi[_j][_qp];
-    else
-      return 0.0;
+    return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * *_phi[_j][_qp];
   }
   else if ((isCoupled("w") && !_w_eq_u) && !isCoupled("v"))
   {
     if (jvar == _w_id)
-      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] *  mult1 * _phi[_j][_qp];
+      return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * mult1 * _phi[_j][_qp];
     else
       return 0.0;
   }
@@ -138,42 +180,47 @@ ProductSecondOrder::computeQpOffDiagJacobian(unsigned int jvar)
       if (_v_id != _w_id)
       {
         if (jvar == _v_id)
-          return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] *  mult2 * _phi[_j][_qp];
+          return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * mult2 *
+                 _phi[_j][_qp];
         else if (jvar == _w_id)
-          return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] *  mult1 * _phi[_j][_qp];
+          return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * mult1 *
+                 _phi[_j][_qp];
         else
           return 0.0;
       }
       else
       {
         if (jvar == _v_id)
-          return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * 2.0 * mult1 * _phi[_j][_qp];
+          return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * 2.0 * mult1 *
+                 _phi[_j][_qp];
         else
           return 0.0;
       }
     }
     else if (_v_eq_u && !_w_eq_u)
     {
       if (jvar == _w_id)
-        return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] *  mult1 * _phi[_j][_qp];
+        return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * mult1 *
+               _phi[_j][_qp];
       else
         return 0.0;
     }
     else if (!_v_eq_u && _w_eq_u)
     {
       if (jvar == _v_id)
-        return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] *  mult2 * _phi[_j][_qp];
+        return -_test[_i][_qp] * _stoichiometric_coeff * _reaction_coeff[_qp] * mult2 *
+               _phi[_j][_qp];
       else
         return 0.0;
     }
     else
       return 0.0;
-
   }
   else
   {
-    // std::cout << getParam<std::string>("reaction") << ": " << _v_eq_u << ", " << _w_eq_u << std::endl;
-    mooseError("ProductSecondOrder, computeQpOffDiagJacobian: this is not yet implemented for the current case.");
-  }
-
+    // std::cout << getParam<std::string>("reaction") << ": " << _v_eq_u << ", " << _w_eq_u <<
+    // std::endl;
+    mooseError("ProductSecondOrder, computeQpOffDiagJacobian: this is not yet implemented for the "
+               "current case.");
+  }*/
 }