diff --git a/framework/include/base/DisplacedSystem.h b/framework/include/base/DisplacedSystem.h index 6ae9e0bff9ce..7e0e49f18819 100644 --- a/framework/include/base/DisplacedSystem.h +++ b/framework/include/base/DisplacedSystem.h @@ -64,6 +64,17 @@ class DisplacedSystem : public SystemTempl */ 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; }; diff --git a/framework/include/base/MooseVariable.h b/framework/include/base/MooseVariable.h index 05496a1184ef..24a76c2fac44 100644 --- a/framework/include/base/MooseVariable.h +++ b/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 diff --git a/framework/include/base/SystemBase.h b/framework/include/base/SystemBase.h index b1504c2c9e52..8fe9ec51712c 100644 --- a/framework/include/base/SystemBase.h +++ b/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 > _var_map; + std::vector _vars_to_be_zeroed_on_residual; + friend class ComputeInitialConditionThread; }; diff --git a/framework/include/bcs/IntegratedBC.h b/framework/include/bcs/IntegratedBC.h index 06d235908dcf..7c07676fafad 100644 --- a/framework/include/bcs/IntegratedBC.h +++ b/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 _local_re; + + /// The aux variables to save the residual contributions to + bool _has_save_in; + std::vector _save_in; + std::vector _save_in_strings; + virtual Real computeQpResidual() = 0; virtual Real computeQpJacobian(); /** diff --git a/framework/include/kernels/Kernel.h b/framework/include/kernels/Kernel.h index f64fb5b41d90..5658b41d6ec6 100644 --- a/framework/include/kernels/Kernel.h +++ b/framework/include/kernels/Kernel.h @@ -126,6 +126,14 @@ class Kernel : MooseArray & _grad_zero; MooseArray & _second_zero; + /// Holds residual entries as their accumulated by this Kernel + DenseVector _local_re; + + /// The aux variables to save the residual contributions to + bool _has_save_in; + std::vector _save_in; + std::vector _save_in_strings; + virtual Real computeQpResidual() = 0; virtual Real computeQpJacobian(); diff --git a/framework/include/postprocessors/NodalSum.h b/framework/include/postprocessors/NodalSum.h new file mode 100644 index 000000000000..32037705ce75 --- /dev/null +++ b/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(); + +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 diff --git a/framework/src/base/FEProblem.C b/framework/src/base/FEProblem.C index b669a03381c8..4b2242db7697 100644 --- a/framework/src/base/FEProblem.C +++ b/framework/src/base/FEProblem.C @@ -1694,6 +1694,10 @@ void FEProblem::computeResidual(NonlinearImplicitSystem & /*sys*/, const NumericVector& soln, NumericVector& residual) { _nl.set_solution(soln); + + _nl.zeroVariables(); + _aux.zeroVariables(); + computePostprocessors(EXEC_RESIDUAL); if (_displaced_problem != NULL) diff --git a/framework/src/base/Moose.C b/framework/src/base/Moose.C index c173fc67c5f6..5fd549b394a8 100644 --- a/framework/src/base/Moose.C +++ b/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); diff --git a/framework/src/base/SystemBase.C b/framework/src/base/SystemBase.C index f0940e586299..fe36599cfe0e 100644 --- a/framework/src/base/SystemBase.C +++ b/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 dof_indices_to_zero = aldit._all_dof_indices; + + NumericVector & solution = this->solution(); + + solution.close(); + + for(std::set::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() { diff --git a/framework/src/bcs/IntegratedBC.C b/framework/src/bcs/IntegratedBC.C index a6b305494392..d65388efd9ca 100644 --- a/framework/src/bcs/IntegratedBC.C +++ b/framework/src/bcs/IntegratedBC.C @@ -21,9 +21,11 @@ template<> InputParameters validParams() { - return validParams(); -} + InputParameters params = validParams(); + params.addParam >("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 >("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 & 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 diff --git a/framework/src/kernels/Kernel.C b/framework/src/kernels/Kernel.C index ae6134d9e7bd..d9ef522b0b12 100644 --- a/framework/src/kernels/Kernel.C +++ b/framework/src/kernels/Kernel.C @@ -19,6 +19,9 @@ #include "SubProblem.h" #include "SystemBase.h" +// libmesh includes +#include "threads.h" + template<> InputParameters validParams() { @@ -26,6 +29,7 @@ InputParameters validParams() params += validParams(); params.addRequiredParam("variable", "The name of the variable that this kernel operates on"); params.addParam >("block", "The list of ids of the blocks (subdomain) that this kernel will be applied to"); + params.addParam >("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("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 >("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 & 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 diff --git a/framework/src/kernels/KernelGrad.C b/framework/src/kernels/KernelGrad.C index b07548249cc9..0d3189499b79 100644 --- a/framework/src/kernels/KernelGrad.C +++ b/framework/src/kernels/KernelGrad.C @@ -15,6 +15,7 @@ #include "KernelGrad.h" #include "SubProblem.h" #include "Moose.h" +#include "SystemBase.h" template<> InputParameters validParams() @@ -39,6 +40,8 @@ KernelGrad::computeResidual() // Moose::perf_log.push("computeResidual()","KernelGrad"); DenseVector & 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; _isys().solution().add_vector(_local_re, _save_in[i]->dofIndices()); + } // Moose::perf_log.pop("computeResidual()","KernelGrad"); } diff --git a/framework/src/kernels/KernelSecond.C b/framework/src/kernels/KernelSecond.C index ea4a97d48435..ee842bddd410 100644 --- a/framework/src/kernels/KernelSecond.C +++ b/framework/src/kernels/KernelSecond.C @@ -14,6 +14,7 @@ #include "KernelSecond.h" #include "SubProblem.h" +#include "SystemBase.h" template<> InputParameters validParams() @@ -37,12 +38,23 @@ void KernelSecond::computeResidual() { DenseVector & re = _assembly.residualBlock(_var.number()); + _local_re.resize(re.size()); + _local_re.zero(); _value.resize(_qrule->n_points()); precalculateResidual(); for (_i=0; _i<_test.size(); _i++) for (_qp=0; _qp<_qrule->n_points(); _qp++) - re(_i) += _JxW[_qp]*_coord[_qp]*(_value[_qp]*_second_test[_i][_qp].tr()); + _local_re(_i) += _JxW[_qp]*_coord[_qp]*(_value[_qp]*_second_test[_i][_qp].tr()); + + 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()); + } } Real diff --git a/framework/src/kernels/KernelValue.C b/framework/src/kernels/KernelValue.C index 3cc52959210e..e672db0cd23f 100644 --- a/framework/src/kernels/KernelValue.C +++ b/framework/src/kernels/KernelValue.C @@ -14,6 +14,7 @@ #include "KernelValue.h" #include "SubProblem.h" +#include "SystemBase.h" template<> InputParameters validParams() @@ -38,12 +39,23 @@ KernelValue::computeResidual() // Moose::perf_log.push("computeResidual()","KernelGrad"); DenseVector & re = _assembly.residualBlock(_var.number()); + _local_re.resize(re.size()); + _local_re.zero(); for (_qp = 0; _qp < _qrule->n_points(); _qp++) { _value = precomputeQpResidual(); for (_i = 0; _i < _test.size(); _i++) - re(_i) += _JxW[_qp] * _coord[_qp] * _value * _test[_i][_qp]; + _local_re(_i) += _JxW[_qp] * _coord[_qp] * _value * _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"); diff --git a/framework/src/postprocessors/NodalSum.C b/framework/src/postprocessors/NodalSum.C new file mode 100644 index 000000000000..7d04638694ef --- /dev/null +++ b/framework/src/postprocessors/NodalSum.C @@ -0,0 +1,59 @@ +/****************************************************************/ +/* 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 */ +/****************************************************************/ + +#include "NodalSum.h" +#include "MooseMesh.h" +#include "SubProblem.h" +// libMesh +#include "boundary_info.h" + +template<> +InputParameters validParams() +{ + InputParameters params = validParams(); + return params; +} + +NodalSum::NodalSum(const std::string & name, InputParameters parameters) : + NodalPostprocessor(name, parameters), + _sum(0) +{ +} + +void +NodalSum::initialize() +{ + _sum = 0; +} + +void +NodalSum::execute() +{ + _sum += _u[_qp]; +} + +Real +NodalSum::getValue() +{ + gatherSum(_sum); + + return _sum; +} + +void +NodalSum::threadJoin(const Postprocessor & y) +{ + const NodalSum & pps = dynamic_cast(y); + _sum += pps._sum; +} diff --git a/test/tests/save_in_test/gold/neumannbc_out.e b/test/tests/save_in_test/gold/neumannbc_out.e new file mode 100644 index 000000000000..4d788702f3b0 Binary files /dev/null and b/test/tests/save_in_test/gold/neumannbc_out.e differ diff --git a/test/tests/save_in_test/gold/out.e b/test/tests/save_in_test/gold/out.e new file mode 100644 index 000000000000..68c55ffb530c Binary files /dev/null and b/test/tests/save_in_test/gold/out.e differ diff --git a/test/tests/save_in_test/save_in_test.i b/test/tests/save_in_test/save_in_test.i new file mode 100644 index 000000000000..fb9004ce00c7 --- /dev/null +++ b/test/tests/save_in_test/save_in_test.i @@ -0,0 +1,74 @@ +[Mesh] + file = square.e + type = GeneratedMesh + dim = 2 + nx = 10 + ny = 10 + nz = 0 + _dimension = 2 +[] + +[Variables] + [./u] + order = FIRST + family = LAGRANGE + [../] +[] + +[AuxVariables] + [./saved] + [../] + [./bc_saved] + [../] + [./accumulated] + [../] +[] + +[Kernels] + [./diff] + type = Diffusion + variable = u + save_in = 'saved accumulated' + [../] +[] + +[BCs] + [./left] + type = DirichletBC + variable = u + boundary = left + value = 0 + [../] + [./nbc] + type = NeumannBC + variable = u + boundary = right + value = 1 + save_in = 'bc_saved accumulated' + [../] +[] + +[Postprocessors] + [./left_flux] + type = NodalSum + variable = saved + boundary = 1 + [../] +[] + +[Executioner] + type = Steady + petsc_options = '-snes_mf_operator' + petsc_options_iname = '-pc_type -pc_hypre_type' + petsc_options_value = 'hypre boomeramg' +[] + +[Output] + file_base = out + output_initial = true + interval = 1 + exodus = true + print_linear_residuals = true + perf_log = true +[] + diff --git a/test/tests/save_in_test/save_in_test.py b/test/tests/save_in_test/save_in_test.py new file mode 100644 index 000000000000..bfdbd25155a6 --- /dev/null +++ b/test/tests/save_in_test/save_in_test.py @@ -0,0 +1,10 @@ +from options import * + +test = { INPUT : 'save_in_test.i', + EXODIFF : ['out.e'], + SCALE_REFINE : 4 + } + + + +