adding the ability for kernels and bcs to save their residual contrib…

…utions to an auxiliary variable closes #1030 #708 #722 #1150

r11203

framework/include/base/DisplacedSystem.h

@@ -64,6 +64,17 @@ class DisplacedSystem : public SystemTempl<TransientExplicitSystem>
    */
   virtual bool currentlyComputingJacobian() { return _undisplaced_system.currentlyComputingJacobian(); }
 
+    /**
+   * Adds this variable to the list of variables to be zeroed during each residual evaluation.
+   * @param var_name The name of the variable to be zeroed.
+   */
+  virtual void addVariableToZeroOnResidual(std::string var_name) { _undisplaced_system.addVariableToZeroOnResidual(var_name); }
+
+  /**
+   * Zero out the variables that have been specified to be zeroed during each residual evaluation.
+   */
+  virtual void zeroVariables() { _undisplaced_system.zeroVariables(); }
+
 protected:
   SystemBase & _undisplaced_system;
 };

framework/include/base/MooseVariable.h

@@ -72,6 +72,11 @@ class MooseVariable
    */
   unsigned int number() { return _moose_var_num; }
 
+  /**
+   * Get the system this variable is part of.
+   */
+  SystemBase & sys() { return _sys; }
+
   /**
    * Get the variable number
    * @return The variable number

framework/include/base/SystemBase.h

@@ -240,6 +240,17 @@ class SystemBase
 
   virtual unsigned int nScalarVariables() = 0;
 
+  /**
+   * Adds this variable to the list of variables to be zeroed during each residual evaluation.
+   * @param var_name The name of the variable to be zeroed.
+   */
+  virtual void addVariableToZeroOnResidual(std::string var_name);
+
+  /**
+   * Zero out the variables that have been specified to be zeroed during each residual evaluation.
+   */
+  virtual void zeroVariables();
+
   /**
    * Get minimal quadrature order needed for integrating variables in this system
    * @return The minimal order of quadrature
@@ -342,6 +353,8 @@ class SystemBase
   /// Map of variables (variable id -> array of subdomains where it lives)
   std::map<unsigned int, std::set<SubdomainID> > _var_map;
 
+  std::vector<std::string> _vars_to_be_zeroed_on_residual;
+
   friend class ComputeInitialConditionThread;
 };
 

framework/include/bcs/IntegratedBC.h

@@ -95,6 +95,14 @@ class IntegratedBC :
   /// the gradient of the unknown variable this BC is acting on
   const VariableGradient & _grad_u;
 
+  /// Holds residual entries as their accumulated by this Kernel
+  DenseVector<Number> _local_re;
+
+  /// The aux variables to save the residual contributions to
+  bool _has_save_in;
+  std::vector<MooseVariable*> _save_in;
+  std::vector<std::string> _save_in_strings;
+
   virtual Real computeQpResidual() = 0;
   virtual Real computeQpJacobian();
   /**

framework/include/kernels/Kernel.h

@@ -126,6 +126,14 @@ class Kernel :
   MooseArray<RealGradient> & _grad_zero;
   MooseArray<RealTensor> & _second_zero;
 
+  /// Holds residual entries as their accumulated by this Kernel
+  DenseVector<Number> _local_re;
+
+  /// The aux variables to save the residual contributions to
+  bool _has_save_in;
+  std::vector<MooseVariable*> _save_in;
+  std::vector<std::string> _save_in_strings;
+
   virtual Real computeQpResidual() = 0;
   virtual Real computeQpJacobian();
 

framework/include/postprocessors/NodalSum.h

@@ -0,0 +1,46 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#ifndef NODALSUM_H
+#define NODALSUM_H
+
+#include "NodalPostprocessor.h"
+
+//Forward Declarations
+class NodalSum;
+class MooseMesh;
+
+template<>
+InputParameters validParams<NodalSum>();
+
+class NodalSum : public NodalPostprocessor
+{
+public:
+  NodalSum(const std::string & name, InputParameters parameters);
+
+  virtual void initialize();
+  virtual void execute();
+
+  /**
+   * This will return the degrees of freedom in the system.
+   */
+  virtual Real getValue();
+
+  void threadJoin(const Postprocessor & y);
+
+protected:
+  Real _sum;
+};
+
+#endif //NODALSUM_H

framework/src/base/FEProblem.C

@@ -1694,6 +1694,10 @@ void
 FEProblem::computeResidual(NonlinearImplicitSystem & /*sys*/, const NumericVector<Number>& soln, NumericVector<Number>& residual)
 {
   _nl.set_solution(soln);
+
+  _nl.zeroVariables();
+  _aux.zeroVariables();
+
   computePostprocessors(EXEC_RESIDUAL);
 
   if (_displaced_problem != NULL)

framework/src/base/Moose.C

@@ -104,6 +104,7 @@
 // PPS
 #include "AverageElementSize.h"
 #include "AverageNodalVariableValue.h"
+#include "NodalSum.h"
 #include "ElementAverageValue.h"
 #include "ElementH1Error.h"
 #include "ElementH1SemiError.h"
@@ -281,6 +282,7 @@ registerObjects()
   // PPS
   registerPostprocessor(AverageElementSize);
   registerPostprocessor(AverageNodalVariableValue);
+  registerPostprocessor(NodalSum);
   registerPostprocessor(ElementAverageValue);
   registerPostprocessor(ElementH1Error);
   registerPostprocessor(ElementH1SemiError);

framework/src/base/SystemBase.C

@@ -18,6 +18,7 @@
 #include "MooseVariable.h"
 #include "Conversion.h"
 #include "Parser.h"
+#include "AllLocalDofIndicesThread.h"
 
 // libMesh
 #include "quadrature_gauss.h"
@@ -96,6 +97,39 @@ SystemBase::getVariableBlocks(unsigned int var_number)
     return & _var_map[var_number];
 }
 
+void
+SystemBase::addVariableToZeroOnResidual(std::string var_name)
+{
+  _vars_to_be_zeroed_on_residual.push_back(var_name);
+}
+
+void
+SystemBase::zeroVariables()
+{
+  if(_vars_to_be_zeroed_on_residual.size() > 0)
+  {
+    AllLocalDofIndicesThread aldit(system(), _vars_to_be_zeroed_on_residual);
+    ConstElemRange & elem_range = *_mesh.getActiveLocalElementRange();
+    Threads::parallel_reduce(elem_range, aldit);
+
+    std::set<unsigned int> dof_indices_to_zero = aldit._all_dof_indices;
+
+    NumericVector<Number> & solution = this->solution();
+
+    solution.close();
+
+    for(std::set<unsigned int>::iterator it = dof_indices_to_zero.begin();
+        it != dof_indices_to_zero.end();
+        ++it)
+      solution.set(*it, 0);
+
+    solution.close();
+
+    // Call update to update the current_local_solution for this system
+    system().update();
+  }
+}
+
 Order
 SystemBase::getMinQuadratureOrder()
 {

framework/src/bcs/IntegratedBC.C

@@ -21,9 +21,11 @@
 template<>
 InputParameters validParams<IntegratedBC>()
 {
-  return validParams<BoundaryCondition>();
-}
+  InputParameters params = validParams<BoundaryCondition>();
+  params.addParam<std::vector<std::string> >("save_in", "The name of auxiliary variables to save this Kernel's residual contributions to.  Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)");
 
+  return params;
+}
 
 IntegratedBC::IntegratedBC(const std::string & name, InputParameters parameters) :
     BoundaryCondition(name, parameters),
@@ -49,20 +51,45 @@ IntegratedBC::IntegratedBC(const std::string & name, InputParameters parameters)
     _grad_test(_var.gradPhiFace()),
 
     _u(_var.sln()),
-    _grad_u(_var.gradSln())
+    _grad_u(_var.gradSln()),
+
+    _save_in_strings(parameters.get<std::vector<std::string> >("save_in"))
 {
+  _save_in.resize(_save_in_strings.size());
+
+  for(unsigned int i=0; i<_save_in_strings.size(); i++)
+  {
+    MooseVariable * var = &_problem.getVariable(_tid, _save_in_strings[i]);
+
+    if(var->feType() != _var.feType())
+      mooseError("Error in " + _name + ". When saving residual values in an Auxiliary variable the AuxVariable must be the same type as the nonlinear variable the object is acting on.");
+
+    _save_in[i] = var;
+    var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
+  }
+
+  _has_save_in = _save_in.size() > 0;
 }
 
 void
 IntegratedBC::computeResidual()
 {
   DenseVector<Number> & re = _assembly.residualBlock(_var.number());
+  _local_re.resize(re.size());
+  _local_re.zero();
 
   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     for (_i = 0; _i < _test.size(); _i++)
-    {
-      re(_i) += _JxW[_qp]*_coord[_qp]*computeQpResidual();
-    }
+      _local_re(_i) += _JxW[_qp]*_coord[_qp]*computeQpResidual();
+
+  re += _local_re;
+
+  if(_has_save_in)
+  {
+    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
+    for(unsigned int i=0; i<_save_in.size(); i++)
+      _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
+  }
 }
 
 void

framework/src/kernels/Kernel.C

@@ -19,13 +19,17 @@
 #include "SubProblem.h"
 #include "SystemBase.h"
 
+// libmesh includes
+#include "threads.h"
+
 template<>
 InputParameters validParams<Kernel>()
 {
   InputParameters params = validParams<MooseObject>();
   params += validParams<TransientInterface>();
   params.addRequiredParam<std::string>("variable", "The name of the variable that this kernel operates on");
   params.addParam<std::vector<SubdomainName> >("block", "The list of ids of the blocks (subdomain) that this kernel will be applied to");
+  params.addParam<std::vector<std::string> >("save_in", "The name of auxiliary variables to save this Kernel's residual contributions to.  Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)");
 
   // testing, dude
   params.addPrivateParam<bool>("use_displaced_mesh", false);
@@ -76,21 +80,46 @@ Kernel::Kernel(const std::string & name, InputParameters parameters) :
     _real_zero(_problem._real_zero[_tid]),
     _zero(_problem._zero[_tid]),
     _grad_zero(_problem._grad_zero[_tid]),
-    _second_zero(_problem._second_zero[_tid])
+    _second_zero(_problem._second_zero[_tid]),
+
+    _save_in_strings(parameters.get<std::vector<std::string> >("save_in"))
 {
+  _save_in.resize(_save_in_strings.size());
+
+  for(unsigned int i=0; i<_save_in_strings.size(); i++)
+  {
+    MooseVariable * var = &_problem.getVariable(_tid, _save_in_strings[i]);
+
+    if(var->feType() != _var.feType())
+      mooseError("Error in " + _name + ". When saving residual values in an Auxiliary variable the AuxVariable must be the same type as the nonlinear variable the object is acting on.");
+
+    _save_in[i] = var;
+    var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
+  }
+
+  _has_save_in = _save_in.size() > 0;
 }
 
 void
 Kernel::computeResidual()
 {
   DenseVector<Number> & re = _assembly.residualBlock(_var.number());
+  _local_re.resize(re.size());
+  _local_re.zero();
 
   precalculateResidual();
   for (_i = 0; _i < _test.size(); _i++)
     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
-    {
-      re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
-    }
+      _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
+
+  re += _local_re;
+
+  if(_has_save_in)
+  {
+    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
+    for(unsigned int i=0; i<_save_in.size(); i++)
+      _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
+  }
 }
 
 void

framework/src/kernels/KernelGrad.C

@@ -15,6 +15,7 @@
 #include "KernelGrad.h"
 #include "SubProblem.h"
 #include "Moose.h"
+#include "SystemBase.h"
 
 template<>
 InputParameters validParams<KernelGrad>()
@@ -39,6 +40,8 @@ KernelGrad::computeResidual()
 //  Moose::perf_log.push("computeResidual()","KernelGrad");
 
   DenseVector<Number> & re = _assembly.residualBlock(_var.number());
+  _local_re.resize(re.size());
+  _local_re.zero();
 
   unsigned int n_qp = _qrule->n_points();
   unsigned int n_test = _test.size();
@@ -51,9 +54,17 @@ KernelGrad::computeResidual()
     Real coord = _coord[_qp];
 
     for (_i=0; _i<n_test; _i++)
-      re(_i) += jxw*coord*_value*_grad_test[_i][_qp];
+      _local_re(_i) += jxw*coord*_value*_grad_test[_i][_qp];
   }
 
+  re += _local_re;
+
+  if(_has_save_in)
+  {
+    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
+    for(unsigned int i=0; i<_save_in.size(); i++)
+      _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
+  }
 //  Moose::perf_log.pop("computeResidual()","KernelGrad");
 }