From a36e89da9e8028ea068c330f516c7cc805d4cdeb Mon Sep 17 00:00:00 2001 From: Derek Gaston Date: Mon, 30 Aug 2010 15:48:42 +0000 Subject: [PATCH] add damper system r2291 --- framework/include/base/ComputeDamping.h | 4 + framework/include/base/MooseSystem.h | 15 +++ framework/include/bcs/BoundaryCondition.h | 1 - framework/include/dampers/ConstantDamper.h | 32 +++++ framework/include/dampers/Damper.h | 84 +++++++++++++ framework/include/dampers/DamperData.h | 89 ++++++++++++++ framework/include/dampers/DamperFactory.h | 83 +++++++++++++ framework/include/dampers/DamperWarehouse.h | 34 ++++++ framework/include/parser/DampersBlock.h | 22 ++++ framework/include/parser/GenericDamperBlock.h | 25 ++++ framework/include/utils/PetscSupport.h | 3 +- framework/src/base/ComputeDamping.C | 113 ++++++++++++++++++ framework/src/base/Moose.C | 10 ++ framework/src/base/MooseSystem.C | 26 ++++ framework/src/dampers/ConstantDamper.C | 26 ++++ framework/src/dampers/Damper.C | 44 +++++++ framework/src/dampers/DamperData.C | 60 ++++++++++ framework/src/dampers/DamperFactory.C | 57 +++++++++ framework/src/dampers/DamperWarehouse.C | 33 +++++ framework/src/parser/DampersBlock.C | 33 +++++ framework/src/parser/GenericDamperBlock.C | 29 +++++ framework/src/utils/PetscSupport.C | 83 ++----------- .../constant_damper_test.i | 72 +++++++++++ .../constant_damper_test.py | 4 + test/tests/constant_damper_test/gold/out.e | Bin 0 -> 2732 bytes test/tests/constant_damper_test/square.e | Bin 0 -> 2812 bytes 26 files changed, 908 insertions(+), 74 deletions(-) create mode 100644 framework/include/base/ComputeDamping.h create mode 100644 framework/include/dampers/ConstantDamper.h create mode 100644 framework/include/dampers/Damper.h create mode 100644 framework/include/dampers/DamperData.h create mode 100644 framework/include/dampers/DamperFactory.h create mode 100644 framework/include/dampers/DamperWarehouse.h create mode 100644 framework/include/parser/DampersBlock.h create mode 100644 framework/include/parser/GenericDamperBlock.h create mode 100644 framework/src/base/ComputeDamping.C create mode 100644 framework/src/dampers/ConstantDamper.C create mode 100644 framework/src/dampers/Damper.C create mode 100644 framework/src/dampers/DamperData.C create mode 100644 framework/src/dampers/DamperFactory.C create mode 100644 framework/src/dampers/DamperWarehouse.C create mode 100644 framework/src/parser/DampersBlock.C create mode 100644 framework/src/parser/GenericDamperBlock.C create mode 100644 test/tests/constant_damper_test/constant_damper_test.i create mode 100644 test/tests/constant_damper_test/constant_damper_test.py create mode 100644 test/tests/constant_damper_test/gold/out.e create mode 100644 test/tests/constant_damper_test/square.e diff --git a/framework/include/base/ComputeDamping.h b/framework/include/base/ComputeDamping.h new file mode 100644 index 000000000000..f91b56b974ad --- /dev/null +++ b/framework/include/base/ComputeDamping.h @@ -0,0 +1,4 @@ +#ifndef COMPUTEPOSTPROCESSORS_H +#define COMPUTEPOSTPROCESSORS_H + +#endif //COMPUTEPOSTPROCESSORS_H diff --git a/framework/include/base/MooseSystem.h b/framework/include/base/MooseSystem.h index 59d3afc10f19..621801822186 100644 --- a/framework/include/base/MooseSystem.h +++ b/framework/include/base/MooseSystem.h @@ -21,6 +21,8 @@ #include "PostprocessorData.h" #include "PostprocessorWarehouse.h" #include "FunctionWarehouse.h" +#include "DamperData.h" +#include "DamperWarehouse.h" //libmesh includes #include "transient_system.h" @@ -224,6 +226,10 @@ class MooseSystem std::string name, InputParameters parameters); + void addDamper(std::string damper_name, + std::string name, + InputParameters parameters); + /** * Computes a block diagonal jacobian for the full system. */ @@ -247,6 +253,8 @@ class MooseSystem void reinitKernels(THREAD_ID tid, const NumericVector& soln, const Elem * elem, DenseVector * Re, DenseMatrix * Ke = NULL); + void reinitDampers(THREAD_ID tid, const NumericVector& increment); + void reinitDGKernels(THREAD_ID tid, const NumericVector& soln, const Elem * elem, const unsigned int side, const Elem * neighbor, DenseVector * Re, bool reinitKe = false); void reinitBCs(THREAD_ID tid, const NumericVector& soln, const Elem * elem, const unsigned int side, const unsigned int boundary_id); @@ -257,6 +265,8 @@ class MooseSystem void compute_postprocessors (const NumericVector& soln); + Real compute_damping(const NumericVector& soln, const NumericVector& update); + void subdomainSetup(THREAD_ID tid, unsigned int block_id); /** @@ -384,6 +394,7 @@ class MooseSystem std::vector _aux_data; std::vector _material_data; std::vector _postprocessor_data; + std::vector _damper_data; DofMap * _dof_map; @@ -427,6 +438,7 @@ class MooseSystem std::vector _ics; std::vector _pps; std::vector _functions; + std::vector _dampers; std::vector _first; @@ -579,6 +591,7 @@ class MooseSystem friend class ComputeInternalResiduals; friend class ComputeInternalPostprocessors; friend class GenericExecutionerBlock; + friend class ComputeInternalDamping; friend class PDEBase; friend class InitialCondition; @@ -592,11 +605,13 @@ class MooseSystem friend class Steady; friend class Postprocessor; friend class FunctionNeumannBC; + friend class Damper; friend class QuadraturePointData; friend class ElementData; friend class FaceData; friend class AuxData; + friend class DamperData; }; /** diff --git a/framework/include/bcs/BoundaryCondition.h b/framework/include/bcs/BoundaryCondition.h index 26a22e8cfd26..631b6441c791 100644 --- a/framework/include/bcs/BoundaryCondition.h +++ b/framework/include/bcs/BoundaryCondition.h @@ -143,7 +143,6 @@ class BoundaryCondition : */ NumericVector * & _current_residual; - /** * Holds the current solution at the current quadrature point on the face. */ diff --git a/framework/include/dampers/ConstantDamper.h b/framework/include/dampers/ConstantDamper.h new file mode 100644 index 000000000000..bcfd5c98736c --- /dev/null +++ b/framework/include/dampers/ConstantDamper.h @@ -0,0 +1,32 @@ +#ifndef CONSTANTDAMPER_H +#define CONSTANTDAMPER_H + +// Moose Includes +#include "Damper.h" + +//Forward Declarations +class ConstantDamper; + +template<> +InputParameters validParams(); + +class ConstantDamper : public Damper +{ +public: + ConstantDamper(std::string name, MooseSystem & moose_system, InputParameters parameters); + +protected: + /** + * This MUST be overriden by a child ConstantDamper. + * + * This is where they actually compute a number between 0 and 1. + */ + virtual Real computeQpDamping(); + + /** + * The constant amount of the newton update to take. + */ + Real _damping; +}; + +#endif //CONSTANTDAMPER_H diff --git a/framework/include/dampers/Damper.h b/framework/include/dampers/Damper.h new file mode 100644 index 000000000000..66f7873f4353 --- /dev/null +++ b/framework/include/dampers/Damper.h @@ -0,0 +1,84 @@ +#ifndef DAMPER_H +#define DAMPER_H + +// Moose Includes +#include "DamperData.h" +#include "PDEBase.h" +#include "MaterialPropertyInterface.h" + +//Forward Declarations +class Damper; + +template<> +InputParameters validParams(); + +class Damper : protected PDEBase, protected MaterialPropertyInterface +{ +public: + Damper(std::string name, MooseSystem & moose_system, InputParameters parameters); + + /** + * Computes this Damper's damping for one element. + */ + Real computeDamping(); + +protected: + /** + * This MUST be overriden by a child damper. + * + * This is where they actually compute a number between 0 and 1. + */ + virtual Real computeQpDamping() = 0; + + /** + * Data associated just with dampers (ie the newton update). + */ + DamperData & _damper_data; + + /** + * Convenience reference to the ElementData object inside of MooseSystem + */ + ElementData & _element_data; + + /** + * The current newton increment. + */ + MooseArray & _u_increment; + + /** + * Holds the current solution at the current quadrature point. + */ + MooseArray & _u; + + /** + * Holds the current solution gradient at the current quadrature point. + */ + MooseArray & _grad_u; + + /** + * Holds the current solution second derivative at the current quadrature point. + */ + MooseArray & _second_u; + + /** + * Holds the previous solution at the current quadrature point. + */ + MooseArray & _u_old; + + /** + * Holds the t-2 solution at the current quadrature point. + */ + MooseArray & _u_older; + + /** + * Holds the previous solution gradient at the current quadrature point. + */ + MooseArray & _grad_u_old; + + /** + * Holds the t-2 solution gradient at the current quadrature point. + */ + MooseArray & _grad_u_older; +}; + +#endif diff --git a/framework/include/dampers/DamperData.h b/framework/include/dampers/DamperData.h new file mode 100644 index 000000000000..7cd891259b78 --- /dev/null +++ b/framework/include/dampers/DamperData.h @@ -0,0 +1,89 @@ +#ifndef DAMPERDATA_H +#define DAMPERDATA_H + +//MOOSE includes +#include "Moose.h" +#include "MooseArray.h" +#include "ElementData.h" + +//Forward Declarations +class MooseSystem; + +namespace libMesh +{ + class QGauss; + class DofMap; + class FEBase; + template class NumericVector; + template class DenseVector; + template class DenseSubVector; + template class DenseMatrix; +} + +/** + * One stop shop for all the data a Kernel class needs. + * + * _One_ of these will get built for each MooseSystem. + */ +class DamperData +{ +public: + DamperData(MooseSystem & moose_system, ElementData & element_data); + + ~DamperData(); + + void init(); + + /** + * Computes the value of the increment for each variable at the quadrature points. + */ + void reinit(const NumericVector& increment_vec); + + /** + * The MooseSystem + */ + MooseSystem & _moose_system; + + MooseArray > _var_increments; + + /** + * A reference to the element data class... so we can reuse that data. + * + * This also means that the element data needs to get reinit BEFORE damper data! + */ + ElementData & _element_data; + + /** + * Dof Maps for all the variables. + */ + std::vector > & _var_dof_indices; + + /** + * quadrature rule. + */ + QGauss * & _qrule; + + /** + * number of quadrature points for current element + */ + unsigned int & _n_qpoints; + + /** + * Map to vector of variable numbers that need to be evaluated + * at the quadrature points + * + * NOTE: This variable is used differently in ElementData and FaceData + * In ElementData, the mapping uses only index zero (0), since all variables lives + * inside the whole domain. + * In FaceData, the mapping goes from boundary id to list of variables active + * on this boundary (not all variables are needed on all boundaries) + */ + std::map > & _var_nums; + + /** + * Shape function. + */ + std::map > *> & _phi; +}; + +#endif //DAMPERDATA_H diff --git a/framework/include/dampers/DamperFactory.h b/framework/include/dampers/DamperFactory.h new file mode 100644 index 000000000000..81a17d004bb9 --- /dev/null +++ b/framework/include/dampers/DamperFactory.h @@ -0,0 +1,83 @@ +#ifndef DAMPERFACTORY_H +#define DAMPERFACTORY_H + +#include "Damper.h" + +// System includes +#include +#include +#include +#include + +// LibMesh includes +#include + +// forward declarations +class MooseSystem; + +/** + * Typedef to make things easier. + */ +typedef Damper * (*damperBuildPtr)(std::string name, MooseSystem & moose_system, InputParameters parameters); + +/** + * Typedef to hide implementation details + */ +typedef std::vector::iterator DamperNamesIterator; + +/** + * Typedef to make things easier. + */ +typedef InputParameters (*damperParamsPtr)(); + +/** + * Templated build function used for generating function pointers to build classes on demand. + */ +template +Damper * buildDamper(std::string name, MooseSystem & moose_system, InputParameters parameters) +{ + return new DamperType(name, moose_system, parameters); +} + +/** + * Responsible for building Dampers on demand and storing them for retrieval + */ +class DamperFactory +{ +public: + static DamperFactory * instance(); + + template + void registerDamper(std::string name) + { + if (_name_to_build_pointer.find(name) == _name_to_build_pointer.end()) + { + _name_to_build_pointer[name] = &buildDamper; + _name_to_params_pointer[name] = &validParams; + } + else + mooseError("Damper '" + name + "' already registered."); + } + + Damper *create(std::string damper_name, std::string name, MooseSystem & moose_system, InputParameters parameters) + { + return (*_name_to_build_pointer[damper_name])(name, moose_system, parameters); + } + + DamperNamesIterator registeredDampersBegin(); + DamperNamesIterator registeredDampersEnd(); + + InputParameters getValidParams(std::string name); + +private: + DamperFactory(); + + virtual ~DamperFactory(); + + std::map _name_to_build_pointer; + std::map _name_to_params_pointer; + + std::vector _registered_damper_names; +}; + +#endif //DAMPERFACTORY_H diff --git a/framework/include/dampers/DamperWarehouse.h b/framework/include/dampers/DamperWarehouse.h new file mode 100644 index 000000000000..576139dfd04b --- /dev/null +++ b/framework/include/dampers/DamperWarehouse.h @@ -0,0 +1,34 @@ +#ifndef DAMPERWAREHOUSE_H +#define DAMPERWAREHOUSE_H + +#include "Damper.h" + +#include +#include +#include + +/** + * Typedef to hide implementation details + */ +typedef std::vector::iterator DamperIterator; + + +/** + * Holds dampers and provides some services + */ +class DamperWarehouse +{ +public: + DamperWarehouse(); + virtual ~DamperWarehouse(); + + DamperIterator dampersBegin(); + DamperIterator dampersEnd(); + + void addDamper(Damper *damper); + +protected: + std::vector _dampers; +}; + +#endif // DAMPERWAREHOUSE_H diff --git a/framework/include/parser/DampersBlock.h b/framework/include/parser/DampersBlock.h new file mode 100644 index 000000000000..59d7c086fd57 --- /dev/null +++ b/framework/include/parser/DampersBlock.h @@ -0,0 +1,22 @@ +#ifndef DAMPERSBLOCK_H +#define DAMPERSBLOCK_H + +#include "ParserBlock.h" + +class DampersBlock; + +template<> +InputParameters validParams(); + +class DampersBlock: public ParserBlock +{ +public: + DampersBlock(std::string name, MooseSystem & moose_system, InputParameters params); + + virtual void execute(); +}; + + + + +#endif //DAMPERSBLOCK_H diff --git a/framework/include/parser/GenericDamperBlock.h b/framework/include/parser/GenericDamperBlock.h new file mode 100644 index 000000000000..16854ef6c628 --- /dev/null +++ b/framework/include/parser/GenericDamperBlock.h @@ -0,0 +1,25 @@ +#ifndef GENERICDAMPERBLOCK_H +#define GENERICDAMPERBLOCK_H + +#include "ParserBlock.h" + +//Forward Declarations +class GenericDamperBlock; + +template<> +InputParameters validParams(); + +class GenericDamperBlock: public ParserBlock +{ +public: + GenericDamperBlock(std::string name, MooseSystem & moose_system, InputParameters params); + + virtual void execute(); + +private: + std::string _type; +}; + + + +#endif //GENERICDAMPERBLOCK_H diff --git a/framework/include/utils/PetscSupport.h b/framework/include/utils/PetscSupport.h index 6de2c11c6623..43103d869038 100644 --- a/framework/include/utils/PetscSupport.h +++ b/framework/include/utils/PetscSupport.h @@ -25,7 +25,8 @@ namespace Moose void petscSetDefaults(MooseSystem &moose_system, Executioner *executioner); - PetscErrorCode petscPhysicsBasedLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w, PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag); +// PetscErrorCode petscPhysicsBasedLineSearch(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w, PetscReal fnorm,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag); + PetscErrorCode dampedCheck(SNES snes, Vec x, Vec y, Vec w, void *lsctx, PetscTruth * changed_y, PetscTruth * changed_w); } } diff --git a/framework/src/base/ComputeDamping.C b/framework/src/base/ComputeDamping.C new file mode 100644 index 000000000000..bcd55692c5c6 --- /dev/null +++ b/framework/src/base/ComputeDamping.C @@ -0,0 +1,113 @@ +//Moose Includes +#include "Moose.h" +#include "MaterialFactory.h" +#include "BoundaryCondition.h" +#include "ParallelUniqueId.h" +#include "MooseSystem.h" +#include "ElementData.h" +#include "DamperWarehouse.h" + +//libMesh includes +#include "numeric_vector.h" +#include "dense_vector.h" +#include "petsc_matrix.h" +#include "dof_map.h" +#include "mesh.h" +#include "boundary_info.h" +#include "elem_range.h" + +#include + +class ComputeInternalDamping +{ +public: + ComputeInternalDamping(MooseSystem &sys, const NumericVector& in_soln, const NumericVector& update) + :_moose_system(sys), + _soln(in_soln), + _update(update), + _damping(1.0) + {} + + // Splitting Constructor + ComputeInternalDamping(ComputeInternalDamping & x, Threads::split) + : _moose_system(x._moose_system), + _soln(x._soln), + _update(x._update), + _damping(1.0) + {} + + // Join Operator + void join(const ComputeInternalDamping & y) + { + if(y._damping < _damping) + _damping = y._damping; + } + + void operator() (const ConstElemRange & range) + { + ParallelUniqueId puid; + + unsigned int tid = puid.id; + + ConstElemRange::const_iterator el = range.begin(); + + DamperIterator damper_begin = _moose_system._dampers[tid].dampersBegin(); + DamperIterator damper_end = _moose_system._dampers[tid].dampersEnd(); + DamperIterator damper_it = damper_begin; + + Real cur_damping; + + for (el = range.begin() ; el != range.end(); ++el) + { + const Elem* elem = *el; + + _moose_system.reinitKernels(tid, _soln, elem, NULL); + _moose_system.reinitDampers(tid, _update); + + unsigned int cur_subdomain = elem->subdomain_id(); + + for(damper_it=damper_begin;damper_it!=damper_end;++damper_it) + { + cur_damping = (*damper_it)->computeDamping(); + if(cur_damping < _damping) + _damping = cur_damping; + } + } + } + + Real _damping; + +protected: + MooseSystem & _moose_system; + const NumericVector & _soln; + const NumericVector & _update; +}; + +Real MooseSystem::compute_damping(const NumericVector& soln, const NumericVector& update) +{ + Moose::perf_log.push("compute_dampers()","Solve"); + + // Default to no damping + Real damping = 1.0; + + // TODO: Make this work with threads! + DamperIterator damper_begin = _dampers[0].dampersBegin(); + DamperIterator damper_end = _dampers[0].dampersEnd(); + DamperIterator damper_it = damper_begin; + + if(damper_begin != damper_end) + { + update_aux_vars(soln); + + ComputeInternalDamping cid(*this, soln, update); + + Threads::parallel_reduce(*getActiveLocalElementRange(), cid); + + damping = cid._damping; + } + + return damping; + + Moose::perf_log.pop("compute_dampers()","Solve"); +} + diff --git a/framework/src/base/Moose.C b/framework/src/base/Moose.C index bbc50381de82..903a1553d195 100644 --- a/framework/src/base/Moose.C +++ b/framework/src/base/Moose.C @@ -104,6 +104,12 @@ #include "PrintNumNodes.h" #include "AverageElementSize.h" +#include "Damper.h" +#include "DamperFactory.h" +#include "DampersBlock.h" +#include "GenericDamperBlock.h" +#include "ConstantDamper.h" + #include "Moose.h" #include "PetscSupport.h" @@ -213,6 +219,8 @@ Moose::registerObjects() ParserBlockFactory::instance()->registerParserBlock("Executioner/Adaptivity"); ParserBlockFactory::instance()->registerParserBlock("Postprocessors"); ParserBlockFactory::instance()->registerParserBlock("Postprocessors/*"); + ParserBlockFactory::instance()->registerParserBlock("Dampers"); + ParserBlockFactory::instance()->registerParserBlock("Dampers/*"); InitialConditionFactory::instance()->registerInitialCondition("ConstantIC"); InitialConditionFactory::instance()->registerInitialCondition("BoundingBoxIC"); @@ -242,6 +250,8 @@ Moose::registerObjects() PostprocessorFactory::instance()->registerPostprocessor("PrintNumElems"); PostprocessorFactory::instance()->registerPostprocessor("PrintNumNodes"); PostprocessorFactory::instance()->registerPostprocessor("AverageElementSize"); + + DamperFactory::instance()->registerDamper("ConstantDamper"); } void diff --git a/framework/src/base/MooseSystem.C b/framework/src/base/MooseSystem.C index e7216ce65e5e..065eef6ebcdc 100644 --- a/framework/src/base/MooseSystem.C +++ b/framework/src/base/MooseSystem.C @@ -14,6 +14,7 @@ #include "ComputeResidual.h" #include "ComputeJacobian.h" #include "ComputeInitialConditions.h" +#include "DamperFactory.h" //libMesh includes #include "numeric_vector.h" @@ -154,6 +155,7 @@ MooseSystem::~MooseSystem() delete _face_data[tid]; delete _neighbor_face_data[tid]; delete _aux_data[tid]; + delete _damper_data[tid]; } if (_preconditioner != NULL) @@ -246,12 +248,15 @@ MooseSystem::sizeEverything() _face_data.resize(n_threads); _neighbor_face_data.resize(n_threads); _aux_data.resize(n_threads); + _damper_data.resize(n_threads); + for (THREAD_ID tid = 0; tid < n_threads; ++tid) { _element_data[tid] = new ElementData(*this, _dof_data[tid]); _face_data[tid] = new FaceData(*this, _dof_data[tid]); _neighbor_face_data[tid] = new FaceData(*this, _neighbor_dof_data[tid]); _aux_data[tid] = new AuxData(*this, _dof_data[tid], *_element_data[tid]); + _damper_data[tid] = new DamperData(*this, *_element_data[tid]); } _kernels.resize(n_threads); @@ -263,6 +268,7 @@ MooseSystem::sizeEverything() _ics.resize(n_threads); _pps.resize(n_threads); _functions.resize(n_threads); + _dampers.resize(n_threads); // Kernels::sizeEverything _bdf2_wei.resize(3); @@ -318,6 +324,7 @@ MooseSystem::init() _face_data[tid]->init(); _neighbor_face_data[tid]->init(); _aux_data[tid]->init(); + _damper_data[tid]->init(); } _t = 0; @@ -925,6 +932,19 @@ MooseSystem::addFunction(std::string func_name, } } +void +MooseSystem::addDamper(std::string damper_name, + std::string name, + InputParameters parameters) +{ + for(THREAD_ID tid=0; tid < libMesh::n_threads(); ++tid) + { + parameters.set("_tid") = tid; + + _dampers[tid].addDamper(DamperFactory::instance()->create(damper_name, name, *this, parameters)); + } +} + void MooseSystem::reinitKernels(THREAD_ID tid, const NumericVector& soln, const Elem * elem, DenseVector * Re, DenseMatrix * Ke) { @@ -1100,6 +1120,12 @@ MooseSystem::reinitAuxKernels(THREAD_ID tid, const NumericVector& soln, _aux_data[tid]->reinit(soln, elem); } +void +MooseSystem::reinitDampers(THREAD_ID tid, const NumericVector& increment) +{ + _damper_data[tid]->reinit(increment); +} + void MooseSystem::updateNewtonStep() { diff --git a/framework/src/dampers/ConstantDamper.C b/framework/src/dampers/ConstantDamper.C new file mode 100644 index 000000000000..c3e753e63b4f --- /dev/null +++ b/framework/src/dampers/ConstantDamper.C @@ -0,0 +1,26 @@ +#include "ConstantDamper.h" + +// Moose includes +#include "MooseSystem.h" + +template<> +InputParameters validParams() +{ + InputParameters params = validParams(); + params.addRequiredParam("damping", "The percentage (between 0 and 1) of the newton update to take."); + return params; +} + +ConstantDamper::ConstantDamper(std::string name, MooseSystem & moose_system, InputParameters parameters) + :Damper(name, moose_system, parameters), + _damping(parameters.get("damping")) +{} + +Real +ConstantDamper::computeQpDamping() +{ + return _damping; +} + + + diff --git a/framework/src/dampers/Damper.C b/framework/src/dampers/Damper.C new file mode 100644 index 000000000000..d2045d9b9065 --- /dev/null +++ b/framework/src/dampers/Damper.C @@ -0,0 +1,44 @@ +#include "Damper.h" + +// Moose includes +#include "MooseSystem.h" + +template<> +InputParameters validParams() +{ + InputParameters params = validParams(); + params.addRequiredParam("variable", "The name of the variable that this damper operates on"); + return params; +} + +Damper::Damper(std::string name, MooseSystem & moose_system, InputParameters parameters) + :PDEBase(name, moose_system, parameters, *moose_system._element_data[parameters.get("_tid")]), + MaterialPropertyInterface(moose_system._material_data[_tid]), + _damper_data(*moose_system._damper_data[_tid]), + _element_data(*moose_system._element_data[_tid]), + _u_increment(_damper_data._var_increments[_var_num]), + _u(_element_data._var_vals[_var_num]), + _grad_u(_element_data._var_grads[_var_num]), + _second_u(_element_data._var_seconds[_var_num]), + _u_old(_element_data._var_vals_old[_var_num]), + _u_older(_element_data._var_vals_older[_var_num]), + _grad_u_old(_element_data._var_grads_old[_var_num]), + _grad_u_older(_element_data._var_grads_older[_var_num]) +{} + +Real +Damper::computeDamping() +{ + Real damping = 1.0; + Real cur_damping = 1.0; + + for (_qp=0; _qp<_qrule->n_points(); _qp++) + { + cur_damping = computeQpDamping(); + if(cur_damping < damping) + damping = cur_damping; + } + + return damping; +} + diff --git a/framework/src/dampers/DamperData.C b/framework/src/dampers/DamperData.C new file mode 100644 index 000000000000..e6698155c469 --- /dev/null +++ b/framework/src/dampers/DamperData.C @@ -0,0 +1,60 @@ +//Moose includes +#include "DamperData.h" +#include "MooseSystem.h" +#include "ComputeQPSolution.h" + +//libmesh includes +#include "numeric_vector.h" +#include "dense_subvector.h" +#include "quadrature_gauss.h" +#include "dof_map.h" +#include "fe_base.h" + +DamperData::DamperData(MooseSystem & moose_system, ElementData & element_data) : + _moose_system(moose_system), + _element_data(element_data), + _var_dof_indices(_element_data._dof_data._var_dof_indices), + _qrule(element_data._qrule), + _n_qpoints(element_data._n_qpoints), + _var_nums(element_data._var_nums), + _phi(element_data._phi) +{} + +DamperData::~DamperData() +{} + +void +DamperData::init() +{ + unsigned int n_vars = _moose_system.getNonlinearSystem()->n_vars(); + _var_increments.resize(n_vars); +} + +void +DamperData::reinit(const NumericVector& increment_vec) +{ + // 0 is for the block id + std::set::iterator var_num_it = _var_nums[0].begin(); + std::set::iterator var_num_end = _var_nums[0].end(); + + for(;var_num_it != var_num_end; ++var_num_it) + { + unsigned int var_num = *var_num_it; + + FEType fe_type = _moose_system._dof_map->variable_type(var_num); + + FEFamily family = fe_type.family; + + _var_increments[var_num].resize(_n_qpoints); + + const std::vector > & static_phi = *_phi[fe_type]; + + std::vector & dof_indices = _var_dof_indices[var_num]; + + MooseArray & increment = _var_increments[var_num]; + + // Compute the increment at each quadrature point + for(unsigned int qp=0; qp<_n_qpoints; qp++) + computeQpSolution(increment[qp], increment_vec, dof_indices, qp, static_phi); + } +} diff --git a/framework/src/dampers/DamperFactory.C b/framework/src/dampers/DamperFactory.C new file mode 100644 index 000000000000..0650a1b377d8 --- /dev/null +++ b/framework/src/dampers/DamperFactory.C @@ -0,0 +1,57 @@ +#include "DamperFactory.h" +#include "MooseSystem.h" + +DamperFactory * +DamperFactory::instance() +{ + static DamperFactory * instance; + if(!instance) + instance=new DamperFactory; + + return instance; +} + +InputParameters +DamperFactory::getValidParams(std::string name) +{ + if( _name_to_params_pointer.find(name) == _name_to_params_pointer.end() ) + { + mooseError(std::string("A _") + name + "_ is not a registered Damper\n\n"); + } + + InputParameters params = _name_to_params_pointer[name](); + return params; +} + +DamperNamesIterator +DamperFactory::registeredDampersBegin() +{ + // Make sure the _registered_damper_names are up to date + _registered_damper_names.clear(); + _registered_damper_names.reserve(_name_to_params_pointer.size()); + + // build a vector of strings from the params pointer map + for (std::map::iterator i = _name_to_params_pointer.begin(); + i != _name_to_params_pointer.end(); + ++i) + { + _registered_damper_names.push_back(i->first); + } + + return _registered_damper_names.begin(); +} + +DamperNamesIterator +DamperFactory::registeredDampersEnd() +{ + return _registered_damper_names.end(); +} + + +DamperFactory::DamperFactory() +{ +} + +DamperFactory:: ~DamperFactory() +{} + diff --git a/framework/src/dampers/DamperWarehouse.C b/framework/src/dampers/DamperWarehouse.C new file mode 100644 index 000000000000..2b12fe16d957 --- /dev/null +++ b/framework/src/dampers/DamperWarehouse.C @@ -0,0 +1,33 @@ +#include "DamperWarehouse.h" + +#include "MooseSystem.h" +#include "Damper.h" + +DamperWarehouse::DamperWarehouse() +{ +} + +DamperWarehouse::~DamperWarehouse() +{ + DamperIterator j; + for (j=_dampers.begin(); j!=_dampers.end(); ++j) + delete *j; +} + +void +DamperWarehouse::addDamper(Damper *damper) +{ + _dampers.push_back(damper); +} + +DamperIterator +DamperWarehouse::dampersBegin() +{ + return _dampers.begin(); +} + +DamperIterator +DamperWarehouse::dampersEnd() +{ + return _dampers.end(); +} diff --git a/framework/src/parser/DampersBlock.C b/framework/src/parser/DampersBlock.C new file mode 100644 index 000000000000..0fb8f6470385 --- /dev/null +++ b/framework/src/parser/DampersBlock.C @@ -0,0 +1,33 @@ +#include "DampersBlock.h" + +#include "DamperFactory.h" + +template<> +InputParameters validParams() +{ + return validParams(); +} + +DampersBlock::DampersBlock(std::string name, MooseSystem & moose_system, InputParameters params) + :ParserBlock(name, moose_system, params) +{ + // Register execution prereqs + addPrereq("Mesh"); + addPrereq("Variables"); + addPrereq("Preconditioning"); + addPrereq("AuxVariables"); + addPrereq("Materials"); +} + +void +DampersBlock::execute() +{ +#ifdef DEBUG + std::cerr << "Inside the DampersBlock Object\n"; +#endif + + // Add the dampers to the system + visitChildren(); +} + + diff --git a/framework/src/parser/GenericDamperBlock.C b/framework/src/parser/GenericDamperBlock.C new file mode 100644 index 000000000000..7ab4a3eb72c2 --- /dev/null +++ b/framework/src/parser/GenericDamperBlock.C @@ -0,0 +1,29 @@ +#include "GenericDamperBlock.h" +#include "DamperFactory.h" +#include "Parser.h" + +template<> +InputParameters validParams() +{ + InputParameters params = validParams(); + return params; +} + +GenericDamperBlock::GenericDamperBlock(std::string name, MooseSystem & moose_system, InputParameters params) + :ParserBlock(name, moose_system, params), + _type(getType()) +{ + setClassParams(DamperFactory::instance()->getValidParams(_type)); +} + +void +GenericDamperBlock::execute() +{ +#ifdef DEBUG + std::cerr << "Inside the GenericDamperBlock Object\n"; + std::cerr << "Damper:" << _type << ":" + << "\tname:" << getShortName() << ":"; +#endif + + _moose_system.addDamper(_type, getShortName(), getClassParams()); +} diff --git a/framework/src/utils/PetscSupport.C b/framework/src/utils/PetscSupport.C index 2c13ad39c14e..0b3ca824fa33 100644 --- a/framework/src/utils/PetscSupport.C +++ b/framework/src/utils/PetscSupport.C @@ -183,7 +183,7 @@ namespace Moose KSPSetPreconditionerSide(ksp, PC_RIGHT); SNESSetMaxLinearSolveFailures(snes, 1000000); -// SNESLineSearchSet(snes, petscPhysicsBasedLineSearch, moose_system.getEquationSystems()); + SNESLineSearchSetPostCheck(snes, dampedCheck, moose_system.getEquationSystems()); #if PETSC_VERSION_LESS_THAN(3,0,0) KSPSetConvergenceTest(ksp, petscConverged, &moose_system); @@ -217,22 +217,21 @@ namespace Moose */ } - PetscErrorCode petscPhysicsBasedLineSearch(SNES snes,void *lsctx,Vec x,Vec /*f*/,Vec g,Vec y,Vec w, PetscReal /*fnorm*/,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag) +// PetscErrorCode petscPhysicsBasedLineSearch(SNES snes,void *lsctx,Vec x,Vec /*f*/,Vec g,Vec y,Vec w, PetscReal /*fnorm*/,PetscReal *ynorm,PetscReal *gnorm,PetscTruth *flag) + PetscErrorCode dampedCheck(SNES snes, Vec x, Vec y, Vec w, void *lsctx, PetscTruth * changed_y, PetscTruth * changed_w) { //w = updated solution = x+ scaling*y //x = current solution //y = updates. // for simple newton use w = x-y - int ierr; - Real facmin = 0.2; - Real max_fraction = 0.5; - - *flag = PETSC_TRUE; + Real damping = 1.0; EquationSystems * equation_systems = static_cast(lsctx); + MooseSystem * moose_system = equation_systems->parameters.get("moose_system"); + MeshBase & mesh = equation_systems->get_mesh(); TransientNonlinearImplicitSystem& system = @@ -245,76 +244,16 @@ namespace Moose PetscVector update_vec_y(y); PetscVector update_vec_w(w); - //vector stores updated solution - update_vec_w.zero(); - update_vec_w.add(1.,update_vec_x); - update_vec_w.add(-1.0,update_vec_y); -/* - MeshBase::const_node_iterator nd = mesh.local_nodes_begin(); - MeshBase::const_node_iterator nd_end = mesh.local_nodes_end(); - MeshBase::const_node_iterator nd_it = mesh.local_nodes_begin();; - - for(nd_it = nd;nd_it != nd_end; ++nd_it) - { - Node * node = *nd_it; - - unsigned int dof = node->dof_number(sys, 0, 0); - -// std::cout<<"x: "<= 1 ) - { - if(update_vec_y(dof)) - { - Real fac = (update_vec_x(dof) - .99995457)/update_vec_y(dof); - if( fac < facmin ) - facmin = fac; - } - } - } - - Parallel::min(facmin); + damping = moose_system->compute_damping(update_vec_w, update_vec_y); - if(facmin < 1.0) + if(damping != 1.0) { - std::cout << std::endl; - std::cout << "facmin..: " << facmin << std::endl; - std::cout << std::endl; + VecScale(y, damping); + *changed_y = PETSC_TRUE; } -*/ - update_vec_w.zero(); - //this is standard newton update with damping parameter - ierr = VecNorm(y,NORM_2,ynorm); - ierr = VecWAXPY(w,-facmin,y,x); -/* - for(nd_it = nd;nd_it != nd_end; ++nd_it) - { - Node * node = *nd_it; - unsigned int dof = node->dof_number(sys, 0, 0); - - if(update_vec_w(dof) <= 0 || update_vec_w(dof) >= 1) - std::cout<<"w: "<DJarGLBUVpz1Fs@LSxmJxC75X`> zD7t|W6DtdyW4({DGN{HQYK+7b!yFEI4Z}-YkB2HNhBTKk656`T^2QMosF+w@<|Hvr z$Jy$jd}&Ebn7iVqn%sC_=4$b#$Z+lyzTsmOnqltPIE=mN}Ht?L-87;9-=sKb2J6qSZ)v8BED zpPypgo`YW-lvQd_F;!QEHtMJvW%#ZW){EcC_%TK@C`#O!+X(i!ziJ9o(gDJIXjNl# z%iC+PuVZvZGjC2%6Y=dAc6H19R!hczHA z_Kk_3&;w_9v5f<_Caj^yNY&#(N0nNP64st?_ww2w}Y3-~i+WqBm(Vlm9Uu@BH z<9Yv?xIVu4|JQ&023U9LHTE#q>^b&6`;Qnk-#GsFTyu|a5;6JYx15KY^B?53aM#IZ;m;L;>bk*Z3lrL7P>R-)(CHtVz&o89f&p)KOTfggbQ zMI1Tu6Z#W+gHl6O4M&og7scn=O=&sD56q1G^(u}A`?bD5l(g(vuW z#A2BU%#Sgu{dAa9WB6JQ1U85aMj?%4I0zy+pgE7>(@xSzC`<66JYpgiB+*aP+QDi3 z&XAaw!7joe7rha-OrDYvVhUaOV~nLcTyqZtv^M3(Olc@W=f)4K;2&$Woc;;4;l4 zobdzio4$D$9){nSG0trlN^6)_YFg@;8z4OUes6dx^|RKALN2YdQ0XXd4TL%us--T5 zVJ-@_A)fE7ciOE^x3$@BZ)4uw+3N0eUKm%EMW{y@O;hK~^Rx)UoE9v_eCzlkM~`_n zZOw2BTOjg?Gq1AH{*2T=K{tD8ng>aU6b*%IJbPJ)9abuJ;+{Aiv57azM9hgj1@{fK zuM|2o)Wf1*m^*vWiSY++48Q|7%hQa7bv)LV^D5JWrtgI#WCC$}HRjywPC!mM2D!wr ztNsz%WA7K8N#A^bPMptg06g34rh{kkzJs?MNu`5+*wdqSH#xI@PFSOL#Qy>9;tl*F1%EH3b`kuE(suF<56e zcrX7#=#T#w-f?yJ*1DXfvFhqUJwJ=*aB%3KTo=RTKZ%pW8nd>JOU)kSzmIdS@t-j2 q-`yKn5T~sLd-5DNZpCW1K91h`x#zFg@AggVNk0$v!wvlJH}DI@7U&iL literal 0 HcmV?d00001