ug4
Loading...
Searching...
No Matches
nl_gauss_seidel_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/*
34 * (main parts are based on the structure of
35 * newton_impl.h and some ideas of Sebastian Reiter & Andreas Vogel)
36 */
37
38#ifndef __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NL_GAUSS_SEIDEL__NL_GAUSS_SEIDEL_IMPL_H_
39#define __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NL_GAUSS_SEIDEL__NL_GAUSS_SEIDEL_IMPL_H_
40
41// extern includes
42#include <iostream>
43
46#include "nl_gauss_seidel.h"
47
48#define PROFILE_NL_GAUSSSEIDEL
49#ifdef PROFILE_NL_GAUSSSEIDEL
50 #define NL_GAUSSSEIDEL_PROFILE_FUNC() PROFILE_FUNC_GROUP("NL GaussSeidel")
51 #define NL_GAUSSSEIDEL_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "NL GaussSeidel")
52 #define NL_GAUSSSEIDEL_PROFILE_END() PROFILE_END()
53#else
54 #define NL_GAUSSSEIDEL_PROFILE_FUNC()
55 #define NL_GAUSSSEIDEL_PROFILE_BEGIN(name)
56 #define NL_GAUSSSEIDEL_PROFILE_END()
57#endif
58
59namespace ug{
60
61
62template <typename TAlgebra>
65{
66 const LocalIndices& ind = lvec.get_indices();
67
68 for(size_t fct=0; fct < lvec.num_all_fct(); ++fct)
69 for(size_t dof=0; dof < lvec.num_all_dof(fct); ++dof)
70 {
71 const size_t index = ind.index(fct,dof);
72 if (index == m_assemblingIndex)
73 {
74 const size_t comp = ind.comp(fct,dof);
75 BlockRef(vec[0], comp) += lvec.value(fct,dof);
76 }
77 }
78
79}
80
81template <typename TAlgebra>
84{
85 const LocalIndices& rowInd = lmat.get_row_indices();
86 const LocalIndices& colInd = lmat.get_col_indices();
87
88 for(size_t fct1=0; fct1 < lmat.num_all_row_fct(); ++fct1)
89 for(size_t dof1=0; dof1 < lmat.num_all_row_dof(fct1); ++dof1)
90 {
91 const size_t rowIndex = rowInd.index(fct1,dof1);
92
93 if (rowIndex == m_assemblingIndex)
94 {
95 const size_t rowComp = rowInd.comp(fct1,dof1);
96 for(size_t fct2=0; fct2 < lmat.num_all_col_fct(); ++fct2)
97 for(size_t dof2=0; dof2 < lmat.num_all_col_dof(fct2); ++dof2)
98 {
99 const size_t colIndex = colInd.index(fct2,dof2);
100 if (colIndex == m_assemblingIndex)
101 {
102 const size_t colComp = colInd.comp(fct2,dof2);
103 BlockRef(mat(0, 0), rowComp, colComp)
104 += lmat.value(fct1,dof1,fct2,dof2);
105 }
106 }
107 }
108 }
109
110}
111
112template <typename TDomain, typename TAlgebra>
116 m_spApproxSpace(spApproxSpace),
117 m_spConvCheck(spConvCheck),
118 m_damp(1.0),
119 m_bProjectedGS(false),
120 m_dgbCall(0)
121{};
122
123template <typename TDomain, typename TAlgebra>
126 m_spApproxSpace(NULL),
127 m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
128 m_damp(1.0),
129 m_bProjectedGS(false),
130 m_dgbCall(0)
131{};
132
133
134template <typename TDomain, typename TAlgebra>
137{
138 m_spConvCheck = spConvCheck;
139 m_spConvCheck->set_offset(3);
140 m_spConvCheck->set_symbol('#');
141 m_spConvCheck->set_name("Nonlinear Gauss Seidel Solver");
142}
143
144template <typename TDomain, typename TAlgebra>
145bool
148{
149 NL_GAUSSSEIDEL_PROFILE_BEGIN(NL_GAUSSSEIDELSolver_init);
150 m_spAssOp = op.template cast_dynamic<AssembledOperator<TAlgebra> >();
151 if(m_spAssOp.invalid())
152 UG_THROW("NLGaussSeidelSolver: currently only works for AssembledDiscreteOperator.");
153
154 m_spAss = m_spAssOp->discretization();
155 if(m_spAss.invalid())
156 UG_THROW("AssembledLinearOperator: Assembling routine not set.");
157
158 // Check for approxSpace
159 if(m_spApproxSpace.invalid())
160 UG_THROW("NLGaussSeidelSolver::prepare: Approximation Space not set.");
161
162 m_gridLevel = m_spAssOp->level();
163
164 // set DoF distribution type
165 if(m_gridLevel.type() == GridLevel::LEVEL)
166 m_spLevDD = m_spApproxSpace->dof_distribution(m_gridLevel);
167 else if (m_gridLevel.type() == GridLevel::SURFACE)
168 m_spSurfDD = m_spApproxSpace->dof_distribution(m_gridLevel);
169 else
170 UG_THROW("Grid Level not recognized.");
171
172 /*TDomain& dom = *m_spApproxSpace->domain();
173 typename TDomain::grid_type& grid = *dom.grid();
174
175 grid.attach_to_vertices(m_aElemList);
176 m_aaElemList.access(grid, m_aElemList);*/
177
178 return true;
179}
180
181
182template <typename TDomain, typename TAlgebra>
184{
185 // In this method a elemList is created for every DoF 'globIndex'.
186 // This element-list incorporates all elements which have contributions
187 // to 'globIndex'.
188
189 // Check for approxSpace
190 if(m_spApproxSpace.invalid())
191 UG_THROW("NLGaussSeidelSolver::apply: Approximation Space not set.");
192
193 // Check for DoF distribution
194 if(m_spLevDD.invalid() && m_spSurfDD.invalid())
195 UG_THROW("NLGaussSeidelSolver::apply: DoFDistribution not set."
196 " 'NLGaussSeidelSolver::init'-call is necessary!");
197
198 // some vars for Element iterators
199 TDomain& dom = *m_spApproxSpace->domain();
200 typename TDomain::grid_type& grid = *dom.grid();
201
203 typedef typename TDomain::grid_type::template traits<grid_base_object>::iterator
204 ElemIter;
205
206 ElemIter iterBegin = grid.template begin<grid_base_object>(grid.top_level());
207 ElemIter iterEnd = grid.template end<grid_base_object>(grid.top_level());
208
209 typedef typename elemList::iterator ListIter;
210
211 /* VERSION �BER ATTACHMENT:
212 typename Grid::traits<grid_base_object>::secure_container elems;
213
214 int vtr_count = 0;
215
216 // loop over all vertices on grid/multigrid
217 for(VertexIterator iter = grid.vertices_begin(); iter != grid.vertices_end(); ++iter)
218 {
219 grid.associated_elements(elems, *iter);
220 for(size_t i = 0; i < elems.size(); ++i) m_aaElemList[*iter].push_back(elems[i]);
221 }*/
222
223 LocalVector locU; LocalIndices ind;
224
225 m_vElemList.resize(u.size());
226
227 // loop over all elements in toplevel of the grid/multigrid
228 for(ElemIter elemIter = iterBegin; elemIter != iterEnd; ++elemIter)
229 {
230 // get element
231 grid_base_object* elem = *elemIter;
232
233 if(m_gridLevel.type() == GridLevel::LEVEL)
234 m_spLevDD->indices(elem, ind);
235 else if (m_gridLevel.type() == GridLevel::SURFACE)
236 m_spSurfDD->indices(elem, ind);
237 else
238 UG_THROW("Grid Level not recognized.");
239
240 // adapt local algebra
241 locU.resize(ind);
242
243 bool elemInList = false;
244
245 for(size_t fct=0; fct < locU.num_all_fct(); ++fct)
246 for(size_t dof=0; dof < locU.num_all_dof(fct); ++dof)
247 {
248 size_t globIndex = ind.index(fct,dof);
249
250 // check if the globIndex-th elemList incorporates the element *iter
251 for(ListIter listIter = m_vElemList[globIndex].begin();
252 listIter != m_vElemList[globIndex].end(); ++listIter)
253 if (*listIter == *elemIter)
254 elemInList = true;
255
256 // only add the element *elemIter to the globIndex-th elemList,
257 // if the list does not incorporate the element already
258 if (!elemInList) m_vElemList[globIndex].push_back(*elemIter);
259
260 } //end (dof)
261 } //end (elem)
262
263 return true;
264}
265
266
267template <typename TDomain, typename TAlgebra>
269{
270 NL_GAUSSSEIDEL_PROFILE_BEGIN(NL_GAUSSSEIDELSolver_apply);
271 // increase call count
272 m_dgbCall++;
273
274 // Check for approxSpace
275 if(m_spApproxSpace.invalid())
276 UG_THROW("NLGaussSeidelSolver::apply: Approximation Space not set.");
277
278 // Check for DoF distribution
279 if(m_spLevDD.invalid() && m_spSurfDD.invalid())
280 UG_THROW("NLGaussSeidelSolver::apply: DoFDistribution not set."
281 " 'NLGaussSeidelSolver::init'-call is necessary!");
282
283 // Jacobian
284 if(m_spJBlock.invalid() || m_spJBlock->discretization() != m_spAss) {
285 m_spJBlock = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
286 m_spJBlock->set_level(m_gridLevel);
287 }
288
289 // create tmp vectors
290 SmartPtr<vector_type> spD = u.clone_without_values();
291 SmartPtr<vector_type> spC = u.clone_without_values();
292 SmartPtr<vector_type> spDBlock = u.clone_without_values();
293
294 // resize vectors, because they are only used as block-vectors with one entry
295 (*spC).resize(1); (*spDBlock).resize(1);
296
297 // Set dirichlet values
298 try{
299 m_spAssOp->prepare(u);
300 }
301 UG_CATCH_THROW("NLGaussSeidelSolver::apply: Prepare of Operator failed.");
302
303 // Compute first Defect d = L(u)
304 try{
305 NL_GAUSSSEIDEL_PROFILE_BEGIN(NL_GAUSSSEIDELComputeDefect1);
306 m_spAssOp->apply(*spD, u);
308 }UG_CATCH_THROW("NLGaussSeidelSolver::apply: "
309 "Computation of Start-Defect failed.");
310
311 // write start defect for debug
312 int loopCnt = 0;
313 char ext[20]; snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
314 std::string name("NLGaussSeidel_Defect");
315 name.append(ext);
316 write_debug(*spD, name.c_str());
317 write_debug(u, "NLGaussSeidel_StartSolution");
318
319 // start convergence check
320 m_spConvCheck->start(*spD);
321
322 // assign selector to grid
323 TDomain& dom = *m_spApproxSpace->domain();
324 typename TDomain::grid_type& grid = *dom.grid();
325 m_sel.assign_grid(grid);
326
327 matrix_type& JBlock = m_spJBlock->get_matrix();
328
329 // loop iteration
330 while(!m_spConvCheck->iteration_ended())
331 {
332 //bool activeSet_changed = false;
333 m_spAss->ass_tuner()->set_mapping(&m_map);
334
335 // loop all indizes
336 for (size_t i = 0; i < u.size(); i++)
337 {
338 // Compute Jacobian J(u) using the updated u-components
339
340 // since we only need J(i,i) and d(i) in every DoF loop,
341 // we access the i-th Element list indicating the neighborhood of DoF i!
342 m_sel.clear();
343 m_sel.select(m_vElemList[i].begin(), m_vElemList[i].end());
344
345 // by passing the selector to the assembling the assemble operators
346 // are build up only by looping over the elements which has been selected
347 m_spAss->ass_tuner()->set_selector(&m_sel);
348
349 // assemble only with respect to index i (causes resizing of matrices/vectors)
350 m_spAss->ass_tuner()->set_single_index_assembling(i);
351 m_map.set_assembling_index(i);
352
353 try{
354 NL_GAUSSSEIDEL_PROFILE_BEGIN(NL_GAUSSSEIDELComputeJacobian);
355 m_spJBlock->init(u);
357 }UG_CATCH_THROW("NLGaussSeidelSolver::apply: "
358 "Initialization of Jacobian failed.");
359
360 // get i-th block of defect d: d(i) =: m_d_block
361 NL_GAUSSSEIDEL_PROFILE_BEGIN(NL_GAUSSSEIDELComputeLastCompDefect);
362 m_spAssOp->apply(*spDBlock, u);
364
365 // get i,i-th block of J: J(i,i)
366 // depending on the AlgebraType J(i,i) is a 1x1, 2x2, 3x3 Matrix
367 // m_c_i = m_damp * d_i /J_ii
368 NL_GAUSSSEIDEL_PROFILE_BEGIN(NL_GAUSSSEIDELInvertJ);
369 InverseMatMult((*spC)[0], m_damp, JBlock(0,0) , (*spDBlock)[0]);
371
372 // update i-th block of solution
373 u[i] -= (*spC)[0];
374
375 // should projected GS be performed?
376 if (m_bProjectedGS)
377 {
378 // call ProjectVectorCorrection
379
380 /*std::vector<DoFIndex> vActiveSet;
381 value_type diff;
382 diff = u[i] - m_ConsVec[i];
383 // get number of unknowns per value_type
384 // (e.g. if CPU == 3 -> nrFcts = 3!)
385 size_t nrFcts = GetSize(diff);
386
387 // loop fcts
388 for (size_t fct = 0; fct < nrFcts; fct++)
389 {
390 bool penetrate = false;
391
392 if (BlockRef(diff,fct) > 0) // i.e.: u > m_ConsVec
393 {
394 bool MultiIndex_is_in_activeSet = false;
395 penetrate = true;
396
397 std::vector<DoFIndex>::iterator iter;
398
399 // adds MultiIndex-pair (i,fct) to vActiveSet
400 for (iter = vActiveSet.begin(); iter != vActiveSet.end(); ++iter)
401 {
402 DoFIndex activeMultiIndex = *iter;
403
404 if (activeMultiIndex[0] == i && activeMultiIndex[1] == fct)
405 MultiIndex_is_in_activeSet = true;
406 // TODO: anstatt 'MultiIndex_is_in_activeSet' hier zu verwenden,
407 // mit continue in loop weiterspringen
408 }
409
410 if (!MultiIndex_is_in_activeSet)
411 {
412 DoFIndex activeMultiIndex(i,fct);
413 vActiveSet.push_back(activeMultiIndex);
414 //activeSet_changed = true;
415 }
416 }
417
418 if (penetrate)
419 BlockRef(u[i],fct) = BlockRef(m_ConsVec[i],fct);
420
421 } //end (fcts)*/
422 } //end(m_bProjectedGS)
423
424 } //end(DoFs)
425
426 // TODO: if active set not changed, iteration converged!
427 //if (!activeSet_changed)
428 // break;
429
430 // set mapping, selector and ass_index to NULL
431 m_spAss->ass_tuner()->set_mapping();
432 m_spAss->ass_tuner()->set_selector();
433 m_spAss->ass_tuner()->disable_single_index_assembling();
434
435 NL_GAUSSSEIDEL_PROFILE_BEGIN(NL_GAUSSSEIDELComputeLastCompDefect);
436 m_spAssOp->prepare(u);
437 m_spAssOp->apply(*spD, u);
439
440 // update counter
441 loopCnt++;
442 snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
443
444 // check convergence
445 m_spConvCheck->update(*spD);
446
447 // write defect for debug
448 std::string name("NLGaussSeidel_Defect"); name.append(ext);
449 write_debug(*spD, name.c_str());
450 }
451
452 return m_spConvCheck->post();
453}
454
455template <typename TDomain, typename TAlgebra>
457{
458// add iter count to name
459 std::string name(filename);
460 char ext[20]; snprintf(ext, sizeof(ext),"_call%03d", m_dgbCall);
461 name.append(ext).append(".vec");
462
463// write
464 base_writer_type::write_debug(vec, name.c_str());
465}
466
467template <typename TDomain, typename TAlgebra>
469{
470// add iter count to name
471 std::string name(filename);
472 char ext[20]; snprintf(ext, sizeof(ext),"_call%03d", m_dgbCall);
473 name.append(ext).append(".mat");
474
475// write
476 base_writer_type::write_debug(mat, name.c_str());
477}
478
479
480}
481
482#endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NL_GAUSS_SEIDEL__NL_GAUSS_SEIDEL_IMPL_H_ */
location name
Definition checkpoint_util.lua:128
Definition smart_pointer.h:296
Definition smart_pointer.h:108
matrix operator based on the assembling of a problem
Definition assembled_linear_operator.h:60
@ SURFACE
Definition grid_level.h:48
@ LEVEL
Definition grid_level.h:48
Definition convergence_check.h:72
describes a mapping X->Y
Definition operator.h:86
Definition local_algebra.h:50
index_type index(size_t fct, size_t dof) const
global algebra index for (fct, dof)
Definition local_algebra.h:139
comp_type comp(size_t fct, size_t dof) const
algebra comp for (fct, dof)
Definition local_algebra.h:153
Definition local_algebra.h:422
size_t num_all_col_fct() const
returns the number of all functions
Definition local_algebra.h:627
const LocalIndices & get_row_indices() const
get current local indices
Definition local_algebra.h:467
size_t num_all_row_fct() const
returns the number of all functions
Definition local_algebra.h:624
size_t num_all_col_dof(size_t fct) const
returns the number of dofs for a function
Definition local_algebra.h:633
const LocalIndices & get_col_indices() const
get current local indices
Definition local_algebra.h:470
size_t num_all_row_dof(size_t fct) const
returns the number of dofs for a function
Definition local_algebra.h:630
number & value(size_t rowFct, size_t rowDoF, size_t colFct, size_t colDoF)
access to coupling (rowFct, rowDoF) x (colFct, colDoF)
Definition local_algebra.h:636
algebra_type::matrix_type matrix_type
Type of algebra matrix.
Definition nl_gauss_seidel.h:58
void add_local_vec_to_global(vector_type &vec, const LocalVector &lvec, ConstSmartPtr< DoFDistribution > dd)
adds a local vector to the global rhs
Definition nl_gauss_seidel_impl.h:63
void add_local_mat_to_global(matrix_type &mat, const LocalMatrix &lmat, ConstSmartPtr< DoFDistribution > dd)
adds a local matrix to the global matrix
Definition nl_gauss_seidel_impl.h:82
algebra_type::vector_type vector_type
Type of algebra vector.
Definition nl_gauss_seidel.h:61
Definition local_algebra.h:198
size_t num_all_fct() const
returns the number of all functions
Definition local_algebra.h:370
const LocalIndices & get_indices() const
get current local indices
Definition local_algebra.h:225
number & value(size_t fct, size_t dof)
access to dof of a fct (unrestricted functions)
Definition local_algebra.h:376
void resize(const LocalIndices &ind)
resize for current local indices
Definition local_algebra.h:214
size_t num_all_dof(size_t fct) const
returns the number of dofs for a function (unrestricted functions)
Definition local_algebra.h:373
virtual bool prepare(vector_type &u)
prepare Operator
Definition nl_gauss_seidel_impl.h:183
virtual bool init(SmartPtr< IOperator< vector_type > > op)
This operator inverts the Operator op: Y -> X.
Definition nl_gauss_seidel_impl.h:147
NLGaussSeidelSolver()
constructor
Definition nl_gauss_seidel_impl.h:125
TAlgebra::vector_type vector_type
Vector type.
Definition nl_gauss_seidel.h:134
TAlgebra::matrix_type matrix_type
Matrix type.
Definition nl_gauss_seidel.h:131
domain_traits< TDomain::dim >::grid_base_object grid_base_object
Type of geometric base object.
Definition nl_gauss_seidel.h:146
virtual bool apply(vector_type &u)
apply Operator, i.e. op^{-1}(0) = u
Definition nl_gauss_seidel_impl.h:268
void set_convergence_check(SmartPtr< IConvergenceCheck< vector_type > > spConvCheck)
Definition nl_gauss_seidel_impl.h:136
void write_debug(const vector_type &vec, const char *filename)
Definition nl_gauss_seidel_impl.h:456
Definition convergence_check.h:179
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
the ug namespace
double & BlockRef(T &vec, size_t i)
Definition blocks.h:66
bool InverseMatMult(number &dest, const double &beta, const TMat &mat, const TVec &vec)
you can implement this function with GetInverse and MatMult
#define NL_GAUSSSEIDEL_PROFILE_BEGIN(name)
Definition nl_gauss_seidel_impl.h:51
#define NL_GAUSSSEIDEL_PROFILE_END()
Definition nl_gauss_seidel_impl.h:52
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
Definition domain_traits.h:53
function ProblemDisc new(problemDesc, dom)