Merge pull request idaholab#4203 from andrsd/ode-fix

Fixing computation of off-diagonal entries in ODEKernel
andrsd · Nov 11, 2014 · 3d3adc7 · 3d3adc7
2 parents 31064b1 + 35c05ff
commit 3d3adc7
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 16 deletions.
diff --git a/framework/src/kernels/ODEKernel.C b/framework/src/kernels/ODEKernel.C
@@ -51,32 +51,28 @@ ODEKernel::computeJacobian()
 
   for (_i = 0; _i < _var.order(); _i++)
     for (_j = 0; _j < _var.order(); _j++)
-    {
-      if (_i == _j)
-        ke(_i, _j) += computeQpJacobian();
-      else
-        ke(_i, _j) += computeQpOffDiagJacobian(_var.number());
-    }
+      ke(_i, _j) += computeQpJacobian();
+
+  // compute off-diagonal jacobians wrt scalar variables
+  const std::vector<MooseVariableScalar *> & scalar_vars = _sys.getScalarVariables(_tid);
+  for (std::vector<MooseVariableScalar *>::const_iterator it = scalar_vars.begin(); it != scalar_vars.end(); ++it)
+  {
+    MooseVariableScalar * jvar = *it;
+    computeOffDiagJacobian(jvar->number());
+  }
 }
 
 void
 ODEKernel::computeOffDiagJacobian(unsigned int jvar)
 {
-  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
   if (_sys.isScalarVariable(jvar))
   {
+    DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
     MooseVariableScalar & var_j = _sys.getScalarVariable(_tid, jvar);
     for (_i = 0; _i < _var.order(); _i++)
       for (_j = 0; _j < var_j.order(); _j++)
       {
-        if (jvar == _var.number())
-        {
-          if (_i == _j)
-            ke(_i, _j) += computeQpJacobian();
-          else
-            ke(_i, _j) += computeQpOffDiagJacobian(_var.number());
-        }
-        else
+        if (jvar != _var.number())
           ke(_i, _j) += computeQpOffDiagJacobian(jvar);
       }
   }

diff --git a/framework/src/kernels/ODETimeDerivative.C b/framework/src/kernels/ODETimeDerivative.C
@@ -35,5 +35,8 @@ ODETimeDerivative::computeQpResidual()
 Real
 ODETimeDerivative::computeQpJacobian()
 {
-  return _du_dot_du[_i];
+  if (_i == _j)
+    return _du_dot_du[_i];
+  else
+    return 0;
 }