ug4
Loading...
Searching...
No Matches
transforming.h
Go to the documentation of this file.
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
39
42
43namespace ug{
44
63template<typename TAlgebra, typename TDerived>
65 public ILinearIterator<typename TAlgebra::vector_type>,
66 public DebugWritingObject<TAlgebra>
67{
68private:
69
71 TDerived& derived()
72 { return *static_cast<TDerived*>(this);}
73
74public:
75
77 typedef TAlgebra algebra_type;
78
80 typedef typename algebra_type::vector_type vector_type;
81
83 typedef typename algebra_type::matrix_type matrix_type;
84
85
88 { UG_THROW("Do you really want to call a copy constructor of this abstract base class?"); }
89
90
91 // Begin CRTP
92
95 { return derived().init(J,u);};
96
99 { return derived().init(L);};
100
103 { return derived().original_operator();}
104
107 { return derived().transform_defect(c,d); }
108
111 { return derived().apply_transformed(c,d); }
112
115 { return derived().untransform_correction(c,d); }
116
117
118 // End CRTP
119
121 virtual bool apply(vector_type& c, const vector_type& d)
122 {
123 return (derived().transform_defect(c,d) &&
124 derived().apply_transformed(c,d) &&
126 }
127
128
131 {
132 // compute correction
133 if (!apply(c,d)) return false;
134
135 // compute updated defect
136 original_operator()->apply_sub(d, c);
137
138 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 * */
173template <typename TDomain, typename TAlgebra>
175 //public ILinearIterator<typename TAlgebra::vector_type>,
176 //public DebugWritingObject<TAlgebra>
177 public ITransformingIteration<TAlgebra, AssembledTransformingSmoother<TDomain,TAlgebra> >
178{
179 public:
181 typedef TDomain domain_type;
182
184 typedef TAlgebra algebra_type;
185
187 typedef typename algebra_type::vector_type vector_type;
188
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:
208
210 virtual const char* name() const {return "Assembled Transform Smoother";}
211
213 virtual bool supports_parallel() const
214 {return (m_spAuxSmoother.valid()) ? m_spAuxSmoother->supports_parallel() : false;}
215
216 // Begin CRTP functions
218
221 { 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);
226
227
230
231
232
233
236
237 protected:
238
241
244
247
250
253
254
257
260
263
264
265
266
267};
268
269
270
271template <typename TDomain, typename TAlgebra>
274{
275 //m_spOriginalSystemOp = J;
276
277 GridLevel gridLevel;
278
280 dynamic_cast<const GridFunction<TDomain, TAlgebra>*>(&u);
281 if(pSol){
282 gridLevel = pSol->dof_distribution()->grid_level();
283 }
284
285 // init aux operator
286 try{
287 m_spAuxSystemAss->assemble_jacobian(*m_spAuxMatrixOperator, u, gridLevel);
288 }
289 UG_CATCH_THROW("AssembledTransformingSmoother: "
290 " Cannot assemble transformed system matrix.");
291
292 write_debug(*m_spAuxMatrixOperator, "TrafoSystem");
293
294 if(!m_spAuxSmoother->init(m_spAuxMatrixOperator, u))
295 UG_THROW("AssembledTransformingSmoother: "
296 " Cannot init smoother for transformed system matrix.");
297
298
299 // init right transform
300 try{
301 m_spRightTrafoAss->assemble_jacobian(*m_spRightTrafoMat, u, gridLevel);
302 }
303 UG_CATCH_THROW("AssembledTransformingSmoother: "
304 " Cannot assemble right-transformation matrix.");
305
306 write_debug(*m_spRightTrafoMat, "RightTrafo");
307
308
309 return true;
310};
311
312
313
314
318template <typename TDomain, typename TAlgebra>
321{
334 // obtain aux correction vector
335 m_spAuxCorrection = c.clone_without_values();
336
337 // apply smoother of transformed system
338 if(!m_spAuxSmoother->apply(*m_spAuxCorrection, d))
339 {
340 UG_LOG("AssembledTransformingSmoother: Smoother applied incorrectly.\n");
341 return false;
342 }
343
344 return true;
345}
346
347template <typename TDomain, typename TAlgebra>
350{
351
352 // apply right-transformation
353 if (m_spRightTrafoSmoother.valid())
354 {
355 // e.g., for Wittum-type smoothers
356 if(!m_spRightTrafoSmoother->apply(*m_spAuxCorrection, d))
357 {
358 UG_LOG("AssembledTransformingSmoother: Right Trafo smoother failed!\n");
359 return false;
360 }
361 }
362 else
363 {
364 // e.g., for Brandt-Dinar-DGS-type smoothers
365 m_spRightTrafoMat->apply(c, *m_spAuxCorrection);
366 }
367
368
369 // release auy correction vecto
370 m_spAuxCorrection = SPNULL;
371
372 #ifdef UG_PARALLEL
373 c.change_storage_type(PST_CONSISTENT);
374 #endif
375
376 const number damp = this->damping()->damping();
377 c *= damp;
378
379 return true;
380}
381
382
383template <typename TDomain, typename TAlgebra>
386clone()
387{
390 m_spAuxSystemAss, m_spAuxSmoother,
391 m_spRightTrafoAss, m_spRightTrafoSmoother));
392
393 return clone;
394}
395
396
397} // end namespace ug
398
399#endif /* __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__TRANSFORMING_SMOOTHER__ */
Definition smart_pointer.h:108
Definition transforming.h:178
SmartPtr< ILinearOperator< vector_type > > m_spOriginalSystemOp
matrix for original system
Definition transforming.h:240
SmartPtr< IAssemble< TAlgebra > > m_spRightTrafoAss
assembling the right-transformation
Definition transforming.h:256
SmartPtr< ILinearIterator< vector_type > > clone()
Clone.
Definition transforming.h:386
SmartPtr< ILinearIterator< vector_type > > m_spRightTrafoSmoother
smoother used on transformed system
Definition transforming.h:262
SmartPtr< IAssemble< TAlgebra > > m_spAuxSystemAss
assembling the transformed system
Definition transforming.h:246
virtual const char * name() const
name
Definition transforming.h:210
algebra_type::matrix_type matrix_type
Matrix type.
Definition transforming.h:190
AssembledTransformingSmoother(SmartPtr< IAssemble< TAlgebra > > spAuxSystemAss, SmartPtr< ILinearIterator< vector_type > > spAuxSmoother, SmartPtr< IAssemble< TAlgebra > > spRightTrafoAss, SmartPtr< ILinearIterator< vector_type > > spRightTrafoSmoother=SPNULL)
constructor setting approximation space
Definition transforming.h:197
bool transform_defect(vector_type &c, const vector_type &d)
Definition transforming.h:223
SmartPtr< vector_type > m_spAuxCorrection
auxiliary correction vector
Definition transforming.h:243
bool init(SmartPtr< ILinearOperator< vector_type > > J, const vector_type &u)
Definition transforming.h:273
bool init(SmartPtr< ILinearOperator< vector_type > > L)
Since we need grid information, linear operators are not supported...
Definition transforming.h:220
bool apply_transformed(vector_type &c, const vector_type &d)
Definition transforming.h:320
SmartPtr< ILinearIterator< vector_type > > m_spAuxSmoother
smoother used on transformed system
Definition transforming.h:252
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition transforming.h:213
SmartPtr< ILinearOperator< vector_type > > original_operator()
Definition transforming.h:234
bool untransform_correction(vector_type &c, const vector_type &d)
Definition transforming.h:349
TAlgebra algebra_type
Algebra type.
Definition transforming.h:184
SmartPtr< MatrixOperator< matrix_type, vector_type > > m_spRightTrafoMat
matrix (operator) of right-transformation
Definition transforming.h:259
algebra_type::vector_type vector_type
Vector type.
Definition transforming.h:187
TDomain domain_type
Domain.
Definition transforming.h:181
SmartPtr< MatrixOperator< matrix_type, vector_type > > m_spAuxMatrixOperator
matrix (operator) of transformed system
Definition transforming.h:249
Definition debug_writer.h:354
void write_debug(const matrix_type &mat, const char *filename)
write debug output for a matrix (if debug writer set)
Definition debug_writer.h:399
represents numerical solutions on a grid using an algebraic vector
Definition grid_function.h:121
SmartPtr< DoFDistribution > dof_distribution()
Definition grid_function.h:236
Definition grid_level.h:42
Interface providing Jacobian and Defect of a discretization.
Definition assemble_interface.h:110
describes a linear iterator
Definition linear_iterator.h:81
SmartPtr< IDamping< TAlgebra::vector_type, TAlgebra::vector_type > > damping()
returns the scaling
Definition linear_iterator.h:173
describes a linear mapping X->Y
Definition linear_operator.h:80
Definition transforming.h:67
virtual bool apply_update_defect(vector_type &c, vector_type &d)
implementation of apply (final, non-CRTP)
Definition transforming.h:130
algebra_type::matrix_type matrix_type
Matrix type.
Definition transforming.h:83
bool init(SmartPtr< ILinearOperator< vector_type > > L)
implementation of init for linear (CRTP)
Definition transforming.h:98
virtual bool apply(vector_type &c, const vector_type &d)
implementation of apply (final, non-CRTP)
Definition transforming.h:121
bool apply_transformed(vector_type &c, const vector_type &d)
map: dtilde -> ctilde (CRTP)
Definition transforming.h:110
SmartPtr< ILinearOperator< vector_type > > original_operator()
original operator (CRTP)
Definition transforming.h:102
TAlgebra algebra_type
Algebra type.
Definition transforming.h:77
ITransformingIteration()
Definition transforming.h:86
bool init(SmartPtr< ILinearOperator< vector_type > > J, const vector_type &u)
implementation of init for non-linear (CRTP)
Definition transforming.h:94
TDerived & derived()
CRTP operator.
Definition transforming.h:71
bool transform_defect(vector_type &c, const vector_type &d)
map: d -> dtilde (CRTP)
Definition transforming.h:106
algebra_type::vector_type vector_type
Vector type.
Definition transforming.h:80
bool untransform_correction(vector_type &c, const vector_type &d)
map: ctilde -> c (CRTP)
Definition transforming.h:114
ITransformingIteration(const ITransformingIteration &ti)
Definition transforming.h:87
Definition matrix_operator.h:49
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
virtual void init()
const NullSmartPtr SPNULL
The equivalent to NULL for smart pointers.
Definition smart_pointer.h:90
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
the ug namespace
function ProblemDisc new(problemDesc, dom)