OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
FMGPoissonSolver.cpp
Go to the documentation of this file.
1//
2// Class FMGPoissonSolver
3// Interface to the Fortran based AMR Poisson multigrid solver of BoxLib.
4// The Fortran solver is part of AMReX till version 18.07
5// (Link: https://github.com/AMReX-Codes/amrex/archive/18.07.tar.gz)
6//
7// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institute, Villigen PSI, Switzerland
8// All rights reserved.
9//
10// Implemented as part of the PhD thesis
11// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
12//
13// This file is part of OPAL.
14//
15// OPAL is free software: you can redistribute it and/or modify
16// it under the terms of the GNU General Public License as published by
17// the Free Software Foundation, either version 3 of the License, or
18// (at your option) any later version.
19//
20// You should have received a copy of the GNU General Public License
21// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
22//
23#include "FMGPoissonSolver.h"
24
25#include "Utility/PAssert.h"
26
28
29#include <AMReX_ParmParse.H>
30#include <AMReX_Interpolater.H>
31
33 : AmrPoissonSolver<AmrBoxLib>(itsAmrObject_p),
34 reltol_m(1.0e-15),
35 abstol_m(1.0e-9)
36{
37 // Dirichlet boundary conditions are default
38 for (int d = 0; d < AMREX_SPACEDIM; ++d) {
39 bc_m[2 * d] = MGT_BC_DIR;
40 bc_m[2 * d + 1] = MGT_BC_DIR;
41 }
42
43 this->initParameters_m();
44}
45
49 unsigned short baseLevel,
50 unsigned short finestLevel,
51 bool /*prevAsGuess*/)
52{
53 const GeomContainer_t& geom = itsAmrObject_mp->Geom();
54
55 if (AmrGeometry_t::isAllPeriodic()) {
56 for (int d = 0; d < AMREX_SPACEDIM; ++d) {
57 bc_m[2 * d] = MGT_BC_PER;
58 bc_m[2 * d + 1] = MGT_BC_PER;
59 }
60 } else if ( AmrGeometry_t::isAnyPeriodic() ) {
61 for (int d = 0; d < AMREX_SPACEDIM; ++d) {
62 if ( AmrGeometry_t::isPeriodic(d) ) {
63 bc_m[2 * d] = MGT_BC_PER;
64 bc_m[2 * d + 1] = MGT_BC_PER;
65 }
66 }
67 }
68
69 amrex::Vector< AmrScalarFieldContainer_t > grad_phi_edge(rho.size());
70
71 amrex::Vector< AmrField_t > tmp(rho.size());
72
73 PAssert(baseLevel <= finestLevel);
74 PAssert(finestLevel < (int)rho.size());
75
76 for (int lev = baseLevel; lev <= finestLevel ; ++lev) {
77 const AmrProcMap_t& dmap = rho[lev]->DistributionMap();
78 grad_phi_edge[lev].resize(AMREX_SPACEDIM);
79 tmp[lev].define(rho[lev]->boxArray(), dmap, AMREX_SPACEDIM, 1);
80
81 for (int n = 0; n < AMREX_SPACEDIM; ++n) {
82 AmrGrid_t ba = rho[lev]->boxArray();
83 grad_phi_edge[lev][n].reset(new AmrField_t(ba.surroundingNodes(n), dmap, 1, 1));
84 grad_phi_edge[lev][n]->setVal(0.0, 1);
85 }
86 }
87
88 for (int i = 0; i <= finestLevel; ++i) {
89 phi[i]->setVal(0.0, 1);
90 tmp[i].setVal(0.0, 1);
91 }
92
93 double residNorm = this->solveWithF90_m(amrex::GetVecOfPtrs(rho),
94 amrex::GetVecOfPtrs(phi),
95 amrex::GetVecOfVecOfPtrs(grad_phi_edge),
96 geom,
97 baseLevel,
98 finestLevel);
99
100 if ( residNorm > abstol_m ) {
101 std::stringstream ss;
102 ss << "Residual norm: " << std::setprecision(16) << residNorm
103 << " > " << abstol_m << " (absolute tolerance)";
104 throw OpalException("FMGPoissonSolver::solve()",
105 "Multigrid solver did not converge. " + ss.str());
106 }
107
108 for (int lev = baseLevel; lev <= finestLevel; ++lev) {
109 amrex::average_face_to_cellcenter(tmp[lev],
110 amrex::GetVecOfConstPtrs(grad_phi_edge[lev]),
111 geom[lev]);
112
113 for (int j = 0; j < AMREX_SPACEDIM; ++j) {
114 efield[lev][j]->setVal(0.0, 1);
115 amr::AmrField_t::Copy(*(efield[lev][j].get()),
116 tmp[lev], j, 0, 1, 1);
117 efield[lev][j]->FillBoundary(0, 1, geom[lev].periodicity());
118 // we need also minus sign due to \vec{E} = - \nabla\phi
119 efield[lev][j]->mult(-1.0, 0, 1);
120 }
121 }
122}
123
124
125double FMGPoissonSolver::getXRangeMin(unsigned short level) {
126 return itsAmrObject_mp->Geom(level).ProbLo(0);
127}
128
129
130double FMGPoissonSolver::getXRangeMax(unsigned short level) {
131 return itsAmrObject_mp->Geom(level).ProbHi(0);
132}
133
134
135double FMGPoissonSolver::getYRangeMin(unsigned short level) {
136 return itsAmrObject_mp->Geom(level).ProbLo(1);
137}
138
139
140double FMGPoissonSolver::getYRangeMax(unsigned short level) {
141 return itsAmrObject_mp->Geom(level).ProbHi(1);
142}
143
144
145double FMGPoissonSolver::getZRangeMin(unsigned short level) {
146 return itsAmrObject_mp->Geom(level).ProbLo(2);
147}
148
149
150double FMGPoissonSolver::getZRangeMax(unsigned short level) {
151 return itsAmrObject_mp->Geom(level).ProbHi(2);
152}
153
154
156 os << "* ************* F M G P o i s s o n S o l v e r ************************************ " << endl
157 << "* relative tolerance " << reltol_m << '\n'
158 << "* absolute tolerance " << abstol_m << '\n' << endl
159 << "* ******************************************************************** " << endl;
160 return os;
161}
162
163
165 /* The information about the parameters is copied from the
166 * BoxLib documentation chapter 5 on linear solvers and from
167 * the code itself.
168 *
169 * Some paramters that have to be given to the solver using
170 * AMReX_ParmParse.
171 * "doc" tells how the variable is called in
172 * - amrex/Src/LinearSolvers/F_MG/cc_mg_cpp.f90
173 * - amrex/Src/LinearSolvers/F_MG/mg_tower.f90 (conains defaults)
174 * - amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
175 *
176 */
177
178 amrex::ParmParse pp_mg("mg");
179
180 // maximum number of multigrid cycles (doc: max_iter)
181 pp_mg.add("maxiter", 200);
182
183 // maximum number of iterations for the bottom solver (doc: bottom_max_iter)
184 pp_mg.add("maxiter_b", 200);
185
186 // number of smoothings at each level on the way DOWN the V-cycle (doc: nu1)
187 pp_mg.add("nu_1", 2);
188
189 // number of smoothings at each level on the way UP the V-cycle (doc: nu2)
190 pp_mg.add("nu_2", 2);
191
192 // number of smoothing before and after the bottom solver (doc: nub)
193 pp_mg.add("nu_b", 0);
194
195 // number of smoothing ... ?
196 pp_mg.add("nu_f", 8);
197
198 // verbosity of the multigrid solver. Higher numbers give more verbosity (doc: verbose)
199 pp_mg.add("v" , 0);
200
201
202 // see amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
203 pp_mg.add("usecg", 1);
204
205 // bottom_solver_eps, see amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
206 pp_mg.add("rtol_b", 0.0001);
207
208 // see amrex/Src/LinearSolvers/C_to_F_MG/AMReX_MGT_Solver.cpp
209 pp_mg.add("cg_solver", 1);
210
211 // maximum number of levels (max_nlevel)
212 pp_mg.add("numLevelsMAX", 1024);
213
214 /* MG_SMOOTHER_GS_RB = 1 (red-black Gauss-Seidel)
215 * MG_SMOOTHER_JACOBI = 2 (Jacobi)
216 * MG_SMOOTHER_MINION_CROSS = 5
217 * MG_SMOOTHER_MINION_FULL = 6
218 * MG_SMOOTHER_EFF_RB = 7
219 */
220 pp_mg.add("smoother", 1);
221
222 /* The type of multigrid to do
223 *
224 * MG_FCycle = 1 (F, full multigrid)
225 * MG_WCycle = 2 (W-cycles)
226 * MG_VCycle = 3 (V-cycles)
227 * MG_FVCycle = 4 (F + V cycles)
228 */
229 pp_mg.add("cycle_type", 1);
230
231 /* the type of bottom solver to use
232 *
233 * BiCG: biconjugate gradient stabilized (bottom_solver = 1)
234 * CG: conjugate gradient method (bottom_solver = 2)
235 * CABiCG: (bottom_solver = 3)
236 * special: bottom_solver = 4
237 *
238 * The special bottom solver extends the range of the multigrid coarsening
239 * by aggregating coarse grids on the original mesh together and further coarsening.
240 *
241 * if use_cg == 1 && cg_solver == 0 --> CG
242 * if use_cg == 1 && cg_solver == 1 --> BiCG
243 * if use_cg == 1 && cg_solver == 2 --> CABiCG
244 * if use_cg == 0 --> CABiCG
245 */
246 pp_mg.add("bottom_solver", 1);
247
248 // verbosity of the bottom solver. Higher numbers give more verbosity (doc: cg_verbose)
249 amrex::ParmParse pp_cg("cg");
250 pp_cg.add("v", 0);
251
252
253 /* Additional parameters that can't be set by ParmParse:
254 *
255 * - max_bottom_nlevel = 3 (additional coarsening if you use bottom_solver == 4)
256 * - min_width = 2 (minimum size of grid at coarsest multigrid level)
257 * - max_L0_growth = -1
258 */
259 amrex::MGT_Solver::def_min_width = 2;
260 amrex::MGT_Solver::def_max_L0_growth = -1.0;
261}
262
263
265 const AmrFieldContainer_pt& phi,
266 const amrex::Vector< AmrFieldContainer_pt > & grad_phi_edge,
267 const GeomContainer_t& geom,
268 int baseLevel,
269 int finestLevel)
270{
271 int nlevs = finestLevel - baseLevel + 1;
272
273 GeomContainer_t geom_p(nlevs);
274 AmrFieldContainer_pt rho_p(nlevs);
275 AmrFieldContainer_pt phi_p(nlevs);
276
277 for (int ilev = 0; ilev < nlevs; ++ilev) {
278 geom_p[ilev] = geom[ilev + baseLevel];
279 rho_p[ilev] = rho[ilev + baseLevel];
280 phi_p[ilev] = phi[ilev + baseLevel];
281 }
282
283 //FIXME Refinement ratio
284 amrex::IntVect crse_ratio = (baseLevel == 0) ?
285 amrex::IntVect::TheZeroVector() : itsAmrObject_mp->refRatio(0);
286
287 amrex::FMultiGrid fmg(geom_p, baseLevel, crse_ratio);
288
289 if (baseLevel == 0)
290 fmg.set_bc(bc_m, *phi_p[baseLevel]);
291 else
292 fmg.set_bc(bc_m, *phi_p[baseLevel-1], *phi_p[baseLevel]);
293
294 /* (alpha * a - beta * (del dot b grad)) phi = rho
295 * (b = 1)
296 *
297 * The function call set_const_gravity_coeffs() sets alpha = 0.0
298 * and beta = -1 (in MGT_Solver::set_const_gravity_coeffs)
299 *
300 * --> (del dot grad) phi = rho
301 */
302 fmg.set_const_gravity_coeffs();
303
304 // order of approximation at Dirichlet boundaries
305 fmg.set_maxorder(3);
306
307 int always_use_bnorm = 0;
308 int need_grad_phi = 1;
309
310 double residNorm = fmg.solve(phi_p, rho_p, reltol_m, abstol_m, always_use_bnorm, need_grad_phi);
311
312 for (int ilev = 0; ilev < nlevs; ++ilev) {
313 int amr_level = ilev + baseLevel;
314 fmg.get_fluxes(grad_phi_edge[amr_level], ilev);
315 }
316
317 return residNorm;
318}
319
320
321/*
322void FMGPoissonSolver::interpolate_m(AmrScalarFieldContainer_t& phi,
323 const GeomContainer_t& geom,
324 double l0norm,
325 int finestLevel)
326{
327 amrex::PhysBCFunct cphysbc, fphysbc;
328 int lo_bc[] = {INT_DIR, INT_DIR, INT_DIR}; // periodic boundaries
329 int hi_bc[] = {INT_DIR, INT_DIR, INT_DIR};
330 amrex::Vector<amrex::BCRec> bcs(1, amrex::BCRec(lo_bc, hi_bc));
331 amrex::PCInterp mapper;
332
333 std::unique_ptr<AmrField_t> tmp;
334
335 for (int lev = 1; lev <= finestLevel; ++lev) {
336 const AmrGrid_t& ba = phi[lev]->boxArray();
337 const AmrProcMap_t& dm = phi[lev]->DistributionMap();
338 tmp.reset(new AmrField_t(ba, dm, 1, 0));
339 tmp->setVal(0.0);
340
341 amrex::InterpFromCoarseLevel(*tmp, 0.0, *phi[lev-1],
342 0, 0, 1, geom[lev-1], geom[lev],
343 cphysbc, fphysbc,
344 itsAmrObject_mp->refRatio(lev-1),
345 &mapper, bcs);
346 phi[lev]->plus(*tmp, 0, 1, 0);
347 phi[lev-1]->mult(l0norm, 0, 1);
348 }
349
350 phi[finestLevel]->mult(l0norm, 0, 1);
351}
352*/
amr::AmrField_t AmrField_t
Definition: PBunchDefs.h:34
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
#define PAssert(c)
Definition: PAssert.h:102
amrex::Vector< AmrBoxLib::AmrField_t * > AmrFieldContainer_pt
double abstol_m
Absolute tolerance for solver.
Inform & print(Inform &os) const
double getZRangeMax(unsigned short level=0)
double reltol_m
Relative tolearance for solver.
double getYRangeMin(unsigned short level=0)
AmrBoxLib::AmrGeomContainer_t GeomContainer_t
void solve(AmrScalarFieldContainer_t &rho, AmrScalarFieldContainer_t &phi, AmrVectorFieldContainer_t &efield, unsigned short baseLevel, unsigned short finestLevel, bool prevAsGuess=true)
AmrBoxLib::AmrScalarFieldContainer_t AmrScalarFieldContainer_t
FMGPoissonSolver(AmrBoxLib *itsAmrObject_p)
AmrBoxLib::AmrVectorFieldContainer_t AmrVectorFieldContainer_t
AmrBoxLib::AmrProcMap_t AmrProcMap_t
double getXRangeMax(unsigned short level=0)
double solveWithF90_m(const AmrFieldContainer_pt &rho, const AmrFieldContainer_pt &phi, const amrex::Vector< AmrFieldContainer_pt > &grad_phi_edge, const GeomContainer_t &geom, int baseLevel, int finestLevel)
double getZRangeMin(unsigned short level=0)
double getXRangeMin(unsigned short level=0)
double getYRangeMax(unsigned short level=0)
int bc_m[2 *AMREX_SPACEDIM]
Boundary conditions.
AmrBoxLib::AmrGrid_t AmrGrid_t
The base class for all OPAL exceptions.
Definition: OpalException.h:28
Definition: Inform.h:42