Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel, Arne Nägel
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__UG__LIB_ALGEBRA__OPERATOR__LINEAR_OPERATOR__TRANSFORMING_ITER__
34 : #define __H__UG__LIB_ALGEBRA__OPERATOR__LINEAR_OPERATOR__TRANSFORMING_ITER__
35 :
36 : #include "lib_algebra/operator/interface/linear_iterator.h"
37 : #include "lib_algebra/operator/interface/matrix_operator.h"
38 : #include "lib_algebra/operator/debug_writer.h"
39 :
40 : #include "lib_disc/assemble_interface.h"
41 : #include "lib_disc/function_spaces/grid_function.h"
42 :
43 : namespace ug{
44 :
45 : /**
46 : * Abstract base class for transforming iterations.
47 : * Supporting both left- and right transformations
48 : *
49 : * Given
50 : * \f{eqnarray*}{
51 : * A = T_L^{-1} \hat{A} T_R
52 : * \f}
53 : * this implements a subspace correction based on a defect correction:
54 : * \f{eqnarray*}{
55 : * x := x + T_R^{-1} {\hat{A}}^{-1} T_L (b-Ax)
56 : * \f}
57 : * If inversion is to expensive, we replace may replace this by a (single step) iterative solver.
58 : *
59 : * In order
60 : *
61 : *
62 : */
63 : template<typename TAlgebra, typename TDerived>
64 : class ITransformingIteration :
65 : public ILinearIterator<typename TAlgebra::vector_type>,
66 : public DebugWritingObject<TAlgebra>
67 : {
68 : private:
69 :
70 : /// CRTP operator
71 : TDerived& derived()
72 : { return *static_cast<TDerived*>(this);}
73 :
74 : public:
75 :
76 : /// Algebra type
77 : typedef TAlgebra algebra_type;
78 :
79 : /// Vector type
80 : typedef typename algebra_type::vector_type vector_type;
81 :
82 : /// Matrix type
83 : typedef typename algebra_type::matrix_type matrix_type;
84 :
85 :
86 0 : ITransformingIteration() : DebugWritingObject<TAlgebra>() {}
87 : ITransformingIteration(const ITransformingIteration &ti)
88 : { UG_THROW("Do you really want to call a copy constructor of this abstract base class?"); }
89 :
90 :
91 : // Begin CRTP
92 :
93 : /// implementation of init for non-linear (CRTP)
94 0 : bool init(SmartPtr<ILinearOperator<vector_type> > J, const vector_type& u)
95 0 : { return derived().init(J,u);};
96 :
97 : /// implementation of init for linear (CRTP)
98 0 : bool init(SmartPtr<ILinearOperator<vector_type> > L)
99 0 : { return derived().init(L);};
100 :
101 : /// original operator (CRTP)
102 : SmartPtr<ILinearOperator<vector_type> > original_operator()
103 : { return derived().original_operator();}
104 :
105 : /// map: d -> dtilde (CRTP)
106 : bool transform_defect(vector_type &c, const vector_type& d)
107 : { return derived().transform_defect(c,d); }
108 :
109 : /// map: dtilde -> ctilde (CRTP)
110 : bool apply_transformed(vector_type &c, const vector_type &d)
111 : { return derived().apply_transformed(c,d); }
112 :
113 : /// map: ctilde -> c (CRTP)
114 : bool untransform_correction(vector_type& c, const vector_type& d)
115 : { return derived().untransform_correction(c,d); }
116 :
117 :
118 : // End CRTP
119 :
120 : /// implementation of apply (final, non-CRTP)
121 0 : virtual bool apply(vector_type& c, const vector_type& d)
122 : {
123 0 : return (derived().transform_defect(c,d) &&
124 0 : derived().apply_transformed(c,d) &&
125 0 : derived().untransform_correction(c,d));
126 : }
127 :
128 :
129 : /// implementation of apply (final, non-CRTP)
130 0 : virtual bool apply_update_defect(vector_type &c, vector_type& d)
131 : {
132 : // compute correction
133 0 : if (!apply(c,d)) return false;
134 :
135 : // compute updated defect
136 0 : original_operator()->apply_sub(d, c);
137 :
138 0 : return true;
139 : }
140 :
141 :
142 :
143 : };
144 :
145 :
146 : /***
147 : * Implementation of Transforming-/DGS smoothers, following Wittum, 1989.
148 : *
149 : *
150 : * A) For DGS type, assemble
151 : *
152 : * \hat A = \begin{pmatrix} Q& W\\
153 : * -div_h & -\triangle_h^N \end{pmatrix}
154 : *
155 : * T_R**{-1} = \begin{pmatrix} I & grad_h\\
156 : * 0 & -Q_h**l \end{pmatrix}
157 : *
158 : * and specify a preconditioner for $\hat A$.
159 : *
160 : *
161 : * B) For Wittum-type transforming smoothers, assemble
162 : *
163 : * \hat A = \begin{pmatrix} Q& W\\
164 : * -div_h & -\triangle_h^N \end{pmatrix}
165 : *
166 : * T_R**{-1} = \begin{pmatrix} A & B \\
167 : * 0 & I \end{pmatrix}
168 : *
169 : *
170 : * and specify a preconditioner for both $\hat A$ and T_R.
171 : *
172 : * */
173 : template <typename TDomain, typename TAlgebra>
174 : class AssembledTransformingSmoother :
175 : //public ILinearIterator<typename TAlgebra::vector_type>,
176 : //public DebugWritingObject<TAlgebra>
177 : public ITransformingIteration<TAlgebra, AssembledTransformingSmoother<TDomain,TAlgebra> >
178 : {
179 : public:
180 : /// Domain
181 : typedef TDomain domain_type;
182 :
183 : /// Algebra type
184 : typedef TAlgebra algebra_type;
185 :
186 : /// Vector type
187 : typedef typename algebra_type::vector_type vector_type;
188 :
189 : /// Matrix type
190 : typedef typename algebra_type::matrix_type matrix_type;
191 :
192 : using DebugWritingObject<TAlgebra>::write_debug;
193 : using ILinearIterator<typename TAlgebra::vector_type>::damping;
194 :
195 : public:
196 : /// constructor setting approximation space
197 0 : AssembledTransformingSmoother(SmartPtr<IAssemble<TAlgebra> > spAuxSystemAss,
198 : SmartPtr<ILinearIterator<vector_type> > spAuxSmoother,
199 : SmartPtr<IAssemble<TAlgebra> > spRightTrafoAss,
200 : SmartPtr<ILinearIterator<vector_type> > spRightTrafoSmoother=SPNULL)
201 : : m_spAuxSystemAss(spAuxSystemAss),
202 0 : m_spAuxMatrixOperator(new MatrixOperator<matrix_type, vector_type>),
203 : m_spAuxSmoother(spAuxSmoother),
204 : m_spRightTrafoAss(spRightTrafoAss),
205 0 : m_spRightTrafoMat(new MatrixOperator<matrix_type, vector_type>),
206 0 : m_spRightTrafoSmoother(spRightTrafoSmoother)
207 0 : {}
208 :
209 : /// name
210 0 : virtual const char* name() const {return "Assembled Transform Smoother";}
211 :
212 : /// returns if parallel solving is supported
213 0 : virtual bool supports_parallel() const
214 0 : {return (m_spAuxSmoother.valid()) ? m_spAuxSmoother->supports_parallel() : false;}
215 :
216 : // Begin CRTP functions
217 : bool init(SmartPtr<ILinearOperator<vector_type> > J, const vector_type& u);
218 :
219 : /// Since we need grid information, linear operators are not supported...
220 0 : bool init(SmartPtr<ILinearOperator<vector_type> > L)
221 0 : { UG_THROW("Not Implemented!!"); }
222 :
223 : bool transform_defect(vector_type& c, const vector_type& d) {return true;}
224 : bool apply_transformed(vector_type &c, const vector_type &d);
225 : bool untransform_correction(vector_type& c, const vector_type& d);
226 :
227 :
228 : /// Clone
229 : SmartPtr<ILinearIterator<vector_type> > clone();
230 :
231 :
232 :
233 :
234 : SmartPtr<ILinearOperator<vector_type> > original_operator()
235 : {return m_spOriginalSystemOp;}
236 :
237 : protected:
238 :
239 : /// matrix for original system
240 : SmartPtr<ILinearOperator<vector_type> > m_spOriginalSystemOp;
241 :
242 : /// auxiliary correction vector
243 : SmartPtr<vector_type> m_spAuxCorrection;
244 :
245 : /// assembling the transformed system
246 : SmartPtr<IAssemble<TAlgebra> > m_spAuxSystemAss;
247 :
248 : /// matrix (operator) of transformed system
249 : SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spAuxMatrixOperator;
250 :
251 : /// smoother used on transformed system
252 : SmartPtr<ILinearIterator<vector_type> > m_spAuxSmoother;
253 :
254 :
255 : /// assembling the right-transformation
256 : SmartPtr<IAssemble<TAlgebra> > m_spRightTrafoAss;
257 :
258 : /// matrix (operator) of right-transformation
259 : SmartPtr<MatrixOperator<matrix_type, vector_type> > m_spRightTrafoMat;
260 :
261 : /// smoother used on transformed system
262 : SmartPtr<ILinearIterator<vector_type> > m_spRightTrafoSmoother;
263 :
264 :
265 :
266 :
267 : };
268 :
269 :
270 :
271 : template <typename TDomain, typename TAlgebra>
272 0 : bool AssembledTransformingSmoother<TDomain, TAlgebra>::
273 : init(SmartPtr<ILinearOperator<vector_type> > J, const vector_type& u)
274 : {
275 : //m_spOriginalSystemOp = J;
276 :
277 : GridLevel gridLevel;
278 :
279 : const GridFunction<TDomain, TAlgebra>* pSol =
280 0 : dynamic_cast<const GridFunction<TDomain, TAlgebra>*>(&u);
281 0 : if(pSol){
282 0 : gridLevel = pSol->dof_distribution()->grid_level();
283 : }
284 :
285 : // init aux operator
286 : try{
287 0 : m_spAuxSystemAss->assemble_jacobian(*m_spAuxMatrixOperator, u, gridLevel);
288 : }
289 0 : UG_CATCH_THROW("AssembledTransformingSmoother: "
290 : " Cannot assemble transformed system matrix.");
291 :
292 0 : write_debug(*m_spAuxMatrixOperator, "TrafoSystem");
293 :
294 0 : if(!m_spAuxSmoother->init(m_spAuxMatrixOperator, u))
295 0 : UG_THROW("AssembledTransformingSmoother: "
296 : " Cannot init smoother for transformed system matrix.");
297 :
298 :
299 : // init right transform
300 : try{
301 0 : m_spRightTrafoAss->assemble_jacobian(*m_spRightTrafoMat, u, gridLevel);
302 : }
303 0 : UG_CATCH_THROW("AssembledTransformingSmoother: "
304 : " Cannot assemble right-transformation matrix.");
305 :
306 0 : write_debug(*m_spRightTrafoMat, "RightTrafo");
307 :
308 :
309 0 : return true;
310 : };
311 :
312 :
313 :
314 :
315 : /**
316 : *
317 : * */
318 : template <typename TDomain, typename TAlgebra>
319 0 : bool AssembledTransformingSmoother<TDomain, TAlgebra>::
320 : apply_transformed(vector_type &c, const vector_type& d)
321 : {
322 : /**
323 : // temporary vector for defect
324 : vector_type dTmp; dTmp.resize(d.size());
325 :
326 : // copy defect
327 : dTmp = d;
328 :
329 : // work on copy
330 : return apply_update_defect(c, dTmp);
331 :
332 : */
333 :
334 : // obtain aux correction vector
335 0 : m_spAuxCorrection = c.clone_without_values();
336 :
337 : // apply smoother of transformed system
338 0 : if(!m_spAuxSmoother->apply(*m_spAuxCorrection, d))
339 : {
340 : UG_LOG("AssembledTransformingSmoother: Smoother applied incorrectly.\n");
341 0 : return false;
342 : }
343 :
344 : return true;
345 : }
346 :
347 : template <typename TDomain, typename TAlgebra>
348 0 : bool AssembledTransformingSmoother<TDomain, TAlgebra>::
349 : untransform_correction(vector_type& c, const vector_type& d)
350 : {
351 :
352 : // apply right-transformation
353 0 : if (m_spRightTrafoSmoother.valid())
354 : {
355 : // e.g., for Wittum-type smoothers
356 0 : if(!m_spRightTrafoSmoother->apply(*m_spAuxCorrection, d))
357 : {
358 : UG_LOG("AssembledTransformingSmoother: Right Trafo smoother failed!\n");
359 0 : return false;
360 : }
361 : }
362 : else
363 : {
364 : // e.g., for Brandt-Dinar-DGS-type smoothers
365 0 : m_spRightTrafoMat->apply(c, *m_spAuxCorrection);
366 : }
367 :
368 :
369 : // release auy correction vecto
370 0 : m_spAuxCorrection = SPNULL;
371 :
372 : #ifdef UG_PARALLEL
373 : c.change_storage_type(PST_CONSISTENT);
374 : #endif
375 :
376 0 : const number damp = this->damping()->damping();
377 : c *= damp;
378 :
379 : return true;
380 : }
381 :
382 :
383 : template <typename TDomain, typename TAlgebra>
384 : SmartPtr<ILinearIterator<typename TAlgebra::vector_type> >
385 0 : AssembledTransformingSmoother<TDomain, TAlgebra>::
386 : clone()
387 : {
388 0 : SmartPtr<AssembledTransformingSmoother<TDomain, TAlgebra> > clone(
389 0 : new AssembledTransformingSmoother<TDomain, TAlgebra>(
390 : m_spAuxSystemAss, m_spAuxSmoother,
391 : m_spRightTrafoAss, m_spRightTrafoSmoother));
392 :
393 0 : return clone;
394 : }
395 :
396 :
397 : } // end namespace ug
398 :
399 : #endif /* __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__TRANSFORMING_SMOOTHER__ */
|
