Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Andreas Vogel, Christian Wehner
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__LIB_ALGEBRA__OPERATOR__DAMPING__
34 : #define __H__LIB_ALGEBRA__OPERATOR__DAMPING__
35 :
36 : #include "common/common.h"
37 : #include "common/util/smart_pointer.h"
38 : #include "interface/linear_operator.h"
39 :
40 : namespace ug{
41 :
42 : /**
43 : * Base class for damping of correction in iterative schemes. An iteration for the
44 : * solution of a matrix problem \f$ A*x = b \f$ is given, for example, by
45 : *
46 : * \f[ x = x + c, \f]
47 : *
48 : * where \f$ c = B*d \f$ is some proposed correction.
49 : *
50 : * The damping class now computes a damping factor \f$ \kappa \f$, that is used to
51 : * damp the correction, i.e.,
52 : *
53 : * \f[ x = x + \kappa c. \f]
54 : *
55 : * In general, the damping may depend on the correction, the (old) defect and
56 : * the operator A itself.
57 : */
58 : template <typename X, typename Y = X>
59 : class IDamping
60 : {
61 : public:
62 : /// returns the damping
63 : /**
64 : * For a given correction, defect and Operator the damping is returned.
65 : *
66 : * @param c the correction
67 : * @param d the defect
68 : * @param spLinOp the operator
69 : * @return the damping
70 : */
71 : virtual number damping(const Y& c, const X& d, ConstSmartPtr<ILinearOperator<Y,X> > spLinOp) const = 0;
72 :
73 : /// returns if the damping is constant
74 : /**
75 : * returns if the damping is constant, i.e. does not depend on correction,
76 : * defect and linear operator.
77 : *
78 : * @return true if constant damping
79 : */
80 : virtual bool constant_damping() const = 0;
81 :
82 : /// returns the constant damping, throws exception if non-constant damping
83 : virtual number damping() const = 0;
84 :
85 : /// virtual destructor
86 : virtual ~IDamping() {}
87 :
88 : /// returns information about configuration parameters
89 : /**
90 : * this should return necessary information about parameters and possibly
91 : * calling config_string of subcomponents.
92 : *
93 : * \returns std::string necessary information about configuration parameters
94 : */
95 : virtual std::string config_string() const = 0;
96 : };
97 :
98 : /// constant damping factor
99 : template <typename X, typename Y = X>
100 : class ConstantDamping : public IDamping<X,Y>
101 : {
102 : public:
103 0 : ConstantDamping(number factor) : m_factor(factor) {}
104 :
105 : /// returns the constant damping factor
106 0 : virtual number damping(const Y& c, const X& d, ConstSmartPtr<ILinearOperator<Y,X> > spLinOp) const
107 : {
108 0 : return m_factor;
109 : }
110 :
111 : /// returns the constant damping factor
112 0 : virtual number damping() const
113 : {
114 0 : return m_factor;
115 : }
116 :
117 : /// returns if damping is constant
118 0 : virtual bool constant_damping() const {return true;};
119 :
120 0 : virtual std::string config_string() const
121 : {
122 0 : std::stringstream ss; ss << "ConstantDamping(" << m_factor << ")"; return ss.str();
123 0 : }
124 :
125 : protected:
126 : number m_factor; ///< constant damping factor
127 : };
128 :
129 : /// damping computed based on the minimal residuum
130 : template <typename X, typename Y = X>
131 0 : class MinimalResiduumDamping : public IDamping<X,Y>
132 : {
133 : public:
134 : /// returns the damping factor
135 0 : virtual number damping(const Y& c, const X& d, ConstSmartPtr<ILinearOperator<Y,X> > spLinOp) const
136 : {
137 0 : SmartPtr<X> spAc = d.clone_without_values();
138 : X& Ac = *spAc;
139 :
140 : try{
141 0 : spLinOp.cast_const()->apply(Ac, c);
142 0 : }UG_CATCH_THROW("MinimalResiduumDamping: Computing Ac failed.")
143 :
144 : // Compute scaling
145 : try{
146 0 : const number kappa = VecProd(Ac, d) / VecProd(Ac, Ac);
147 :
148 0 : if (kappa<0.3) return 0.3;
149 :
150 : // return result
151 : return kappa;
152 : }UG_CATCH_THROW("MinimalResiduumDamping: Computing (d,Ac)/(Ac,Ac) failed.")
153 : }
154 :
155 : /// returns if damping is constant
156 0 : virtual bool constant_damping() const {return false;};
157 :
158 : /// returns the constant damping factor
159 0 : virtual number damping() const
160 : {
161 0 : UG_THROW("MinimalResiduumDamping: non-constant damping.");
162 : }
163 :
164 0 : virtual std::string config_string() const
165 : {
166 0 : return "MinimalResiduumDamping";
167 : }
168 : };
169 :
170 : /// damping computed based on the minimal energy
171 : template <typename X, typename Y = X>
172 0 : class MinimalEnergyDamping : public IDamping<X,Y>
173 : {
174 : public:
175 : /// returns the damping factor
176 0 : virtual number damping(const Y& c, const X& d, ConstSmartPtr<ILinearOperator<Y,X> > spLinOp) const
177 : {
178 0 : SmartPtr<X> spAc = d.clone_without_values();
179 : X& Ac = *spAc;
180 0 : spLinOp.cast_const()->apply(Ac, c);
181 :
182 : // Compute scaling
183 0 : const number kappa = VecProd(d,c) / VecProd(Ac, c);
184 :
185 0 : if (kappa<0.3) return 0.3;
186 :
187 : // return result
188 : return kappa;
189 : }
190 :
191 : /// returns if damping is constant
192 0 : virtual bool constant_damping() const {return false;};
193 :
194 : /// returns the constant damping factor
195 0 : virtual number damping() const
196 : {
197 0 : UG_THROW("MinimalEnergyDamping: non-constant damping.");
198 : }
199 :
200 0 : virtual std::string config_string() const
201 : {
202 0 : return "MinimalEnergyDamping";
203 : }
204 : };
205 :
206 : } // end namespace ug
207 :
208 : #endif /* __H__LIB_ALGEBRA__OPERATOR__DAMPING__ */
|
