Line data Source code
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 :
36 : #include "proj_gauss_seidel_interface.h"
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 :
49 : namespace ug{
50 :
51 : template <typename TDomain, typename TAlgebra>
52 : void
53 : IProjGaussSeidel<TDomain,TAlgebra>::
54 : truncateVec(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 :
70 : template <typename TDomain, typename TAlgebra>
71 : void
72 : IProjGaussSeidel<TDomain,TAlgebra>::
73 : truncateMat(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 :
91 : template <typename TDomain, typename TAlgebra>
92 : bool
93 0 : IProjGaussSeidel<TDomain,TAlgebra>::
94 : init(SmartPtr<ILinearOperator<vector_type> > J, const vector_type& u)
95 : {
96 : PROFILE_FUNC_GROUP("IProjGaussSeidel");
97 :
98 : // call GaussSeidelBase-init
99 0 : base_type::init(J,u);
100 :
101 : // remember solution
102 0 : m_spSol = u.clone();
103 :
104 : // remember, that operator has been initialized
105 0 : 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 0 : for( ; iter != iterEnd; iter++)
114 0 : (*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 0 : return true;
125 : }
126 :
127 : template <typename TDomain, typename TAlgebra>
128 : void
129 0 : IProjGaussSeidel<TDomain,TAlgebra>::
130 : project_correction(value_type& c_i, const size_t i)
131 : {
132 0 : 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 0 : 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 0 : 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 0 : for( ; iter != iterEnd; iter++)
151 : {
152 : // check, if the dof lies in an obstacle subset: if not -> continue!
153 0 : if (!((*iter)->is_obs_dof(dof)))
154 0 : continue;
155 :
156 : //UG_LOG("IS IN OBS SUBSET \n");
157 : dofIsObsDoF = true;
158 :
159 0 : (*iter)->adjust_sol_and_cor((*m_spSol)[i], c_i, dofIsActive, dof);
160 : }
161 :
162 0 : if (dofIsObsDoF && (!dofIsActive))
163 : {
164 : // dof is admissible -> do regular solution update intern of the
165 : // Projected GaussSeidel class
166 0 : BlockRef((*m_spSol)[i], comp) += BlockRef(c_i, comp);
167 : }
168 : }
169 : }
170 :
171 : template <typename TDomain, typename TAlgebra>
172 : bool
173 0 : IProjGaussSeidel<TDomain,TAlgebra>::
174 : apply(vector_type &c, const vector_type& d)
175 : {
176 : PROFILE_FUNC_GROUP("IProjGaussSeidel");
177 : // Check that operator is initialized
178 0 : if(!m_bInit)
179 : {
180 : UG_LOG("ERROR in 'IProjGaussSeidel::apply': Iterator not initialized.\n");
181 0 : return false;
182 : }
183 :
184 : // loop all obstacle constraints, which are set & reset its active dofs
185 0 : 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 0 : for( ; iter != iterEnd; iter++)
192 : (*iter)->reset_active_dofs();
193 : }
194 :
195 0 : base_type::apply(c, d);
196 :
197 : // TODO: in case of using the projected GaussSeidel as smoother in a GMG: TRUNCATION
198 0 : const GF& def = dynamic_cast<const GF&>(d);
199 0 : ConstSmartPtr<GF> spD = def.clone_without_values();
200 0 : if(spD.valid())
201 : {
202 : // check if the DofDistribution is a MGDofDistribution
203 : //if (spD->dof_distribution()->multi_grid().valid())
204 0 : int surfaceLev = spD->dof_distribution()->grid_level().level();
205 :
206 0 : UG_LOG("NumIndices :" <<spD->dof_distribution()->num_indices() << "\n");
207 0 : UG_LOG("numLevels: " << spD->approx_space()->num_levels() << "\n");
208 :
209 0 : if (spD->dof_distribution()->multi_grid().valid())
210 : {
211 0 : size_t topLev = spD->dof_distribution()->multi_grid()->top_level();
212 0 : UG_LOG("surfaceLev: " << surfaceLev << "!\n");
213 : UG_LOG("topLev: " << topLev << "!\n");
214 :
215 0 : if((size_t)surfaceLev == topLev)
216 : {
217 : UG_LOG("topLev gleich surfaceLev!\n");
218 : //TODO: for all obstacle constraints:
219 0 : 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 0 : 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 0 : for( ; dofIter != dofIterEnd; dofIter++)
236 0 : 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 :
295 : template <typename TDomain, typename TAlgebra>
296 : bool
297 0 : IProjGaussSeidel<TDomain,TAlgebra>::
298 : apply_update_defect(vector_type &c, vector_type& d)
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 0 : base_type::apply_update_defect(c, d);
307 :
308 : // adjust defect for the active dofs
309 0 : 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 0 : for( ; iter != iterEnd; iter++)
316 0 : (*iter)->adjust_defect_to_constraint(d);
317 : }
318 :
319 0 : 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__ */
|
