1 | /*
|
---|
2 | * vmg - a versatile multigrid solver
|
---|
3 | * Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
|
---|
4 | *
|
---|
5 | * vmg is free software: you can redistribute it and/or modify
|
---|
6 | * it under the terms of the GNU General Public License as published by
|
---|
7 | * the Free Software Foundation, either version 3 of the License, or
|
---|
8 | * (at your option) any later version.
|
---|
9 | *
|
---|
10 | * vmg is distributed in the hope that it will be useful,
|
---|
11 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
12 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
13 | * GNU General Public License for more details.
|
---|
14 | *
|
---|
15 | * You should have received a copy of the GNU General Public License
|
---|
16 | * along with this program. If not, see <http://www.gnu.org/licenses/>.
|
---|
17 | */
|
---|
18 |
|
---|
19 | /**
|
---|
20 | * @file level_operator_fas.cpp
|
---|
21 | * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
|
---|
22 | * @date Mon Apr 18 13:00:23 2011
|
---|
23 | *
|
---|
24 | * @brief VMG::LevelOperatorFAS
|
---|
25 | *
|
---|
26 | */
|
---|
27 |
|
---|
28 | #ifdef HAVE_CONFIG_H
|
---|
29 | #include <libvmg_config.h>
|
---|
30 | #endif
|
---|
31 |
|
---|
32 | #include <limits>
|
---|
33 |
|
---|
34 | #include "base/discretization.hpp"
|
---|
35 | #include "base/helper.hpp"
|
---|
36 | #include "base/index.hpp"
|
---|
37 | #include "comm/comm.hpp"
|
---|
38 | #include "grid/grid_index_translations.hpp"
|
---|
39 | #include "grid/multigrid.hpp"
|
---|
40 | #include "grid/tempgrid.hpp"
|
---|
41 | #include "level/level_operator_fas.hpp"
|
---|
42 | #include "mg.hpp"
|
---|
43 |
|
---|
44 | using namespace VMG;
|
---|
45 |
|
---|
46 | void LevelOperatorFAS::Restrict(Multigrid& sol, Multigrid& rhs)
|
---|
47 | {
|
---|
48 | Grid::iterator iter_f, iter_c;
|
---|
49 | Stencil::iterator stencil_iter;
|
---|
50 |
|
---|
51 | Comm& comm = *MG::GetComm();
|
---|
52 |
|
---|
53 | Grid& sol_f = sol(sol.Level());
|
---|
54 | Grid& rhs_f = rhs(rhs.Level());
|
---|
55 |
|
---|
56 | Grid& sol_c_dist = sol(sol.Level()-1);
|
---|
57 | Grid& rhs_c_dist = rhs(rhs.Level()-1);
|
---|
58 |
|
---|
59 | Grid& sol_c_undist = comm.GetCoarserGrid(sol);
|
---|
60 | Grid& rhs_c_undist = comm.GetCoarserGrid(rhs);
|
---|
61 |
|
---|
62 | GridIteratorSet bounds_c, bounds_f;
|
---|
63 | GridIndexTranslations::GetGridAlignment(rhs_c_undist, bounds_c, rhs_f, bounds_f);
|
---|
64 |
|
---|
65 | const Stencil& op_res = OperatorRestrict();
|
---|
66 | const Stencil& op_sol = OperatorSol();
|
---|
67 | const Stencil& op_pde = MG::GetDiscretization()->GetStencil();
|
---|
68 | const vmg_float prefactor = MG::GetDiscretization()->OperatorPrefactor(sol_c_undist);
|
---|
69 |
|
---|
70 | // Communication step
|
---|
71 | comm.CommToGhosts(sol_f);
|
---|
72 | comm.CommSubgrid(sol_c_dist, sol_c_undist, 0);
|
---|
73 | comm.CommSubgrid(rhs_c_dist, rhs_c_undist, 0);
|
---|
74 |
|
---|
75 | MG::GetDiscretization()->SetInnerBoundary(sol_f, rhs_f, sol_c_undist);
|
---|
76 |
|
---|
77 | // Compute residual
|
---|
78 | VMG::TempGrid *res_grid = MG::GetTempGrid();
|
---|
79 | res_grid->SetProperties(rhs_f);
|
---|
80 | res_grid->ImportFromResidual(sol_f, rhs_f);
|
---|
81 |
|
---|
82 | for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
|
---|
83 | Helper::AssertVectorsEqual(sol_c_undist.GetSpatialPos(*iter_c), sol_f.GetSpatialPos(*iter_f), 1e6 * std::numeric_limits<vmg_float>::epsilon());
|
---|
84 | sol_c_undist(*iter_c) = op_sol.Apply(sol_f, *iter_f);
|
---|
85 | }
|
---|
86 |
|
---|
87 | comm.CommToGhosts(sol_c_undist);
|
---|
88 |
|
---|
89 | for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
|
---|
90 | Helper::AssertVectorsEqual(rhs_c_undist.GetSpatialPos(*iter_c), (*res_grid).GetSpatialPos(*iter_f), 1e6 * std::numeric_limits<vmg_float>::epsilon());
|
---|
91 | rhs_c_undist(*iter_c) = op_res.Apply(*res_grid, *iter_f) + prefactor * op_pde.Apply(sol_c_undist, *iter_c);
|
---|
92 | }
|
---|
93 |
|
---|
94 | /*
|
---|
95 | if (rhs_c_undist.Global().BoundaryType() == GlobalCoarsened)
|
---|
96 | rhs_c_undist.SetBoundary(sol_c_undist);
|
---|
97 | */
|
---|
98 |
|
---|
99 | comm.CommSubgrid(sol_c_undist, sol_c_dist, 1);
|
---|
100 | comm.CommSubgrid(rhs_c_undist, rhs_c_dist, 1);
|
---|
101 |
|
---|
102 | VMG::TempGrid* sol_old = this->GetTempGrid(sol_c_dist.Level());
|
---|
103 | sol_old->SetProperties(sol_c_dist);
|
---|
104 | sol_old->SetGrid(sol_c_dist);
|
---|
105 |
|
---|
106 | // comm.CommToGhosts(rhs_c_dist);
|
---|
107 |
|
---|
108 | sol.ToCoarserLevel();
|
---|
109 | rhs.ToCoarserLevel();
|
---|
110 | }
|
---|
111 |
|
---|
112 | void LevelOperatorFAS::Prolongate(Multigrid& sol, Multigrid& rhs)
|
---|
113 | {
|
---|
114 | Grid::iterator iter_f, iter_c;
|
---|
115 | Stencil::iterator stencil_iter;
|
---|
116 | vmg_float val;
|
---|
117 |
|
---|
118 | Comm& comm = *MG::GetComm();
|
---|
119 |
|
---|
120 | Grid& sol_c = sol(sol.Level());
|
---|
121 | Grid& sol_f_dist = sol(sol.Level()+1);
|
---|
122 | Grid& rhs_f_dist = rhs(rhs.Level()+1);
|
---|
123 | Grid& sol_f_undist = comm.GetFinerGrid(sol);
|
---|
124 |
|
---|
125 | const Stencil& op = OperatorProlongate();
|
---|
126 |
|
---|
127 | TempGrid *sol_old = this->GetTempGrid(sol_c.Level());
|
---|
128 |
|
---|
129 | GridIteratorSet bounds_c, bounds_f;
|
---|
130 | GridIndexTranslations::GetGridAlignment(sol_c, bounds_c, sol_f_undist, bounds_f);
|
---|
131 |
|
---|
132 | comm.CommSubgrid(sol_f_dist, sol_f_undist, 0);
|
---|
133 | sol_f_undist.ClearHalo();
|
---|
134 |
|
---|
135 | for (iter_f=bounds_f.Begin(), iter_c=bounds_c.Begin(); iter_c!=bounds_c.End(); iter_f+=2, ++iter_c) {
|
---|
136 | Helper::AssertVectorsEqual(sol_c.GetSpatialPos(*iter_c), sol_f_undist.GetSpatialPos(*iter_f), 1e6 * std::numeric_limits<vmg_float>::epsilon());
|
---|
137 | val = sol_c.GetVal(*iter_c) - sol_old->GetVal(*iter_c);
|
---|
138 | sol_f_undist(*iter_f) += op.GetDiag() * val;
|
---|
139 | for (stencil_iter = op.begin(); stencil_iter != op.end(); ++stencil_iter)
|
---|
140 | sol_f_undist(*iter_f + stencil_iter->Disp()) += stencil_iter->Val() * val;
|
---|
141 | }
|
---|
142 |
|
---|
143 | comm.CommFromGhosts(sol_f_undist);
|
---|
144 | comm.CommSubgrid(sol_f_undist, sol_f_dist, 1);
|
---|
145 |
|
---|
146 | if (sol_f_dist.Global().BoundaryType() == LocallyRefined)
|
---|
147 | MG::GetDiscretization()->SetInnerBoundary(sol_f_dist, rhs_f_dist, sol_c);
|
---|
148 |
|
---|
149 | sol.ToFinerLevel();
|
---|
150 | rhs.ToFinerLevel();
|
---|
151 | }
|
---|