# multigrid package¶

This is the pyro multigrid solver. THere are several versions.

MG implements a second-order discretization of a constant-coefficient Helmholtz equation:

(alpha - beta L) phi = f


variable_coeff_MG implements a variable-coefficient Poisson equation:

div { eta grad phi } = f


general_MG implements a more general elliptic equation:

alpha phi + div { beta grad phi } + gamma . grad phi = f


All use pure V-cycles to solve elliptic problems

## multigrid.MG module¶

The multigrid module provides a framework for solving elliptic problems. A multigrid object is just a list of grids, from the finest mesh down (by factors of two) to a single interior zone (each grid has the same number of guardcells).

The main multigrid class is setup to solve a constant-coefficient Helmholtz equation:

(alpha - beta L) phi = f


where L is the Laplacian and alpha and beta are constants. If alpha = 0 and beta = -1, then this is the Poisson equation.

We support Dirichlet or Neumann BCs, or a periodic domain.

The general usage is as follows:

a = multigrid.CellCenterMG2d(nx, ny, verbose=1, alpha=alpha, beta=beta)


this creates the multigrid object a, with a finest grid of nx by ny zones and the default boundary condition types. alpha and beta are the coefficients of the Helmholtz equation. Setting verbose = 1 causing debugging information to be output, so you can see the residual errors in each of the V-cycles.

Initialization is done as:

a.init_zeros()


this initializes the solution vector with zeros (this is not necessary if you just created the multigrid object, but it can be used to reset the solution between runs on the same object).

Next:

a.init_RHS(zeros((nx, ny), numpy.float64))


this initializes the RHS on the finest grid to 0 (Laplace’s equation). Any RHS can be set by passing through an array of (nx, ny) values here.

Then to solve, you just do:

a.solve(rtol = 1.e-10)


where rtol is the desired tolerance (residual norm / source norm)

to access the final solution, use the get_solution method:

v = a.get_solution()


For convenience, the grid information on the solution level is available as attributes to the class,

a.ilo, a.ihi, a.jlo, a.jhi are the indices bounding the interior of the solution array (i.e. excluding the ghost cells).

a.x and a.y are the coordinate arrays a.dx and a.dy are the grid spacings

class multigrid.MG.CellCenterMG2d(nx, ny, ng=1, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, xl_BC_type='dirichlet', xr_BC_type='dirichlet', yl_BC_type='dirichlet', yr_BC_type='dirichlet', xl_BC=None, xr_BC=None, yl_BC=None, yr_BC=None, alpha=0.0, beta=-1.0, nsmooth=10, nsmooth_bottom=50, verbose=0, aux_field=None, aux_bc=None, true_function=None, vis=0, vis_title='')[source]

Bases: object

The main multigrid class for cell-centered data.

We require that nx = ny be a power of 2 and dx = dy, for simplicity

get_solution(grid=None)[source]

Return the solution after doing the MG solve

If a grid object is passed in, then the solution is put on that grid – not the passed in grid must have the same dx and dy

Returns: out : ndarray
get_solution_gradient(grid=None)[source]

Return the gradient of the solution after doing the MG solve. The x- and y-components are returned in separate arrays.

If a grid object is passed in, then the gradient is computed on that grid. Note: the passed-in grid must have the same dx, dy

Returns: out : ndarray, ndarray
get_solution_object()[source]

Return the full solution data object at the finest resolution after doing the MG solve

Returns: out : CellCenterData2d object
grid_info(level, indent=0)[source]

Report simple grid information

init_RHS(data)[source]

Initialize the right hand side, f, of the Helmholtz equation (alpha - beta L) phi = f

Parameters: data : ndarray An array (of the same size as the finest MG level) with the values to initialize the solution to the elliptic problem.
init_solution(data)[source]

Initialize the solution to the elliptic problem by passing in a value for all defined zones

Parameters: data : ndarray An array (of the same size as the finest MG level) with the values to initialize the solution to the elliptic problem.
init_zeros()[source]

Set the initial solution to zero

smooth(level, nsmooth)[source]

Use red-black Gauss-Seidel iterations to smooth the solution at a given level. This is used at each stage of the V-cycle (up and down) in the MG solution, but it can also be called directly to solve the elliptic problem (although it will take many more iterations).

Parameters: level : int The level in the MG hierarchy to smooth the solution nsmooth : int The number of r-b Gauss-Seidel smoothing iterations to perform
solve(rtol=1e-11)[source]

The main driver for the multigrid solution of the Helmholtz equation. This controls the V-cycles, smoothing at each step of the way and uses simple smoothing at the coarsest level to perform the bottom solve.

Parameters: rtol : float The relative tolerance (residual norm / source norm) to solve to. Note that if the source norm is 0 (e.g. the righthand side of our equation is 0), then we just use the norm of the residual.
v_cycle(level)[source]

Perform a V-cycle for a single 2-level solve. This is applied recursively do V-cycle through the entire hierarchy.

## multigrid.edge_coeffs module¶

class multigrid.edge_coeffs.EdgeCoeffs(g, eta, empty=False)[source]

Bases: object

a simple container class to hold edge-centered coefficients and restrict them to coarse levels

restrict()[source]

restrict the edge values to a coarser grid. Return a new EdgeCoeffs object

## multigrid.general_MG module¶

This multigrid solver is build from multigrid/MG.py and implements a more general solver for an equation of the form:

alpha phi + div { beta grad phi } + gamma . grad phi = f


where alpha, beta, and gamma are defined on the same grid as phi. These should all come in as cell-centered quantities. The solver will put beta on edges. Note that gamma is a vector here, with x- and y-components.

A cell-centered discretization for phi is used throughout.

class multigrid.general_MG.GeneralMG2d(nx, ny, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, xl_BC_type='dirichlet', xr_BC_type='dirichlet', yl_BC_type='dirichlet', yr_BC_type='dirichlet', xl_BC=None, xr_BC=None, yl_BC=None, yr_BC=None, nsmooth=10, nsmooth_bottom=50, verbose=0, coeffs=None, true_function=None, vis=0, vis_title='')[source]

this is a multigrid solver that supports our general elliptic equation.

we need to accept a coefficient CellCenterData2d object with fields defined for alpha, beta, gamma_x, and gamma_y on the fine level.

We then restrict this data through the MG hierarchy (and average beta to the edges).

we need a new compute_residual() and smooth() function, that understands these coeffs.

smooth(level, nsmooth)[source]

Use red-black Gauss-Seidel iterations to smooth the solution at a given level. This is used at each stage of the V-cycle (up and down) in the MG solution, but it can also be called directly to solve the elliptic problem (although it will take many more iterations).

Parameters: level : int The level in the MG hierarchy to smooth the solution nsmooth : int The number of r-b Gauss-Seidel smoothing iterations to perform

## multigrid.mg_test_general_alphabeta_only module¶

Test the general MG solver with a variable coeffcient Helmholtz problem. This ensures we didn’t screw up the base functionality here.

Here we solve:

alpha phi + div . ( beta grad phi ) = f


with:

alpha = 1.0
beta = 2.0 + cos(2.0*pi*x)*cos(2.0*pi*y)

f = (-16.0*pi**2*cos(2*pi*x)*cos(2*pi*y) - 16.0*pi**2 + 1.0)*sin(2*pi*x)*sin(2*pi*y)


This has the exact solution:

phi = sin(2.0*pi*x)*sin(2.0*pi*y)


on [0,1] x [0,1]

We use Dirichlet BCs on phi. For beta, we do not have to impose the same BCs, since that may represent a different physical quantity. Here we take beta to have Neumann BCs. (Dirichlet BCs for beta will force it to 0 on the boundary, which is not correct here)

multigrid.mg_test_general_alphabeta_only.alpha(x, y)[source]
multigrid.mg_test_general_alphabeta_only.beta(x, y)[source]
multigrid.mg_test_general_alphabeta_only.f(x, y)[source]
multigrid.mg_test_general_alphabeta_only.gamma_x(x, y)[source]
multigrid.mg_test_general_alphabeta_only.gamma_y(x, y)[source]
multigrid.mg_test_general_alphabeta_only.test_general_poisson_dirichlet(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]

test the general MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark

multigrid.mg_test_general_alphabeta_only.true(x, y)[source]

## multigrid.mg_test_general_beta_only module¶

Test the general MG solver with a variable coeffcient Poisson problem (in essence, we are making this solver act like the variable_coefficient_MG.py solver). This ensures we didn’t screw up the base functionality here.

Here we solve:

div . ( beta grad phi ) = f


with:

beta = 2.0 + cos(2.0*pi*x)*cos(2.0*pi*y)

f = -16.0*pi**2*(cos(2*pi*x)*cos(2*pi*y) + 1)*sin(2*pi*x)*sin(2*pi*y)


This has the exact solution:

phi = sin(2.0*pi*x)*sin(2.0*pi*y)


on [0,1] x [0,1]

We use Dirichlet BCs on phi. For beta, we do not have to impose the same BCs, since that may represent a different physical quantity. Here we take beta to have Neumann BCs. (Dirichlet BCs for beta will force it to 0 on the boundary, which is not correct here)

multigrid.mg_test_general_beta_only.alpha(x, y)[source]
multigrid.mg_test_general_beta_only.beta(x, y)[source]
multigrid.mg_test_general_beta_only.f(x, y)[source]
multigrid.mg_test_general_beta_only.gamma_x(x, y)[source]
multigrid.mg_test_general_beta_only.gamma_y(x, y)[source]
multigrid.mg_test_general_beta_only.test_general_poisson_dirichlet(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]

test the general MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark

multigrid.mg_test_general_beta_only.true(x, y)[source]

## multigrid.mg_test_general_constant module¶

Test the general MG solver with a CONSTANT coefficient problem – the same one from the multigrid class test. This ensures we didn’t screw up the base functionality here.

We solve:

u_xx + u_yy = -2[(1-6x**2)y**2(1-y**2) + (1-6y**2)x**2(1-x**2)]
u = 0 on the boundary


this is the example from page 64 of the book A Multigrid Tutorial, 2nd Ed.

The analytic solution is u(x,y) = (x**2 - x**4)(y**4 - y**2)

multigrid.mg_test_general_constant.alpha(x, y)[source]
multigrid.mg_test_general_constant.beta(x, y)[source]
multigrid.mg_test_general_constant.f(x, y)[source]
multigrid.mg_test_general_constant.gamma_x(x, y)[source]
multigrid.mg_test_general_constant.gamma_y(x, y)[source]
multigrid.mg_test_general_constant.test_general_poisson_dirichlet(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]

test the general MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark

multigrid.mg_test_general_constant.true(x, y)[source]

## multigrid.mg_test_general_dirichlet module¶

Test the general MG solver with Dirichlet boundary conditions.

Here we solve:

alpha phi + div{beta grad phi} + gamma . grad phi = f


with:

alpha = 1.0
beta = cos(2*pi*x)*cos(2*pi*y) + 2.0
gamma_x = sin(2*pi*x)
gamma_y = sin(2*pi*y)

f = (-16.0*pi**2*cos(2*pi*x)*cos(2*pi*y) + 2.0*pi*cos(2*pi*x) +
2.0*pi*cos(2*pi*y) - 16.0*pi**2 + 1.0)*sin(2*pi*x)*sin(2*pi*y)


This has the exact solution:

phi = sin(2.0*pi*x)*sin(2.0*pi*y)


on [0,1] x [0,1]

We use Dirichlet BCs on phi.

For the coefficients we do not have to impose the same BCs, since that may represent a different physical quantity. beta is the one that really matters since it must be brought to the edges. Here we take beta to have Neumann BCs. (Dirichlet BCs for beta will force it to 0 on the boundary, which is not correct here)

multigrid.mg_test_general_dirichlet.alpha(x, y)[source]
multigrid.mg_test_general_dirichlet.beta(x, y)[source]
multigrid.mg_test_general_dirichlet.f(x, y)[source]
multigrid.mg_test_general_dirichlet.gamma_x(x, y)[source]
multigrid.mg_test_general_dirichlet.gamma_y(x, y)[source]
multigrid.mg_test_general_dirichlet.test_general_poisson_dirichlet(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]

test the general MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark

multigrid.mg_test_general_dirichlet.true(x, y)[source]

## multigrid.mg_test_general_inhomogeneous module¶

Test the general MG solver with inhomogeneous Dirichlet
boundary conditions.

Here we solve:

alpha phi + div{beta grad phi} + gamma . grad phi = f


with:

alpha = 10.0
beta = x*y + 1  (note: x*y alone doesn't work)
gamma_x = 1
gamma_y = 1

f =  -(pi/2)*(x + 1)*sin(pi*y/2)*cos(pi*x/2)
- (pi/2)*(y + 1)*sin(pi*x/2)*cos(pi*y/2) +
(-pi**2*(x*y+1)/2 + 10)*cos(pi*x/2)*cos(pi*y/2)


This has the exact solution:

phi = cos(pi*x/2)*cos(pi*y/2)


on [0,1] x [0,1], with Dirichlet boundary conditions:

phi(x=0) = cos(pi*y/2)
phi(x=1) = 0
phi(y=0) = cos(pi*x/2)
phi(y=1) = 0


For the coefficients we do not have to impose the same BCs, since that may represent a different physical quantity. beta is the one that really matters since it must be brought to the edges. Here we take beta to have Neumann BCs. (Dirichlet BCs for beta will force it to 0 on the boundary, which is not correct here)

multigrid.mg_test_general_inhomogeneous.alpha(x, y)[source]
multigrid.mg_test_general_inhomogeneous.beta(x, y)[source]
multigrid.mg_test_general_inhomogeneous.f(x, y)[source]
multigrid.mg_test_general_inhomogeneous.gamma_x(x, y)[source]
multigrid.mg_test_general_inhomogeneous.gamma_y(x, y)[source]
multigrid.mg_test_general_inhomogeneous.test_general_poisson_inhomogeneous(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]

test the general MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark

multigrid.mg_test_general_inhomogeneous.true(x, y)[source]
multigrid.mg_test_general_inhomogeneous.xl_func(y)[source]
multigrid.mg_test_general_inhomogeneous.yl_func(x)[source]

## multigrid.mg_test_simple module¶

an example of using the multigrid class to solve Laplace’s equation. Here, we solve:

u_xx + u_yy = -2[(1-6x**2)y**2(1-y**2) + (1-6y**2)x**2(1-x**2)]
u = 0 on the boundary


this is the example from page 64 of the book A Multigrid Tutorial, 2nd Ed.

The analytic solution is u(x,y) = (x**2 - x**4)(y**4 - y**2)

multigrid.mg_test_simple.f(x, y)[source]
multigrid.mg_test_simple.test_poisson_dirichlet(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]
multigrid.mg_test_simple.true(x, y)[source]

## multigrid.mg_test_vc_constant module¶

Test the variable coefficient MG solver with a CONSTANT coefficient problem – the same one from the multigrid class test. This ensures we didn’t screw up the base functionality here.

We solve:

u_xx + u_yy = -2[(1-6x**2)y**2(1-y**2) + (1-6y**2)x**2(1-x**2)]
u = 0 on the boundary


this is the example from page 64 of the book A Multigrid Tutorial, 2nd Ed.

The analytic solution is u(x,y) = (x**2 - x**4)(y**4 - y**2)

multigrid.mg_test_vc_constant.alpha(x, y)[source]
multigrid.mg_test_vc_constant.f(x, y)[source]
multigrid.mg_test_vc_constant.test_vc_constant(N)[source]
multigrid.mg_test_vc_constant.true(x, y)[source]

## multigrid.mg_test_vc_dirichlet module¶

Test the variable-coefficient MG solver with Dirichlet boundary conditions.

Here we solve:

div . ( alpha grad phi ) = f


with:

alpha = 2.0 + cos(2.0*pi*x)*cos(2.0*pi*y)

f = -16.0*pi**2*(cos(2*pi*x)*cos(2*pi*y) + 1)*sin(2*pi*x)*sin(2*pi*y)


This has the exact solution:

phi = sin(2.0*pi*x)*sin(2.0*pi*y)


on [0,1] x [0,1]

We use Dirichlet BCs on phi. For alpha, we do not have to impose the same BCs, since that may represent a different physical quantity. Here we take alpha to have Neumann BCs. (Dirichlet BCs for alpha will force it to 0 on the boundary, which is not correct here)

multigrid.mg_test_vc_dirichlet.alpha(x, y)[source]
multigrid.mg_test_vc_dirichlet.f(x, y)[source]
multigrid.mg_test_vc_dirichlet.test_vc_poisson_dirichlet(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]

test the variable-coefficient MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark

multigrid.mg_test_vc_dirichlet.true(x, y)[source]

## multigrid.mg_test_vc_periodic module¶

Test the variable-coefficient MG solver with periodic data.

Here we solve:

div . ( alpha grad phi ) = f


with:

alpha = 2.0 + cos(2.0*pi*x)*cos(2.0*pi*y)

f = -16.0*pi**2*(cos(2*pi*x)*cos(2*pi*y) + 1)*sin(2*pi*x)*sin(2*pi*y)


This has the exact solution:

phi = sin(2.0*pi*x)*sin(2.0*pi*y)


on [0,1] x [0,1]

We use Dirichlet BCs on phi. For alpha, we do not have to impose the same BCs, since that may represent a different physical quantity. Here we take alpha to have Neumann BCs. (Dirichlet BCs for alpha will force it to 0 on the boundary, which is not correct here)

multigrid.mg_test_vc_periodic.alpha(x, y)[source]
multigrid.mg_test_vc_periodic.f(x, y)[source]
multigrid.mg_test_vc_periodic.test_vc_poisson_periodic(N, store_bench=False, comp_bench=False, make_plot=False, verbose=1)[source]

test the variable-coefficient MG solver. The return value here is the error compared to the exact solution, UNLESS comp_bench=True, in which case the return value is the error compared to the stored benchmark

multigrid.mg_test_vc_periodic.true(x, y)[source]

## multigrid.mg_vis module¶

an example of using the multigrid class to solve Laplace’s equation. Here, we solve:

u_xx + u_yy = -2[(1-6x**2)y**2(1-y**2) + (1-6y**2)x**2(1-x**2)]
u = 0 on the boundary


this is the example from page 64 of the book A Multigrid Tutorial, 2nd Ed.

The analytic solution is u(x,y) = (x**2 - x**4)(y**4 - y**2)

multigrid.mg_vis.doit(nx, ny)[source]
multigrid.mg_vis.f(x, y)[source]
multigrid.mg_vis.true(x, y)[source]

## multigrid.project_periodic module¶

multigrid.project_periodic.doit(nx, ny)[source]

## multigrid.prolong_restrict_demo module¶

multigrid.prolong_restrict_demo.doit()[source]

## multigrid.variable_coeff_MG module¶

This multigrid solver is build from multigrid/MG.py and implements a variable coefficient solver for an equation of the form:

div { eta grad phi } = f


where eta is defined on the same grid as phi.

A cell-centered discretization is used throughout.

class multigrid.variable_coeff_MG.VarCoeffCCMG2d(nx, ny, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, xl_BC_type='dirichlet', xr_BC_type='dirichlet', yl_BC_type='dirichlet', yr_BC_type='dirichlet', nsmooth=10, nsmooth_bottom=50, verbose=0, coeffs=None, coeffs_bc=None, true_function=None, vis=0, vis_title='')[source]

this is a multigrid solver that supports variable coefficients

we need to accept a coefficient array, coeffs, defined at each level. We can do this at the fine level and restrict it down the MG grids once.

we need a new compute_residual() and smooth() function, that understands coeffs.

smooth(level, nsmooth)[source]

Use red-black Gauss-Seidel iterations to smooth the solution at a given level. This is used at each stage of the V-cycle (up and down) in the MG solution, but it can also be called directly to solve the elliptic problem (although it will take many more iterations).

Parameters: level : int The level in the MG hierarchy to smooth the solution nsmooth : int The number of r-b Gauss-Seidel smoothing iterations to perform