Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
multigrid.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: multigrid.templates.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__multigrid_templates_h
18 #define __deal2__multigrid_templates_h
19 #include <deal.II/multigrid/multigrid.h>
20 
21 #include <deal.II/base/logstream.h>
22 
23 #include <iostream>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 
28 template <class VECTOR>
29 Multigrid<VECTOR>::Multigrid (const unsigned int minlevel,
30  const unsigned int maxlevel,
31  const MGMatrixBase<VECTOR> &matrix,
32  const MGCoarseGridBase<VECTOR> &coarse,
33  const MGTransferBase<VECTOR> &transfer,
34  const MGSmootherBase<VECTOR> &pre_smooth,
35  const MGSmootherBase<VECTOR> &post_smooth,
36  typename Multigrid<VECTOR>::Cycle cycle)
37  :
38  cycle_type(cycle),
39  minlevel(minlevel),
40  maxlevel(maxlevel),
41  defect(minlevel,maxlevel),
42  solution(minlevel,maxlevel),
43  t(minlevel,maxlevel),
44  matrix(&matrix, typeid(*this).name()),
45  coarse(&coarse, typeid(*this).name()),
46  transfer(&transfer, typeid(*this).name()),
47  pre_smooth(&pre_smooth, typeid(*this).name()),
48  post_smooth(&post_smooth, typeid(*this).name()),
49  edge_out(0, typeid(*this).name()),
50  edge_in(0, typeid(*this).name()),
51  edge_down(0, typeid(*this).name()),
52  edge_up(0, typeid(*this).name()),
53  debug(0)
54 {}
55 
56 
57 
58 template <class VECTOR>
59 void
60 Multigrid<VECTOR>::reinit(const unsigned int min_level,
61  const unsigned int max_level)
62 {
63  minlevel=min_level;
64  maxlevel=max_level;
65  defect.resize(minlevel, maxlevel);
66 }
67 
68 
69 template <class VECTOR>
70 void
71 Multigrid<VECTOR>::set_maxlevel (const unsigned int l)
72 {
73  Assert (l <= maxlevel, ExcIndexRange(l,minlevel,maxlevel+1));
74  Assert (l >= minlevel, ExcIndexRange(l,minlevel,maxlevel+1));
75  maxlevel = l;
76 }
77 
78 
79 template <class VECTOR>
80 void
81 Multigrid<VECTOR>::set_minlevel (const unsigned int l,
82  const bool relative)
83 {
84  Assert (l <= maxlevel, ExcIndexRange(l,minlevel,maxlevel+1));
85  minlevel = (relative)
86  ? (maxlevel-l)
87  : l;
88 }
89 
90 
91 template <class VECTOR>
92 void
94 {
95  cycle_type = c;
96 }
97 
98 
99 template <class VECTOR>
100 void
101 Multigrid<VECTOR>::set_debug (const unsigned int d)
102 {
103  debug = d;
104 }
105 
106 
107 template <class VECTOR>
108 void
110  const MGMatrixBase<VECTOR> &up)
111 {
112  edge_out = &down;
113  edge_in = &up;
114 }
115 
116 
117 template <class VECTOR>
118 void
120  const MGMatrixBase<VECTOR> &up)
121 {
122  edge_down = &down;
123  edge_up = &up;
124 }
125 
126 
127 template <class VECTOR>
128 void
129 Multigrid<VECTOR>::level_v_step(const unsigned int level)
130 {
131  if (debug>0)
132  deallog << "V-cycle entering level " << level << std::endl;
133  if (debug>2)
134  deallog << "V-cycle Defect norm " << defect[level].l2_norm()
135  << std::endl;
136 
137  if (level == minlevel)
138  {
139  if (debug>0)
140  deallog << "Coarse level " << level << std::endl;
141  (*coarse)(level, solution[level], defect[level]);
142  return;
143  }
144  if (debug>1)
145  deallog << "Smoothing on level " << level << std::endl;
146  // smoothing of the residual by
147  // modifying s
148 // defect[level].print(std::cout, 2,false);
149 // std::cout<<std::endl;
150  pre_smooth->smooth(level, solution[level], defect[level]);
151 // solution[level].print(std::cout, 2,false);
152 
153  if (debug>2)
154  deallog << "Solution norm " << solution[level].l2_norm()
155  << std::endl;
156 
157  if (debug>1)
158  deallog << "Residual on level " << level << std::endl;
159  // t = A*solution[level]
160  matrix->vmult(level, t[level], solution[level]);
161 
162  if (debug>2)
163  deallog << "Residual norm " << t[level].l2_norm()
164  << std::endl;
165 // std::cout<<std::endl;
166 // t[level].print(std::cout, 2,false);
167 
168  // make t rhs of lower level The
169  // non-refined parts of the
170  // coarse-level defect already
171  // contain the global defect, the
172  // refined parts its restriction.
173  for (unsigned int l = level; l>minlevel; --l)
174  {
175  t[l-1] = 0.;
176  if (l==level && edge_out != 0)
177  {
178  edge_out->vmult_add(level, t[level], solution[level]);
179  if (debug>2)
180  deallog << "Norm t[" << level << "] " << t[level].l2_norm() << std::endl;
181  }
182 
183  if (l==level && edge_down != 0)
184  edge_down->vmult(level, t[level-1], solution[level]);
185 
186  transfer->restrict_and_add (l, t[l-1], t[l]);
187  if (debug>3)
188  deallog << "restrict t[" << l-1 << "] " << t[l-1].l2_norm() << std::endl;
189  defect[l-1] -= t[l-1];
190  if (debug>3)
191  deallog << "defect d[" << l-1 << "] " << defect[l-1].l2_norm() << std::endl;
192  }
193 
194  // do recursion
195  solution[level-1] = 0.;
196  level_v_step(level-1);
197 
198  // reset size of the auxiliary
199  // vector, since it has been
200  // resized in the recursive call to
201  // level_v_step directly above
202  t[level] = 0.;
203 
204  // do coarse grid correction
205  transfer->prolongate(level, t[level], solution[level-1]);
206  if (debug>2)
207  deallog << "Prolongate norm " << t[level].l2_norm() << std::endl;
208 
209  solution[level] += t[level];
210 
211  if (edge_in != 0)
212  {
213  edge_in->Tvmult(level, t[level], solution[level]);
214  defect[level] -= t[level];
215  }
216 
217  if (edge_up != 0)
218  {
219  edge_up->Tvmult(level, t[level], solution[level-1]);
220  defect[level] -= t[level];
221  }
222 
223  if (debug>2)
224  deallog << "V-cycle Defect norm " << defect[level].l2_norm()
225  << std::endl;
226 
227  if (debug>1)
228  deallog << "Smoothing on level " << level << std::endl;
229  // post-smoothing
230 
231 // std::cout<<std::endl;
232 // defect[level].print(std::cout, 2,false);
233  post_smooth->smooth(level, solution[level], defect[level]);
234 // solution[level].print(std::cout, 2,false);
235 // std::cout<<std::endl;
236 
237  if (debug>2)
238  deallog << "Solution norm " << solution[level].l2_norm()
239  << std::endl;
240 
241  if (debug>1)
242  deallog << "V-cycle leaving level " << level << std::endl;
243 }
244 
245 
246 
247 template <class VECTOR>
248 void
249 Multigrid<VECTOR>::level_step(const unsigned int level,
250  Cycle cycle)
251 {
252  char cychar = '?';
253  switch (cycle)
254  {
255  case v_cycle:
256  cychar = 'V';
257  break;
258  case f_cycle:
259  cychar = 'F';
260  break;
261  case w_cycle:
262  cychar = 'W';
263  break;
264  }
265 
266  if (debug>0)
267  deallog << cychar << "-cycle entering level " << level << std::endl;
268 
269  // Not actually the defect yet, but
270  // we do not want to spend yet
271  // another vector.
272  if (level>minlevel)
273  {
274  defect2[level-1] = 0.;
275  transfer->restrict_and_add (level, defect2[level-1], defect2[level]);
276  }
277 
278  // We get an update of the defect
279  // from the previous level in t and
280  // from two levels above in
281  // defect2. This is subtracted from
282  // the original defect.
283  t[level].equ(-1.,defect2[level],1.,defect[level]);
284 
285  if (debug>2)
286  deallog << cychar << "-cycle defect norm " << t[level].l2_norm()
287  << std::endl;
288 
289  if (level == minlevel)
290  {
291  if (debug>0)
292  deallog << cychar << "-cycle coarse level " << level << std::endl;
293 
294  (*coarse)(level, solution[level], t[level]);
295  return;
296  }
297  if (debug>1)
298  deallog << cychar << "-cycle smoothing level " << level << std::endl;
299  // smoothing of the residual by
300  // modifying s
301  pre_smooth->smooth(level, solution[level], t[level]);
302 
303  if (debug>2)
304  deallog << cychar << "-cycle solution norm "
305  << solution[level].l2_norm() << std::endl;
306 
307  if (debug>1)
308  deallog << cychar << "-cycle residual level " << level << std::endl;
309  // t = A*solution[level]
310  matrix->vmult(level, t[level], solution[level]);
311  // make t rhs of lower level The
312  // non-refined parts of the
313  // coarse-level defect already
314  // contain the global defect, the
315  // refined parts its restriction.
316  if (edge_out != 0)
317  edge_out->vmult_add(level, t[level], solution[level]);
318 
319  if (edge_down != 0)
320  edge_down->vmult_add(level, defect2[level-1], solution[level]);
321 
322  transfer->restrict_and_add (level, defect2[level-1], t[level]);
323  // do recursion
324  solution[level-1] = 0.;
325  // Every cycle need one recursion,
326  // the V-cycle, which is included
327  // here for the sake of the
328  // F-cycle, needs only one,
329  level_step(level-1, cycle);
330  // If we solve exactly, then we do
331  // not need a second coarse grid
332  // step.
333  if (level>minlevel+1)
334  {
335  // while the W-cycle repeats itself
336  if (cycle == w_cycle)
337  level_step(level-1, cycle);
338  // and the F-cycle does a V-cycle
339  // after an F-cycle.
340  else if (cycle == f_cycle)
341  level_step(level-1, v_cycle);
342  }
343 
344  // reset size of the auxiliary
345  // vector, since it has been
346  // resized in the recursive call to
347  // level_v_step directly above
348  t[level] = 0.;
349  // do coarse grid correction
350  transfer->prolongate(level, t[level], solution[level-1]);
351 
352  solution[level] += t[level];
353 
354 
355  if (edge_in != 0)
356  edge_in->Tvmult(level, t[level], solution[level]);
357 
358  if (edge_up != 0)
359  edge_up->Tvmult(level, t[level], solution[level-1]);
360 
361  t[level].sadd(-1.,-1.,defect2[level],1.,defect[level]);
362 
363  if (debug>2)
364  deallog << cychar << "-cycle Defect norm " << t[level].l2_norm()
365  << std::endl;
366 
367  if (debug>1)
368  deallog << cychar << "-cycle smoothing level " << level << std::endl;
369  // post-smoothing
370  post_smooth->smooth(level, solution[level], t[level]);
371 
372  if (debug>1)
373  deallog << cychar << "-cycle leaving level " << level << std::endl;
374 }
375 
376 
377 template <class VECTOR>
378 void
380 {
381  // The defect vector has been
382  // initialized by copy_to_mg. Now
383  // adjust the other vectors.
384  solution.resize(minlevel, maxlevel);
385  t.resize(minlevel, maxlevel);
386  if (cycle_type != v_cycle)
387  defect2.resize(minlevel, maxlevel);
388 
389  for (unsigned int level=minlevel; level<=maxlevel; ++level)
390  {
391  solution[level].reinit(defect[level]);
392  t[level].reinit(defect[level]);
393  if (cycle_type != v_cycle)
394  defect2[level].reinit(defect[level]);
395  }
396 
397  if (cycle_type == v_cycle)
398  level_v_step (maxlevel);
399  else
400  level_step (maxlevel, cycle_type);
401 }
402 
403 
404 template <class VECTOR>
405 void
407 {
408  // The defect vector has been
409  // initialized by copy_to_mg. Now
410  // adjust the other vectors.
411  solution.resize(minlevel, maxlevel);
412  t.resize(minlevel, maxlevel);
413 
414  for (unsigned int level=minlevel; level<=maxlevel; ++level)
415  {
416  solution[level].reinit(defect[level]);
417  t[level].reinit(defect[level]);
418  }
419  level_v_step (maxlevel);
420 }
421 
422 
423 template <class VECTOR>
424 void
425 Multigrid<VECTOR>::vmult(VECTOR &dst, const VECTOR &src) const
426 {
427  Multigrid<VECTOR> &mg = const_cast<Multigrid<VECTOR>&>(*this);
428  mg.defect[maxlevel] = src;
429  for (unsigned int level=maxlevel; level>minlevel; --level)
430  {
431  mg.defect[level-1] = 0.;
432  mg.transfer->restrict_and_add (level,
433  mg.defect[level-1],
434  mg.defect[level]);
435  }
436 
437  mg.cycle();
438  dst = mg.solution[maxlevel];
439 }
440 
441 
442 template <class VECTOR>
443 void
444 Multigrid<VECTOR>::vmult_add(VECTOR &dst, const VECTOR &src) const
445 {
446  Multigrid<VECTOR> &mg = const_cast<Multigrid<VECTOR>&>(*this);
447  mg.defect[maxlevel] = src;
448  mg.cycle();
449  dst += mg.solution[maxlevel];
450 }
451 
452 
453 template <class VECTOR>
454 void
455 Multigrid<VECTOR>::Tvmult(VECTOR &, const VECTOR &) const
456 {
457  Assert(false, ExcNotImplemented());
458 }
459 
460 
461 template <class VECTOR>
462 void
463 Multigrid<VECTOR>::Tvmult_add(VECTOR &, const VECTOR &) const
464 {
465  Assert(false, ExcNotImplemented());
466 }
467 
468 DEAL_II_NAMESPACE_CLOSE
469 
470 #endif
Multigrid(const DoFHandler< dim > &mg_dof_handler, const MGMatrixBase< VECTOR > &matrix, const MGCoarseGridBase< VECTOR > &coarse, const MGTransferBase< VECTOR > &transfer, const MGSmootherBase< VECTOR > &pre_smooth, const MGSmootherBase< VECTOR > &post_smooth, Cycle cycle=v_cycle)
void Tvmult_add(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
void set_edge_matrices(const MGMatrixBase< VECTOR > &edge_out, const MGMatrixBase< VECTOR > &edge_in)
MGLevelObject< VECTOR > solution
Definition: multigrid.h:411
void Tvmult(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
void level_step(const unsigned int level, Cycle cycle)
void level_v_step(const unsigned int level)
void vmult_add(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
SmartPointer< const MGTransferBase< VECTOR >, Multigrid< VECTOR > > transfer
Definition: multigrid.h:440
#define Assert(cond, exc)
Definition: exceptions.h:299
void set_cycle(Cycle)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void vmult(VECTOR &dst, const VECTOR &src) const DEAL_II_DEPRECATED
void set_maxlevel(const unsigned int)
void set_minlevel(const unsigned int level, bool relative=false)
void reinit(const unsigned int minlevel, const unsigned int maxlevel)
void set_edge_flux_matrices(const MGMatrixBase< VECTOR > &edge_down, const MGMatrixBase< VECTOR > &edge_up)
::ExceptionBase & ExcNotImplemented()
void set_debug(const unsigned int)
MGLevelObject< VECTOR > defect
Definition: multigrid.h:405