29#include <AMReX_ParmParse.H>
30#include <AMReX_Interpolater.H>
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;
49 unsigned short baseLevel,
50 unsigned short finestLevel,
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;
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;
69 amrex::Vector< AmrScalarFieldContainer_t > grad_phi_edge(rho.size());
71 amrex::Vector< AmrField_t > tmp(rho.size());
73 PAssert(baseLevel <= finestLevel);
74 PAssert(finestLevel < (
int)rho.size());
76 for (
int lev = baseLevel; lev <= finestLevel ; ++lev) {
78 grad_phi_edge[lev].resize(AMREX_SPACEDIM);
79 tmp[lev].define(rho[lev]->boxArray(), dmap, AMREX_SPACEDIM, 1);
81 for (
int n = 0;
n < AMREX_SPACEDIM; ++
n) {
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);
88 for (
int i = 0; i <= finestLevel; ++i) {
89 phi[i]->setVal(0.0, 1);
90 tmp[i].setVal(0.0, 1);
94 amrex::GetVecOfPtrs(phi),
95 amrex::GetVecOfVecOfPtrs(grad_phi_edge),
101 std::stringstream ss;
102 ss <<
"Residual norm: " << std::setprecision(16) << residNorm
103 <<
" > " <<
abstol_m <<
" (absolute tolerance)";
105 "Multigrid solver did not converge. " + ss.str());
108 for (
int lev = baseLevel; lev <= finestLevel; ++lev) {
109 amrex::average_face_to_cellcenter(tmp[lev],
110 amrex::GetVecOfConstPtrs(grad_phi_edge[lev]),
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());
119 efield[lev][j]->mult(-1.0, 0, 1);
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'
159 <<
"* ******************************************************************** " <<
endl;
178 amrex::ParmParse pp_mg(
"mg");
181 pp_mg.add(
"maxiter", 200);
184 pp_mg.add(
"maxiter_b", 200);
187 pp_mg.add(
"nu_1", 2);
190 pp_mg.add(
"nu_2", 2);
193 pp_mg.add(
"nu_b", 0);
196 pp_mg.add(
"nu_f", 8);
203 pp_mg.add(
"usecg", 1);
206 pp_mg.add(
"rtol_b", 0.0001);
209 pp_mg.add(
"cg_solver", 1);
212 pp_mg.add(
"numLevelsMAX", 1024);
220 pp_mg.add(
"smoother", 1);
229 pp_mg.add(
"cycle_type", 1);
246 pp_mg.add(
"bottom_solver", 1);
249 amrex::ParmParse pp_cg(
"cg");
259 amrex::MGT_Solver::def_min_width = 2;
260 amrex::MGT_Solver::def_max_L0_growth = -1.0;
266 const amrex::Vector< AmrFieldContainer_pt > & grad_phi_edge,
271 int nlevs = finestLevel - baseLevel + 1;
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];
284 amrex::IntVect crse_ratio = (baseLevel == 0) ?
287 amrex::FMultiGrid fmg(geom_p, baseLevel, crse_ratio);
290 fmg.set_bc(
bc_m, *phi_p[baseLevel]);
292 fmg.set_bc(
bc_m, *phi_p[baseLevel-1], *phi_p[baseLevel]);
302 fmg.set_const_gravity_coeffs();
307 int always_use_bnorm = 0;
308 int need_grad_phi = 1;
310 double residNorm = fmg.solve(phi_p, rho_p,
reltol_m,
abstol_m, always_use_bnorm, need_grad_phi);
312 for (
int ilev = 0; ilev < nlevs; ++ilev) {
313 int amr_level = ilev + baseLevel;
314 fmg.get_fluxes(grad_phi_edge[amr_level], ilev);
amr::AmrField_t AmrField_t
Inform & endl(Inform &inf)
AmrBoxLib * itsAmrObject_mp
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.