ug4
Loading...
Searching...
No Matches
proj_gauss_seidel_interface_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 * Author: Raphael Prohl
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__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_INTERFACE_IMPL__
34#define __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_INTERFACE_IMPL__
35
37
38#define PROFILE_PROJ_GAUSS_SEIDEL
39#ifdef PROFILE_PROJ_GAUSS_SEIDEL
40 #define PROJ_GAUSS_SEIDEL_PROFILE_FUNC() PROFILE_FUNC()
41 #define PROJ_GAUSS_SEIDEL_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "IProjGaussSeidel")
42 #define PROJ_GAUSS_SEIDEL_PROFILE_END() PROFILE_END()
43#else
44 #define PROJ_GAUSS_SEIDEL_PROFILE_FUNC()
45 #define PROJ_GAUSS_SEIDEL_PROFILE_BEGIN(name)
46 #define PROJ_GAUSS_SEIDEL_PROFILE_END()
47#endif
48
49namespace ug{
50
51template <typename TDomain, typename TAlgebra>
52void
54truncateVec(vector_type& vec, vector<DoFIndex>& vInd)
55{
56 typedef typename vector<DoFIndex>::iterator iter_type;
57 iter_type dofIter = vInd.begin();
58 iter_type dofIterEnd = vInd.end();
59 for( ; dofIter != dofIterEnd; dofIter++)
60 {
61 if (vec.size() <= (*dofIter)[0])
62 UG_THROW("vec size is to small in IProjGaussSeidel::truncateVec \n");
63
64 UG_LOG("truncateVec: " <<*dofIter<<"\n");
65 //vec[(*dofIter)[0]] = 0.0;
66 DoFRef(vec, *dofIter) = 0.0;
67 }
68}
69
70template <typename TDomain, typename TAlgebra>
71void
73truncateMat(matrix_type& mat, vector<DoFIndex>& vInd)
74{
75 typedef typename vector<DoFIndex>::iterator iter_type;
76 iter_type dofIter = vInd.begin();
77 iter_type dofIterEnd = vInd.end();
78 for( ; dofIter != dofIterEnd; dofIter++)
79 {
80 UG_LOG("activeDof : " <<*dofIter<< "\n");
81
82 // set row to zero (for dof '(*dofIter)[0]' and its first function)
83 SetRow(mat, (*dofIter)[0], 0.0); //(*dofIter)[1]);
84 // set col to zero
85 //BlockRef(mat, , (*dofIter)[0]) = 0.0;
86 UG_LOG("truncateMat: mat(" <<(*dofIter)[1]<<","<<(*dofIter)[0]<<") \n");
87 mat((*dofIter)[1], (*dofIter)[0]) = 0.0;
88 }
89}
90
91template <typename TDomain, typename TAlgebra>
92bool
95{
96 PROFILE_FUNC_GROUP("IProjGaussSeidel");
97
98// call GaussSeidelBase-init
99 base_type::init(J,u);
100
101// remember solution
102 m_spSol = u.clone();
103
104// remember, that operator has been initialized
105 m_bInit = true;
106
107 UG_LOG("In IProjGaussSeidel::init u hat "<<(*m_spSol).size()<<"Eintraege \n");
108 UG_LOG("\n");
109
110 typedef typename vector<SmartPtr<IObstacleConstraint<TDomain,TAlgebra> > >::iterator iter_type;
111 iter_type iter = m_spvObsConstraint.begin();
112 iter_type iterEnd = m_spvObsConstraint.end();
113 for( ; iter != iterEnd; iter++)
114 (*iter)->preprocess();
115
116// (ugly) hint, that usual damping (x += damp * c) does not make sense for the projected
117// GaussSeidel-method.
118/* const number kappa = this->damping()->damping(u, u, m_spOperator);
119 if(kappa != 1.0){
120 UG_THROW("IProjGaussSeidel::set_damp': Ususal damping is not possible "
121 "for IProjGaussSeidel! Use 'set_sor_relax' instead!");
122 }*/
123
124 return true;
125}
126
127template <typename TDomain, typename TAlgebra>
128void
130project_correction(value_type& c_i, const size_t i)
131{
132 if(!m_bObsCons)
133 return;
134
135 typedef typename vector<SmartPtr<IObstacleConstraint<TDomain,TAlgebra> > >::iterator iter_type;
136 iter_type iterEnd = m_spvObsConstraint.end();
137
138 for(size_t comp = 0; comp < GetSize(c_i); comp++)
139 {
140 DoFIndex dof = DoFIndex(i, comp);
141
142 // loop all obstacle constraint, which are set
143 // & perform a projection: check whether the temporary solution u_{s-1/2}
144 // fulfills the underlying constraint(s) or not
145 bool dofIsActive = false;
146 bool dofIsObsDoF = false;
147 //UG_LOG("dof "<<dof<<"\n");
148 // set iterator to the first obstacle constraint
149 iter_type iter = m_spvObsConstraint.begin();
150 for( ; iter != iterEnd; iter++)
151 {
152 // check, if the dof lies in an obstacle subset: if not -> continue!
153 if (!((*iter)->is_obs_dof(dof)))
154 continue;
155
156 //UG_LOG("IS IN OBS SUBSET \n");
157 dofIsObsDoF = true;
158
159 (*iter)->adjust_sol_and_cor((*m_spSol)[i], c_i, dofIsActive, dof);
160 }
161
162 if (dofIsObsDoF && (!dofIsActive))
163 {
164 // dof is admissible -> do regular solution update intern of the
165 // Projected GaussSeidel class
166 BlockRef((*m_spSol)[i], comp) += BlockRef(c_i, comp);
167 }
168 }
169}
170
171template <typename TDomain, typename TAlgebra>
172bool
174apply(vector_type &c, const vector_type& d)
175{
176 PROFILE_FUNC_GROUP("IProjGaussSeidel");
177// Check that operator is initialized
178 if(!m_bInit)
179 {
180 UG_LOG("ERROR in 'IProjGaussSeidel::apply': Iterator not initialized.\n");
181 return false;
182 }
183
184 // loop all obstacle constraints, which are set & reset its active dofs
185 if(m_bObsCons)
186 {
187 typedef typename vector<SmartPtr<IObstacleConstraint<TDomain,TAlgebra> > >::iterator iter_type;
188 iter_type iter = m_spvObsConstraint.begin();
189 iter_type iterEnd = m_spvObsConstraint.end();
190
191 for( ; iter != iterEnd; iter++)
192 (*iter)->reset_active_dofs();
193 }
194
195 base_type::apply(c, d);
196
197 // TODO: in case of using the projected GaussSeidel as smoother in a GMG: TRUNCATION
198 const GF& def = dynamic_cast<const GF&>(d);
199 ConstSmartPtr<GF> spD = def.clone_without_values();
200 if(spD.valid())
201 {
202 // check if the DofDistribution is a MGDofDistribution
203 //if (spD->dof_distribution()->multi_grid().valid())
204 int surfaceLev = spD->dof_distribution()->grid_level().level();
205
206 UG_LOG("NumIndices :" <<spD->dof_distribution()->num_indices() << "\n");
207 UG_LOG("numLevels: " << spD->approx_space()->num_levels() << "\n");
208
209 if (spD->dof_distribution()->multi_grid().valid())
210 {
211 size_t topLev = spD->dof_distribution()->multi_grid()->top_level();
212 UG_LOG("surfaceLev: " << surfaceLev << "!\n");
213 UG_LOG("topLev: " << topLev << "!\n");
214
215 if((size_t)surfaceLev == topLev)
216 {
217 UG_LOG("topLev gleich surfaceLev!\n");
218 //TODO: for all obstacle constraints:
219 if(m_bObsCons)
220 {
221 typedef typename vector<SmartPtr<IObstacleConstraint<TDomain,TAlgebra> > >::iterator iter_type;
222 iter_type iter = m_spvObsConstraint.begin();
223 iter_type iterEnd = m_spvObsConstraint.end();
224
225 for( ; iter != iterEnd; iter++)
226 {
227 // 1. get all active indices
228 vector<DoFIndex> vActiveDoFs;
229 (*iter)->active_dofs(vActiveDoFs);
230 UG_LOG("vActiveDoFs.size() : " <<vActiveDoFs.size()<< "\n");
231
232 typedef typename vector<DoFIndex>::iterator dof_iter_type;
233 dof_iter_type dofIter = vActiveDoFs.begin();
234 dof_iter_type dofIterEnd = vActiveDoFs.end();
235 for( ; dofIter != dofIterEnd; dofIter++)
236 UG_LOG("activeDof : " <<*dofIter<< "\n");
237
238 /*UG_LOG("\n");
239 // 2. truncation call!
240 vector<DoFIndex> testVec;
241 DoFIndex dofIndex1(0,1); DoFIndex dofIndex2(1,0); DoFIndex dofIndex3(2,1);
242 testVec.push_back(dofIndex1);
243 testVec.push_back(dofIndex2);
244 testVec.push_back(dofIndex3);
245
246 vector_type& d_top = const_cast<vector_type&>(d);
247 d_top[0] = 1.0; d_top[1] = 3.0; d_top[2] = 5.0;
248 UG_LOG("d_top(0): "<<d_top[0]<< "\n");
249 UG_LOG("d_top(1): "<<d_top[1]<< "\n");
250 UG_LOG("d_top(2): "<<d_top[2]<< "\n");
251 UG_LOG("d_top(3): "<<d_top[3]<< "\n");
252 truncateVec(d_top, testVec);
253 UG_LOG("d_top(0): "<<d_top[0]<< "\n");
254 UG_LOG("d_top(1): "<<d_top[1]<< "\n");
255 UG_LOG("d_top(2): "<<d_top[2]<< "\n");
256 UG_LOG("d_top(3): "<<d_top[3]<< "\n");
257
258 matrix_type mat_top;
259 mat_top.resize_and_clear(4, 4);
260 UG_LOG("#rows(mat): "<<mat_top.num_rows()<<"\n");
261 for(size_t i=0; i < mat_top.num_rows(); i++)
262 {
263 size_t num_connections = mat_top.num_connections(i);
264 UG_LOG("#connections in "<<i<<"-th row: "<<num_connections<<"\n");
265 }
266 mat_top(0,0) = 1.0; mat_top(0,1) = 1.5; mat_top(1,1) = 3.0; mat_top(2,2) = 5.0;
267 for(size_t i=0; i < mat_top.num_rows(); i++)
268 for(typename matrix_type::row_iterator conn = mat_top.begin_row(i); conn != mat_top.end_row(i); ++conn)
269 {
270 size_t num_connections = mat_top.num_connections(i);
271 UG_LOG("#connections in "<<i<<"-th row: "<<num_connections<<"\n");
272 UG_LOG("mat: "<<i<<"-th row: "<< conn.value() << "\n");
273 }
274
275 truncateMat(mat_top, testVec);
276
277 for(size_t i=0; i < mat_top.num_rows(); i++)
278 for(typename matrix_type::row_iterator conn = mat_top.begin_row(i); conn != mat_top.end_row(i); ++conn)
279 {
280 UG_LOG("mat: "<<i<<"-th row: "<< conn.value() << "\n");
281 }
282 UG_LOG("\n");*/
283
284 }
285 } //end(if(m_bObsCons))
286
287
288 }
289 }
290 }
291
292 return true;
293}
294
295template <typename TDomain, typename TAlgebra>
296bool
299{
300 PROFILE_FUNC_GROUP("IProjGaussSeidel");
301
302 // by calling 'apply_update_defect' the projected Gauss Seidel preconditioner cannot be
303 // a smoother within a multigrid method
304 //bIsASmoother = false;
305
306 base_type::apply_update_defect(c, d);
307
308 // adjust defect for the active dofs
309 if(m_bObsCons)
310 {
311 typedef typename vector<SmartPtr<IObstacleConstraint<TDomain,TAlgebra> > >::iterator iter_type;
312 iter_type iter = m_spvObsConstraint.begin();
313 iter_type iterEnd = m_spvObsConstraint.end();
314
315 for( ; iter != iterEnd; iter++)
316 (*iter)->adjust_defect_to_constraint(d);
317 }
318
319 return true;
320}
321
322
323} // end namespace ug
324
325#endif /* __H__UG__LIB_ALGEBRA__OPERATOR__PRECONDITIONER__PROJECTED_GAUSS_SEIDEL__PROJ_GAUSS_SEIDEL_INTERFACE_IMPL__ */
Definition smart_pointer.h:296
bool valid() const
returns true if the pointer is valid, false if not.
Definition smart_pointer.h:414
Definition smart_pointer.h:108
describes a linear mapping X->Y
Definition linear_operator.h:80
vector_type::value_type value_type
Value type.
Definition proj_gauss_seidel_interface.h:89
virtual bool init(SmartPtr< ILinearOperator< vector_type > > J, const vector_type &u)
Prepare for Operator J(u) and linearization point u (current solution)
Definition proj_gauss_seidel_interface_impl.h:94
void truncateMat(matrix_type &mat, vector< DoFIndex > &vInd)
Definition proj_gauss_seidel_interface_impl.h:73
void project_correction(value_type &c_i, const size_t i)
projects the correction on the underlying constraints set by the obstacleConstraints
Definition proj_gauss_seidel_interface_impl.h:130
GridFunction< TDomain, TAlgebra > GF
Grid Function type.
Definition proj_gauss_seidel_interface.h:92
algebra_type::matrix_type matrix_type
Matrix type.
Definition proj_gauss_seidel_interface.h:83
algebra_type::vector_type vector_type
Vector type.
Definition proj_gauss_seidel_interface.h:86
virtual bool apply_update_defect(vector_type &c, vector_type &d)
Compute new correction c = B*d and return new defect d := d - A*c.
Definition proj_gauss_seidel_interface_impl.h:298
void truncateVec(vector_type &vec, vector< DoFIndex > &vInd)
Definition proj_gauss_seidel_interface_impl.h:54
virtual bool apply(vector_type &c, const vector_type &d)
Compute new correction c = B*d.
Definition proj_gauss_seidel_interface_impl.h:174
Definition multi_index.h:50
void SetRow(TSparseMatrix &A, size_t i, size_t alpha, number val=0.0)
Definition sparsematrix_util.h:710
#define UG_THROW(msg)
Definition error.h:57
#define UG_LOG(msg)
Definition log.h:367
the ug namespace
double & BlockRef(T &vec, size_t i)
Definition blocks.h:66
MultiIndex< 2 > DoFIndex
type of DoF-Index used to identify an DoF in the Algebra
Definition multi_index.h:272
size_t GetSize(const T &t)
number & DoFRef(TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
Definition multi_index.h:276
#define PROFILE_FUNC_GROUP(groups)
Definition profiler.h:258